Planetary companions orbiting the M dwarfs GJ 724 and GJ 3988 A CARMENES and IRD collaboration ⋆

We report the discovery of two exoplanets around the M dwarfs GJ 724 and GJ 3988 using the radial velocity (RV) method. We obtained a total of 153 3.5m Calar Alto / CARMENES spectra for both targets and measured their RVs and activity indicators. We also added archival ESO / HARPS data for GJ 724 and infrared RV measurements from Subaru / IRD for GJ 3988. We searched for periodic and stable signals to subsequently construct Keplerian models, considering di ff erent numbers of planets, and we selected the best models based on their Bayesian evidence. Gaussian process (GP) regression was included in some models to account for activity signals. For both systems, the best model corresponds to one single planet. The minimum masses are 10 . 75 + 0 . 96 − 0 . 87 and 3 . 69 + 0 . 42 − 0 . 41 Earth-masses for GJ 724b and GJ 3988b, respectively. Both planets have short periods ( P < 10d) and, therefore, they orbit their star closely ( a < 0 . 05au). GJ 724b has an eccentric orbit ( e = 0 . 577 + 0 . 055 − 0 . 052 ), whereas the orbit of GJ 3988b is circular. The high eccentricity of GJ 724b makes it the most eccentric single exoplanet (to this date) around an M dwarf. Thus, we suggest a further analysis to understand its configuration in the context of planetary formation and architecture. In contrast, GJ 3988b is an example of a common type of planet around mid-M dwarfs.


Introduction
The number of detected exoplanets is growing thanks to the advances in astronomical instrumentation.High-precision spectrographs are capable of detecting radial velocity (RV) signals on the order of 1 m s −1 (e.g.HARPS, Mayor et al. 2003;CARMENES, Quirrenbach et al. 2014) and even lower (e.g.ESPRESSO, Pepe et al. 2021;MAROON-X,Seifahrt et al. 2018).In this sense, the detection of terrestrial exoplanets around low-mass stars, such as M dwarfs, is relatively easily attainable and, therefore, advantageous.This is due to the fact that these types of stars have lower masses (0.08-0.60 M ⊙ ) and are smaller (0.1-0.6 R ⊙ ) when compared to Sun-like stars, thus producing larger RV amplitude variations for a given planetary mass.
However, planetary signals can be mimicked by stellar activity (e.g.Robertson et al. 2015;Lubin et al. 2021;Gorrini et al. 2022), especially in M dwarfs, as they tend to be magnetically active (e.g.Reiners et al. 2012Reiners et al. , 2022;;Mignon et al. 2023).Therefore, it is necessary to treat this phenomenon when searching for exoplanets.As stellar activity can manifest itself in the stellar rotation period and its harmonics (Boisse et al. 2011), it is necessary to have a clear determination of the stellar rotation period and analyse the signals that appear on its integer fractions in depth.One successful method for treating stellar activity is ⋆ The RV data used in this work are only available in electronic form at the CDS via anonymous ftp to cdsarc.u-strasbg.fr(130.79.128.5) or via http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/ to model it with Gaussian process (GP) regression (e.g.Haywood et al. 2014;Rajpaul et al. 2015;Angus et al. 2018;Stock & Kemmer 2020;Kossakowski et al. 2022;Kemmer et al. 2022), but this must be used with caution, that is, by ensuring a proper (physical) training to avoid biased results, such as over-fitting or absorbing planetary signals (see e.g.Cabot et al. 2021;Rajpaul et al. 2021;Stock et al. 2023).
Another aspect that has to be taken into consideration when looking for planetary signals is a possible mean motion resonance (MMR) for systems with highly eccentric orbits, such as a 2:1 MMR can mimic an RV signal of a single high-eccentricity planet (e.g.Anglada-Escudé et al. 2010;Wittenmyer et al. 2013;Kürster et al. 2015;Nagel et al. 2019), leading to erroneous conclusions regarding the number of planets in the system.This is important when characterising planetary systems, as the planet's orbital solution can change drastically, coupled with our current knowledge of the incidence of eccentric and circular planetary orbits.
In this paper, we report planets orbiting two stars, GJ 724 (M1.0 V) and GJ 3988 (M4.5 V), observed with CARMENES as a part of the guaranteed time observations to search for exoplanets around M dwarfs (Ribas et al. 2023).GJ 3988 was also observed with the InfraRed Doppler (IRD) instrument at the Subaru telescope.In Sect.2, we summarise the stellar properties of the host stars, while in Sect.3, we describe our data set and measurements.We explain the methods used in the analysis in Sect. 4. In Sect.5, we determine the mean stellar rotation period for each star.We analyse all datasets in Sect.6 and discuss the implications of our study in Sect.7. Finally, we summarise our conclusions in Sect.8.

Stellar properties
In spite of the difference in spectral types, which translate into distinct stellar effective temperatures (T eff ), radii (R ⋆ ), and masses (M ⋆ ), the two planet-hosting stars have many parameters in common.While GJ 724 was already tabulated in the Bonner Durchmusterung (Schönfeld 1886), this star and GJ 3988 were (re-)discovered in high proper-motion surveys by Ross (1928) and Giclas et al. (1971), respectively.The farther distance of the former star partially compensates its earlier spectral type and, therefore, it is not very much brighter than the latter.However, their distances are so short (from 9.9 pc to 17 pc) that their classification as solar neighbours is assumed, supported by the wealth of public data available.
In Table 1, we provide a summary of the most important stellar properties for each star.First, we retrieved the spectral types, coordinates, magnitudes (G and J), parallactic distances, proper motions, and radial and rotation velocities from the literature for both GJ 724 and GJ 3988.There are numerous independent spectral type determinations that are consistent with the tabulated ones, between M0 and M1 for GJ 724 (e.g.Bidelman 1985;Gaidos et al. 2014) and between M4 and M5 for GJ 3988 (e.g.Lépine et al. 2013;Terrien et al. 2015).Something similar happens to the stellar atmospheric parameters (e.g.Passegger et al. 2018;Rajpurohit et al. 2018;Kuznetsov et al. 2019;Khata et al. 2021;Ishikawa et al. 2022).In our case, we used the values of T eff , log g, and [Fe/H] from CARMENES VIS and NIR template spectra determined by Marfil et al. 2021 with the SteParSyn code1 (Tabernero et al. 2022).The bolometric luminosities (L ⋆ ) were measured from the integration of the spectral energy distribution following Cifuentes et al. (2020) but with updated Gaia DR3 photometric and astrometric data (Gaia Collaboration et al. 2022).The radii and masses followed from L ⋆ , Stefan-Boltzmann's law, and the M ⋆ -R ⋆ relation by Schweitzer et al. (2019).We also computed galactocentric space velocities from the available data, which put both GJ 724 and GJ 3988 in the thin disc galactic population and measured mean rotation periods by fitting GPs to the photometric data as described in Sect. 5.The kinematics and long rotation periods are consistent with the pseudoequivalent width of the Hα line in faint absorption (Fuhrmeister et al. 2020), stringent upper limits to the rotational velocities imposed by the CARMENES spectral resolution (Reiners et al. 2018), no detection of X-ray flux (Voges et al. 1999;Stelzer et al. 2013), and (only for GJ 724) faintly measured magnetic fields and Ca ii emission (Astudillo-Defru et al. 2017;Perdelwitz et al. 2021;Reiners et al. 2022).The ages of both stars can be poorly constrained from these data, but are very likely of the order of 4 Gyr (Dungee et al. 2022, although gyrochonology at these late spectral types is not reliable; e.g Díez Alonso et al. 2019).Further details on stellar properties determination were provided by Caballero et al. (2022).
Finally, we address the potential multiplicity of GJ 724 and GJ 3988.Both stars were observed with the lucky imager Fast-Cam at the 1.5 m Telescopio Carlos Sánchez (Jódar et al. 2013;Cortés-Contreras et al. 2017).Besides, since GJ 3988 is located at d < 10 pc, it has also been the subject of numerous highresolution imaging surveys for companions at close and very close separations (down to ∼50 mas) with other lucky imagers (AstraLux at 2.2 m Calar Alto; Janson et al. 2014), coronographs in the near infrared (2.3 m Bok telescope at Steward Observatory; McCarthy & Zuckerman 2004), adaptive optics systems (Robo-AO; Lamman et al. 2020), and even the Hubble Space Telescope (with NICMOS and the F180M band;Dieterich et al. 2012).None of these surveys identified a close or very close companion candidate to GJ 724 and GJ 3988.Recently, Cifuentes (2023) looked for additional companions with Gaia DR3 data at even closer separations, with the RUWE parameter and other close-binarity astrophotometric indicators, and at wide and very wide separations of up to 2 × 10 5 au with a common proper-motion and parallax survey.Again, no companions were reported, so both GJ 724 and GJ 3988 appear to be single M dwarfs.

Spectroscopic data
The spectral data were obtained with CARMENES2 (Quirrenbach et al. 2014) for both GJ 724 and GJ 3988.We also incorporated RV data from the High Accuracy Radial velocity Planet Searcher spectrograph (HARPS; Mayor et al. 2003) for GJ 724.For GJ 3988, we used data from the InfraRed Doppler (IRD; Tamura et al. 2012;Kotani et al. 2018) at the Subaru Telescope.The main properties of both RV data sets are summarised in Table 2.

CARMENES
GJ 724 and GJ 3988 were observed with CARMENES, installed at the 3.5 m telescope at the Calar Alto observatory.It consists of two high-resolution spectrographs, covering the optical 5200-9600 Å (VIS) and near-infrared 9600-17100 Å (NIR) ranges, with spectral resolutions of R = 94 600 and R = 80 400, respectively.For GJ 724, we made use of 83 VIS observations over a time span of 2277 d (∼6 yr), from May 2016 to August 2022.As for GJ 3988, we used 70 VIS observations with a time span of 2313 d, from April 2016 to August 2022.
We did not include the NIR CARMENES RVs in our analysis as the errors and RMS were much higher (for GJ 724 by a factor of 4 and 2, respectively; whereas for GJ 3988 this was a factor of 4 and 3, respectively) than the data in the VIS channel given the same time sampling.Therefore, adding this data set would not contribute in any way to our models and would, in fact, decrease the quality of the fits.We refer to Bauer et al. (2020) for a comparative of VIS and NIR CARMENES RV data of a nearby M dwarf.
All raw spectra were processed with the CARMENES pipeline caracal (Caballero et al. 2016b).The RVs were computed using the serval code (Zechmeister et al. 2018), which uses a high signal-to-noise (S/N) stellar template to perform a least-square fit to each observation, including corrections of the instrumental drift from simultaneous measurements of the Fabry-Pérot interferometer and the Earth's barycentric velocity.As described by Tal-Or et al. (2019) and, especially, Ribas et al. (2023, see their Sect. 4.4), we also corrected RVs for nightly zero points.

HARPS
We also incorporated the data from observations of GJ 724 with HARPS, installed at the 3.6 m telescope located at La Silla observatory.It covers wavelengths between 3780 Å and 6910 Å and has a spectral resolution of 115 000.We made use of the RV data refined by Trifonov et al. (2020), which were re-processed using serval (Zechmeister et al. 2018), obtaining a better precision of the RV measurements when compared to the original ESO-DRS pipeline.The data were directly retrieved from their online catalog (Trifonov et al. 2020).The number of RV data points is 27, covering a time span of 1424 d (∼4 yr), from April 2008 to March 2012, and, therefore, with no overlap with the CARMENES data.

IRD
We observed GJ 3988 with the IRD instrument atop the Subaru 8.2 m telescope as part of the Subaru Strategic Program (SSP), which searches for exoplanets around mid-to-late M dwarfs by a blind Doppler survey.A total of 98 near-IR spectra (9500-17300 Å) were secured between February 2019 and July 2021.The integration time for each spectrum varied from 220 s to 1800 s, depending on the weather condition.We also injected light from the laser-frequency comb (LFC) into a reference fiber, which was used to trace the instantaneous instrumental profile of the spectrograph.
Wavelength-calibrated one-dimensional spectra for both stellar and LFC fibers were extracted with the standard IRD reduction pipeline (Kuzuhara et al. 2018).The typical S/N for the stellar spectra was 65-100 pix −1 at 10 000 Å. We measured RVs by fitting the individual spectra to a template spectrum of GJ 3988, which was also extracted from the observed IRD data, with the forward modeling technique (Hirano et al. 2020).The typical RV internal error was 2-3 m s −1 for each data point.
There are systematic offsets in the IRD RV measurements that are specific to each observing sequence (IRD observations are generally conducted on bright nights, so there are gaps between each sequence of bright nights).The offsets may be due to brightness variations of the LFC or other instrumental instabilities.We plotted the observing epochs versus the RV measurements of GJ 3988 in order to examine the offsets, to later correct them with the following post-processing.We grouped RV data points in each epoch with multiple targets being combined and computed the median of RV measurements in each group, allowing us to cancel signals from potential planets and stellar activity, leaving only instrumental signals.Finally, to apply the correction, we subtracted these median RV values from the individual RVs in each group.For this target, the central 80th percentile range of the offset values is −2 m s −1 to +2 m s −1 .
The observing strategy from IRD is to typically take multiple (short) exposures on a single night, with relatively low S/N values.Since we are not interested in short periods, but in getting observations comparable with CARMENES, we binned the data nightly.The binning resulted in a total of 50 IRD RV data points with a mean error and RMS of 2.81 m s −1 and 5.99 m s −1 , respectively.

Photometric data
We used photometric data from multiple sources for our analysis.An overview of the available data for each star can be found in Table 3. Outliers were removed by applying a 3σ clipping to all photometric datasets, which were afterwards binned daily.In the following, we describe the individual instruments and how the data were taken.

OSN
Photometric CCD observations for GJ 724 were collected at the Observatorio de Sierra Nevada (OSN) in Granada, Spain, with the T90 telescope, which is a 90 cm Ritchey-Chr'etien telescope equipped with a CCD camera VersArray 2k×2k with a field of view of 13.2×13.2arcmin 2 The camera is based on a high quantum efficiency back-illuminated CCD chip, type Marconi-EEV CCD42-4, with optimised response in the ultraviolet (Amado et al. 2021).The observations were collected in both Johnson V and R filters on 57 nights, during two runs: run 1 consisted of 24 epochs obtained during the period August-October 2020, while run 2 consisted of 33 epochs obtained during the period June-August 2021.Each epoch typically consisted of 20 exposures in each filter per night, of 30 s and 20 s in the V and R filters, respectively.
All CCD measurements were obtained by the method of synthetic aperture photometry using no binning.Each CCD frame was corrected in a standard way for bias and flat fielding.Different aperture sizes were also tested in order to choose the best one for our observations.A number of nearby and relatively bright stars within the frames were selected as check stars in order to choose the best ones to be used as reference stars.

TJO
We include photometric observations of GJ 724 and GJ 3988 from the 0.8 m Telescopi Joan Oró (TJO) at the Parc Astrònomic del Montsec, Lleida, Spain.We made use of the Johnson R filter of the LAIA imager, a 4k×4k CCD with a field of view of 30×30 arcmin 2 and a scale of 0.4 arcsec pixel −1 .The images were calibrated with darks, bias and flat fields with the ICAT pipeline (Colome & Ribas 2006) of the TJO.The differential photometry was extracted with AstroImageJ (Collins et al. 2017) using the aperture size that minimised the root mean square (RMS) of the resulting relative fluxes, and a selection of the 30 brightest comparison stars in the field that did not show variability.Then, we used our own pipelines to remove outliers and measurements affected by poor observing conditions or presenting a low S/N.

MEarth
For GJ 3988, we included publicly available data from the MEarth Project (MEarth; Berta et al. 2012;Irwin et al. 2015), corresponding to the Data Release 11 (DR11).MEarth has two different facilities: MEarth-North, located at the Fred Lawrence Whipple Observatory in Arizona, USA, and MEarth-South, located at Cerro Tololo Inter-American Obervatory in Chile.We made use of data from MEarth-North, which consists of eight near-identical telescope systems, with 0.4 m f /9 Ritchey-Chrétien Cassegrain design on a German equatorial mount.The pixel scale is approximately 0.76 arcsec pixel −1 .The filter is fixed and consists of 2×1.5 mm plates of RG715 Schott glass built into the camera housing itself.
Because of possible offsets and other unknown effects, we divided the MEarth light curves into the individual instrument sections denoted in the light curve files.Specifically, we used telescope 7 (version 6, 7, and 9 to 12), and telescope 8 (versions 5, 8, 10), as the other data sets contained only a few data points (< 15) or had short baselines (<50 d) and, therefore, the inclusion of these few points may result in higher uncertainties due to unknown offsets between the instrument versions.

SuperWASP
GJ 3988 was also monitored by the Super-Wide Angle Search for Planets (SuperWASP; Pollacco et al. 2006).The survey comprises two identical robotic telescopes, at La Palma, Spain, and Sutherland, South Africa, each with eight cameras observing through broadband filters (4000-7000 Å).We used data spanning April 2006 to August 2008 from the first public data release (Butters et al. 2010), totalling 260 daily-binned epochs over three seasons.The light curves were extracted via aperture photometry by the SuperWASP team, and corrected for instrumental systematics using the method of Tamuz et al. 2005.

Periodogram analysis
To search for periodic signals in the available RVs, activity indicators, and photometry, we used the generalised Lomb-Scargle periodogram (GLS; Zechmeister & Kürster 2009) implemented in exostriker (Trifonov 2019) and astropy (Astropy Collaboration et al. 2013Collaboration et al. , 2018Collaboration et al. , 2022)).The periodograms were normalised by the residuals of the data following Zechmeister & Kürster (2009) and the false-alarm probabilities (FAPs) were calculated using the analytic expression from Baluev (2008).Our plots show the FAP levels of 10 %, 1 %, and 0.1 %; we consider a signal to be significant if it has a FAP lower than 1 %.The temporal stability of signals was assessed using the stacked Bayesian GLS peridogram (s-BGLS; Mortier et al. 2015;Mortier & Collier Cameron 2017).

Modelling
For all modelling tasks, we used the juliet package (Espinoza et al. 2019).It combines the publicly available packages radvel (Fulton et al. 2018) for fitting RVs and batman (Kreidberg 2015) for fits of photometry.For the sampling within juliet, we chose the dynamic nested sampling algorithm (Higson et al. 2019) from dynesty (Speagle 2020; Koposov et al. 2022), which allowed us to determine the Bayesian evidence for our models.We assumed that a model performs significantly better if the difference in ln Z compared to the competing model was greater than 5. Models with ∆ ln Z > 2.5 were considered to be moderately preferred, and all models below that threshold were considered indistinguishable (following e.g.Trotta 2008).In addition to the static Keplerian model components, juliet provides GP as a red noise component using the celerite (Foreman-Mackey et al. 2017) and george (Ambikasaran et al. 2015) packages.As GP models can flexibly account for correlated noise, we compared the models with a GP separately from those without a GP (only Keplerian).
We used three different kernels to model stellar activity present in our data on the basis of a case-by-case examination: the double simple harmonic oscillator kernel (dSHO; David et al. 2019;Gillen et al. 2020), as described by Kossakowski et al. (2021); the quasi-periodic kernel (QP; Haywood et al. 2014;Rajpaul et al. 2015), following the parametrisation outlined by Espinoza et al. (2019), and lastly the quasi-periodic cosine kernel (QPC; Perger et al. 2021), which is the sum of the QP kernel and a cosine function.The latter is parameterised by the same hyperparameters as the QP kernel but also including the amplitudes of the QP component and the cosine function.An overview of the parameters and standard priors used for the GP modelling is given in Table A.1.

Spectroscopic rotation period
In Fig. 1 we show the periodograms of the photometric data and the activity indicators of GJ 724 that reach at least a FAP level below 10 %.The most significant signal (FAP <0.1 %) from the activity periodograms corresponds to ∼180 d in the TiO 8858Å indicator, which is also present in the CRX, CaH3, Ca IRT c, TiO 7048 Å, Na D 1 , and WFB indices (albeit not at a significant level).This signal is the second harmonic of the period of one year and is therefore likely to be caused by the contamination of the spectra.
Other signals with a FAP below 0.1 % are found around 122 d and 56 d in the TiO 8858Å index.The former signal is also a harmonic of a year (P 1/3 ), whereas the 56 d period ap-pears as well in other indicies (Na D 1 , Ca IRT b and Ca IRT c, TiO 7050, He λ10830 Å, and WFB) with a FAP level between 10 % and 1 %.Half of this periodicity, that is 28 d, is also present in other activity indicators (dLW, Ca IRT a, Ca IRT b, Ca IRT c, H α , and He λ10830 Å), being only significant in the Ca IRT b and Ca IRT c indices.While this signal is close to the moon cycle, it matches also the second harmonic of the 56 d signal.
There are also many other significant peaks between 20 d and 40 d (dLW, Ca IRT b, Ca IRT c, and He 10830Å).In order to check if they originate from a common period related to stellar activity, we performed a GP fit to the activity indicators (separately) using the dSHO kernel with uniform priors on the period from 2 d to 100 d.We only get a clear detection of the Ca IRT c index, with a period around ∼50 d.The Ca IRT a and Ca IRT b indices resulted in bimodal distributions of the period with values of ∼28 d and ∼56 d.By constraining the period of the TiO λ8858Å index to 40 d to 60 d, where the most significant signals are for this indicator (excluding the one year harmonic), we obtained a period of around 57 d.In this sense, we cannot retrieve an accurate rotation period from the activity indicators.

Photometric rotation period
The GLS periodogram of the (combined) OSN and TJO photometric data is also displayed in Fig. 1.The strongest peak is at ∼28 d with a FAP of 2.9 × 10 −6 %, followed by two slightly weaker signals (FAP > 1 × 10 −5 %) at around 147 d and 32 d.We performed a GP on the photometric data in order to retrieve the stellar rotation period of GJ 724 using the dSHO kernel.All the GP parameters were shared between the different datasets, except for the amplitude.We searched for periodicities between 2 and 200 days, using a uniform prior on the periodic component.There is a detection on the posterior distribution of the GP rotation period parameter of 28.4 ± 0.1 d, but for the parameters Q0 and dQ there is a broader distribution for periods below 10 d.By constraining both of these parameters from 10 to 10 5 d, the GP parameter of the rotation period resulted in 57.0 ± 1.0 d.
Thus far, both the spectroscopic and photometric determinations of the stellar rotation period have shown bimodal results (∼ 28 d and ∼ 56 d), so we decided to implement the other two kernels (see Sect. 4.2).The QPC kernel detects a rotation period of 60.0 ± 3.0 d, while the QP kernel prefers a period of 29.8 ± 1.1 d.This latter result is consistent as this kernel is missing the second component at half of the rotation period (P rot /2), which is a mode of the periodic variability that often appears for activity due to the fact that any spot distribution can be to first order approximated by two spots on the star (Jeffers & Keller 2009).
In Fig. 2, we show the GP hyper-parameters (dQ and Q) of the dSHO kernel against the periodic component, indicating that the likelihood and number of posterior samples favour a rotation period of ∼57 d instead of its second harmonic.It is therefore reasonable to adopt the value of the stellar rotation period from the GP model with the (constrained) dSHO kernel: 57.0 ± 1.0 d, but we cannot confidently rule out the ∼ 28 d as the actual rotation period.We note that the uncertainties we obtain are small, for which we interpret our result as the mean rotation period over all latitudes.Therefore, the calculated errors correspond to the errors in the mean and not to the uncertainties in the actual stellar rotation period.

Spectroscopic rotation period
There are only a few significant signals in the GLS periodograms of the activity indicators obtained for GJ 3988 (see Fig. 3).The strongest signal is apparent in the TiO 7048Å band index with a period of ∼114 d and a FAP < 0.1 %.Equivalently, the TiO 9960 Å and TiO 8428 Å indices show strong signals at ∼95 d and ∼90 d, which can be related to the ∼114 d signal through aliasing caused by a roughly yearly sampling frequency apparent in the window function of the CARMENES observations.However, both are also close to the fourth harmonic of one year  Since the signals of the TiO indices are the strongest in the GLS periodograms, we used a GP fit in the next step to determine a more precise estimate for the stellar rotation period from them.
All three indices originate from the same molecule, which is why we performed a joint fit in which we shared the GP's rotation period between the indices.In doing so, we set the prior for the period uniform between 50 d and 150 d to cover the range of the signals visible in the GLS.The result is a spectroscopic stellar rotation period of 120 ± 10 d.

Photometric rotation period
We used the available MEarth, SuperWasp and TJO photometry to determine the photometric rotation period of GJ 3988 as described in Section Sect.5.1.2for GJ 724.A first fit was implemented using the dSHO kernel with a shared rotation period varying freely between 10 d and 200 d.The results yielded to a split posterior distribution for the rotation period, indicating an unstable signal with a low-quality factor and long-term periodicity with a period larger than ∼165 d (peaked at ∼175 d), along with a more stable signal at ∼115 d.Since the 115-day signal matches the period determined from the TiO band indices, we performed a second fit with constrained priors on the difference in quality factor (10 -1 × 10 10 ), from which we determined a photometric stellar rotation period of 116 ± 3 d.As this determination is robust and straightforward, we did not incorporate other kernels like in the case of GJ 724.Moreover, we also interpret the small errors as the uncertainties over the mean value of the rotation period of the star.

RV signal search
The most significant signal present in the GLS periodogram of the RVs is at 5.1 d (FAP = 6.0 × 10 −7 %), as seen in Fig. 4.There are also two other significant signals (FAP of the order of 10 −5 % -10 −4 %) with short periods (P = 0.83 d and P = 1.24 d), due to the 1 d alias.Using the AliasFinder (Stock & Kemmer 2020;Stock et al. 2020) we found 5.1 d to be most likely the underlying true period.The 5.1 d signal is not significant in any periodogram of the activity indicators and it is stable over time (as displayed in Fig. 5), unlike the mean stellar rotation period at 57 d or its harmonic at 28 d.There are no significant signal left in the periodogram after subtracting the 5.1 d signal with a single Keplerian model, but the eccentricity results in a high value (e = 0.577 +0.055 −0.052 ).Consequently, we ran a circular model, but the evidence clearly favours the eccentric one (∆ ln Z = −21.6).

RV modelling
As seen in the GLS periodogram (Fig. 4), there are no other significant peaks besides the 5.1 d signal (and its aliases at P < 2 d).
Since the single Keplerian model results in a high eccentricity, we checked different configurations.As two planets in 2:1 MMR can be misinterpreted as a single high-eccentricity planet (Anglada-Escudé et al. 2010), we also tested a circular model with two planets in a 2:1 mean motion resonance.To prevent conflict with the priors, we use uniform priors on the period between 4 d and 6 d for the outer planet and between 1.5 d to 3.5 d for the inner planet.As stellar activity can affect RV data even though it does not appear significant in the GLS (e.g.Barnes et al. 2011;Boisse et al. 2011;Haywood et al. 2014), we also included a GP component in some models in order to test this scenario.We use normal priors on the periodic component cen-  tred at the mean rotation period or at its second harmonic.The priors of the fits are displayed in Table A.2. Table 4 displays the comparison between the models based on their Bayesian evidence.For the Keplerian only fits (without a GP), the model with the best evidence corresponds to one planet in an eccentric orbit, as ∆ ln Z > 5 when compared to the other models.For the models including Keplerians and a GP, the model with the best evidence is one eccentric planet with a GP fitted at half of the mean rotation period.This model is highly favoured over the same model but with the GP tuned to the mean rotation period (∆ ln Z > 5).As shown in Fig. 1, many activity indicators displayed significant signals at 28 d (Sect.5.1.1),for which this result is reasonable.Since the best model comparison shows that the model including a GP with a period of 56 d is not as good as the model with the GP with a period of 28 d, we adhere to the latter.It is worth mentioning that the planetary parameters do not change between both models, as seen in All the models favour the 5.1 d signal.As this is a significant signal, and there is no other counterpart in the activity indica-tors and photometric time series, we conclude that it is due to a planet.We adopt the model with a GP with the highest evidence, that is the 1P ( 5

RV signal search
The GLS periodogram of the CARMENES and IRD data for GJ 3988 (presented in Fig. 7) shows a highly significant signal with a period of 6.94 d (FAP < 1 × 10 −6 %) accompanied by slightly less significant aliases, due to the dominant daily sampling frequency, at 1.16 d and 0.87 d.If the signal is subtracted from the data using a circular Keplerian model, the highest peak in the GLS periodogram is at 115 d, which is consistent with our measurement of GJ 3988's mean rotation period from Sect.5.2.The s-BGLS of the 6.94-day signal, as a measure of its stability over time, is shown in the left panel of Fig. 8.We do not detect any variability in the signal but only a steady increase in probability, which suggests a static nature for it.On the other hand, the s-BGLS focused on 115 d (right panel) shows slight variations in the probability over time in combination with shifts of the period, which indicates a signal of only quasi-periodic nature, as we would expect for signals related to stellar activity.

RV modelling
In the previous section, we detected two significant signals in the RVs of GJ 3988: a planet candidate signal with a period of 6.94 d and an imprint of the mean stellar rotation period at ∼116 d.We consequently compared a model considering pure stellar activity against models that additionally included the planet candidate.In doing so, we considered both, circular and eccentric orbits.The priors of the fits are listed in Table A.3.The results of the model comparison using the Bayesian evidence are presented in Table 6.We found that the models including the planet candidate are highly significant compared to the activity-only model.Based on this, in combination with the results from the analysis of the activity indicators and the s-BGLS, we qualified the signal as a genuine planetary signal.For the 'only Keplerian' models, the best one is the two-Keplerian in circular orbit, which corresponds to the planet plus a static sinusoidal component used to model the stellar activity (P rot =116 d).Moreover, for the models with Keplerian plus GP, the model with one planet in a circular orbit plus a GP is strongly preferred (∆ln Z > 5.0).Assuming that the circular model for the planet in combination with the quasi-periodic GP for the stellar activity are better physically motivated compared to two sinusoidal components, we chose the circular 1P (7 d-circ) + dSHO-GP (116 d) as our best fit model.
The posteriors from this fit are listed in Table 5 together with the derived parameters for the planet.In Fig. B.4 and Fig. 9, we show the RVs together with the best-fit model and the phasefolded plot of the planet, respectively.

Transit search and detection limits
In this section, we aim to confirm or refute the transiting nature of the planets GJ 724 b and GJ 3988 b, by exploring the public data provided by the TESS mission (Ricker et al. 2015).On the one hand, GJ 724 b has not been observed yet, but will be in sector 80 during the second TESS extended mission.Moreover, the current ephemeris for this planet provides observational windows with an uncertainty of about 6 h, which makes any time-critical observation using ground-based facilities very challenging.Then, unfortunately, we can not reach any conclusion concerning its transiting nature.On the other hand, GJ 3988 b has been observed in multiple sectors and cadences: during the primary mission in sectors 24, 25, and 26 using the 2 min and 30 min cadences, and during the extended mission in sectors 50 and 51 using the 20 sec, 2 min and 10 min cadences, and in sectors 52 and 53 using the 20 sec and the 2 min.For our study, we used all the data corresponding to the 2 min cadence, which cover seven sectors in total.
Neither Science Processing Operations Center (SPOC; Jenkins et al. 2016) nor the Quick-Look Pipeline (QLP; Huang et al. 2020) have issued an alert for GJ 3988, which is a hint of the nontransiting nature of this system.However, these pipelines may not detect some periodic transits if they are shallow and their S/N is below their detection thresholds.Then, the community is often conducting complementary planetary searches either using ground-based facilities (see, e.g., Delrez et al. 2022) or using alternative custom detection pipelines.In this context, we explored the TESS data using the SHERLOCK pipeline 4 , presented initially by Pozuelos et al. (2020) and Demory et al. (2020), and used in several studies (see, e.g., Wells et al. 2021 Parameter (a)  GJ 724 GJ 3988  σ IRD (m s −1 ) . . .3.17 +0.77   −0.72   Notes. (a) Error bars denote the 68% posterior credibility intervals. (b) These derived parameters were computed using the stellar parameters that were sampled from a normal distribution centred around the values from Table 1 with the width of the uncertainty. (c) Assuming a zero Bond albedo. (d) Computed using Eq. ( 4) from Quirrenbach (2022Quirrenbach ( ). et al. 2021;;Schanche et al. 2022).This pipeline allows the user to explore TESS data to recover known planets and alerts and to search for new periodic signals, which may hint at the existence of extra-transiting planets.In short, the pipeline combines six modules to (1) download and prepare the light curves from the MAST, ( 2) search for planetary candidates, (3) perform a semiautomatic vetting of the interesting signals, ( 4) compute a statistical validation, (5) model the signals to refine their ephemerides, and ( 6) compute observational windows from ground-based observatories to trigger a follow-up campaign.We refer the reader to Delrez et al. (2022) and Pozuelos et al. (2023) for recent SHERLOCK applications and further details.
During our search, we explored an orbital range from 0.5 to 40 d, employing a density grid of 247260 periods distributed in a logarithmic distribution, with shorter time steps for shorter orbital periods and larger time steps for longer orbital periods.The optimum period grid to explore was computed based on the stellar mass and radius and the time span covered by the data set (Hippke & Heller 2019).In our case, the minimum and maximum time steps were 4.65 × 10 −6 and 0.0016 d, respectively.A minimum of two transits was required to claim a detection.Fig. 8. Same as Fig. 5 but for GJ 3988.In the instrument coloured bar, CARMENES is shown in orange-red, while IRD in bluish-green.Left panel is centred on the 6.94-day signal, after subtracting the stellar activity using a GP model.Right panel is centred on the imprint of the mean stellar rotation period after subtracting the 6.94-day signal from the data.
We did not find any hint of a periodic signal which might correspond to GJ 3988 b or any other transiting planet in the system.All the signals detected were attributable to systematics or noise.Then, in this scenario, we wondered if the photometric precision provided by TESS is enough to detect the presence of GJ 3988 b if we assume that it is a transiting planet.To answer this question, we conducted an injection-and-recovery experiment using the MATRIX5 code (Dévora-Pajares & Pozuelos 2022).MATRIX injects synthetic planets over the 2 min cadence light curves corresponding to the seven sectors used in this study.
If we assume that GJ 3988 b is a transiting planet, employing the mass-radius relationship given by Chen & Kipping (2017), its radius would be 1.72 +0.5  −0.7 R ⊕ .Then, we explored the R planet -P planet parameter space in the ranges of 0.5-3.0R ⊕ with steps of 0.13 R ⊕ , and 0.5-8.0days with steps of 0.20 days.Moreover, for each combination of R planet -P planet MATRIX explores three different phases, that is, different values of T 0 .In total, we explored 1350 scenarios.Once the synthetic planets are injected, MATRIX detrends the light curves using a bi-weight filter with a window size of 0.5 d, which was the optimal value during the SHERLOCK search.MATRIX considers a synthetic planet as recovered when its epoch matches the injected epoch with 1 hour accuracy, and its period is within 5 % of the injected period.
It is worth noting that for simplicity, the injected planets have impact parameters and eccentricities equal to zero and that since   we injected the synthetic signals in the PDCSAP light curve, these signals were not affected by the PDCSAP systematic corrections; hence, the detection limits that we find can be considered as the most optimistic scenario (see e.g.Pozuelos et al. 2020;Eisner et al. 2020).
The detectability map resulting from this injection-andrecovery experiment is shown in Fig. 10.We found that any planet within the range of the orbital period and radius explored in this experiment would be easily detected, yielding a recovery rate of nearly 100%.Hence, we can robustly confirm the non-transiting nature of GJ 3988 b, and the non-existence of any other transiting planet with an orbital period shorter than 8 d and radius larger than 0.5 R ⊕ .

GJ 724
GJ 724 b has a very high eccentricity (e ∼ 0.6).We tested the orbital configuration of two planets in a circular 2:1 MMR, with the eccentric model being the preferred one (∆ ln Z > 5).In Fig. 11 we plot GJ 724 b alongside all confirmed exoplanets around sin-Fig. 10.Detectability map of GJ 3988 b using the TESS 2 min data corresponding to sectors 24, 25, 26, 50, 51, 52, and 53.We explored a total of 1350 different scenarios.Larger recovery rates are presented in yellow and green colours, while lower recovery rates are in blue and darker hues.In the parameters space explored, each combination of R planet -P planet yielded a recovery rate of nearly 100%, which allows us to confirm the non-transiting nature of this system robustly.The blue dot refers to the planet GJ 3988 b assuming its transiting nature.By analysing the periodiogram of the RVs residuals of GJ 724, we do not find any signal that could be attributed to a second planet.If there is another companion, it is more likely that it has a longer orbital period, perturbing the orbit of the innermost planet by gravitational interaction and therefore explaining the high eccentricity of GJ 724 b.With the aim to check if there is a possible planet that we are not able to detect, we computed a detection limit map using the same method as in Sabotta et al. (2021) (displayed in Fig. 12).GJ 724 b lies in the region with 100% detection probability.A possible companion with a larger minimum mass could have been detected up to an orbital period of ∼100 d.Any potential disturber must have a lower minimum planetary mass than GJ 724 b in order to escape detection with our dataset.The population of synthetic planets orbiting a 0.5 M ⊙ star, as studied by Burn et al. (2021), can be used to investigate the prevalence of such systems as predicted by conventional planet formation models (Emsenhuber et al. 2021).We proceed by randomly drawing multiple inclinations for each system and applying the detection sensitivity for GJ 724 shown in Fig. 12.In this manner, synthetic planet detections can be generated with repeated inclusion of individual synthetic planets for improved statistical significance.This approach is described in detail in Schlecker et al. (2022).In Fig. 13 we show the resulting synthetic population of planets and their eccentricities resulting from theoretical N-body integration over 20 Myr and individual planet evolution including mass-loss and tidal in-spiral of 5 Gyr.It was applied to the synthetic planets with assigned inclinations drawn from an isotropically uniform distribution.We can see that the model generally fails to reproduce GJ 724 b-like planets.However, those eccentricity values do not include long-term tidal damping, which is discussed later in this paper.
Despite the failure to match the exact properties of GJ 724 b, it is insightful to analyse systems with eccentric planets.We selected four systems containing planets closest to GJ 724 b in mass, period, and eccentricity space and show the full evolution of the systems in Fig. 14.Two of the systems (Sim 445 and 897) have a tightly packed group of similar-mass planets, which can be ruled out by observations.Interestingly, the other two simulations (Sim 778 and 153), which only have one observable planet, show a very similar formation history.In both systems, while the disk was present, there were two planets with ∼10 M ⊕ locked in a mean-motion resonance close to the inner edge of the disk.However, over Gyr timescales, the innermost planet spirals into the star due to tidal interactions leaving only the outer planet which experiences weaker tides due to its lower mass, more compact structure, and larger orbital period.This is a promising pathway for systems like the one around GJ 724; however, more research into the exact dynamical evolution of this configuration is required.
Alternatively, close-in planets in eccentric orbits could undergo high-eccentricity tidal migration and circularisation (e.g.Rasio & Ford 1996;Rice et al. 2012;Giacalone et al. 2017;Dong et al. 2021;Yu et al. 2021Yu et al. , 2022)), a process not included in the aforementioned models (Emsenhuber et al. 2021).Similarly to the synthetic planets discussed above, eccentric inner planets can form during the early stages of the planetary system evolution by interacting with massive counterparts (e.g.Weidenschilling & Marzari 1996;Chatterjee et al. 2008;Jurić & Tremaine 2008;Raymond et al. 2009;Carrera et al. 2019).In order to check these two possibilities, we computed an estimate of the circularisation timescale.We did this by following Eqs. 1 and 2 of Jackson et al. (2009).
Since we do not know the radius of the planet, we have to rely on an estimate based on its characteristics.We measured only the minimum planetary mass of GJ 724 b, which is 10.75 +0.96 −0.87 M ⊕ , for which we can assume that it has around Neptune's mass.Earth radii.The stellar age is not known, but with a nearly solar metallicity, we can assume roughly solar age or 5 Gyr.The circularisation timescale is therefore shorter than the age of the system, which means that to explain our observation the planet would have to have suffered from a recent perturbation event, since the planet would remain within the 1 σ uncertainty interval of the determined eccentricity for only about 0.5 Gyr.Alternatively, the eccentricity and semi-major axis can be evolved back in time by changing the sign of Eqs. 1 and 2 from Jackson et al. (2009).With the same assumptions, 5 Gyr ago the planet would have had an eccentricity of 0.9, which means that it would have suffered from a scattering event at an early stage of the planetary system.This seems to be a more plausible explanation.Unfortunately, due to the uncertainty in the Q values, we cannot give a more precise result as we only have an estimate of from the solar system.
We also considered the possibility of a Kozai effect in the system induced by another companion, as in the case of GJ 436 b (Bourrier et al. 2018).The perturber could be another planet, a brown dwarf, or, under certain circumstances, a very low-mass star.As mentioned in Sect.2, GJ 724 does not seem to be part of a binary or multi-stellar system.Whereas for the planetary companion, we already discarded the possibility of a Jupiter-like planet in a Jupiter-like orbit with the detection map (Fig. 12).Nevertheless, a possible Kozai effect cannot be completely ruled out as an explanation for the high eccentricity of GJ 724 b.In this sense, more data is desirable to search for possible companions that had not been detected with our dataset.Future studies based on dynamical interactions could also check if such an (undetected) companion can cause such high eccentricity, for which we strongly suggest future works on this system.

GJ 3988
The orbital fit of GJ 3988 system is unambiguous as it is consistent with a circular orbit.Based on the planet's minimum mass (M b sin i = 3.69 +0.42  −0.41 M ⊕ ), we expect GJ 3988 b to be a super-Earth or mini-Neptune, as we do not know its radius and bulk density.Due to its proximity to the star (a = 0.0405 +0.0011 −0.0012 au), this planet (if it indeed Neptune-like) is likely to undergo atmospheric mass loss due to the emission received from its host star, or to lose its atmosphere if is a rocky planet (e.g.Lammer et al. 2003;Lecavelier des Etangs et al. 2004;Lammer et al. 2009;Erkaev et al. 2016;Owen & Wu 2016;Lehmer & Catling 2017;Locci et al. 2019).
There are other exoplanets with characteristics (orbital period and minimum mass) similar to those of GJ 3899 b.For instance, G 264-012 c (Amado et al. 2021) is a terrestrial planet (M c sin i = 3.75 +0.48  −0.47 M ⊕ ) orbiting in a circular trajectory around an M4.0 dwarf.Its orbital period is 8.5 d and its semi-major axis is 0.02279 ± 0.00061 au; these parameters are comparable to those of GJ 3988 b.Another example is GJ 3138 b (Astudillo-Defru et al. 2017), orbiting in a nearly-circular orbit (e = 0.11 +0.11  −0.07 ) around an M0 star.Its minimum mass (M b sin i = 4.18 +0.61  −0.59 M ⊕ ) is in the regime of super-Earths and mini-Neptunes.It orbits its host star with a short period (P = 5.974 ± 0.001 d, a = 0.057 ± 0.001 au), similarly to GJ 3988 b.Additionally, TOI-776 b (Luque et al. 2021), is a transiting planet orbiting a bright M1 V star.Its mass is M b = 4.0 ± 0.9 M ⊕ and its orbital period is 8.25 d.As the other examples, this planet orbits near to its host star (a = 0.0652 ± 0.0014 au).In Fig. 15 we display a (minimum) mass against orbital period plot of small (M < 15 M ⊕ ) and close-in (P < 100 d) confirmed exoplan-  ets, demonstrating that GJ 3988 b represents a common type of planet around M dwarfs.

Summary
We discovered GJ 724 b and GJ 3988 b using RV measurements from CARMENES, HARPS, and IRD.We measured the mean stellar rotation periods for both stars using photometric data, obtaining 116 ± 2 d for GJ 3988, and 57 ± 1 d for GJ 724, although we cannot completely exclude its second harmonic (∼ 28 d) as a possible rotation period.We also examined spectroscopic activity indicators and computed their GLS periodograms to look for stellar activity signals.Both planetary signals are dominant in the RV periodograms and thus easy to detect.We analysed their stability and concluded that they are not caused by stellar activity.We tested different configurations for both systems and found that the best models with the strongest evidences corresponded to single planets.The minimum masses of GJ 724 b and GJ 3988 b are 10.75 +0.96  −0.87 M ⊕ and 3.69 +0.42 −0.41 M ⊕ , while their orbital periods are 5.10 d and 6.94 d.The separations from their host stars are less than 0.05 au.GJ 724 b has a very eccentric orbit (e = 0.577 +0.055 −0.052 ), having the highest eccentricities of all single planets around M dwarfs.We discussed some possibilities for this configuration, such as another (undetectable) companion in the system, tidal in-spiral migration, perturbation or scattering event in the early stage of the system, or a Kozai effect.With the aim to constrain our current knowledge of planetary formation and architecture, we suggest to obtain more data to verify our measurement itself and to search for possible companions.We also encourage future dynamical interaction studies, which would also help to improve our understanding of the high eccentricity of GJ 724 b.On the other hand, GJ 3988 b orbits around its host star in a circular trajectory, as do a fair number of planets around M dwarfs, making it a common planet around these type of stars.
Acknowledgements.This publication was based on observations collected under the CARMENES Legacy+ project.CARMENES is an instrument at the Centro Astronómico Hispano en Andalucía (CAHA) at Calar Alto (Almería, Spain), operated jointly by the Junta de Andalucía and the Instituto de Astrofísica de Andalucía (CSIC).CARMENES was funded by the Max-Planck-Gesellschaft (MPG), the Consejo Superior de Investigaciones Científicas (CSIC), the Ministerio de Economía y Competitividad (MINECO) and the European Regional Development Fund (ERDF) through projects FICTS-2011-02, ICTS-2017-07-CAHA-4, and CAHA16-CE-3978, and the members of the CARMENES Consortium (Max-Planck-Institut für Astronomie, Instituto de Astrofísica de Andalucía, Landessternwarte Königstuhl, Institut de Ciències de l'Espai, Institut für Astrophysik Göttingen, Universidad Complutense de Madrid, Thüringer Landessternwarte Tautenburg, Instituto de Astrofísica de Canarias, Hamburger Sternwarte, Centro de Astrobiología and Centro Astronómico Hispano-Alemán), with additional contributions by the MINECO, the Deutsche Forschungsgemeinschaft (DFG) through the Major Research Instrumentation Programme and Research Unit FOR2544 "Blue Planets around Red Stars", the Klaus Tschira Stiftung, the states of Baden-Württemberg and Niedersachsen, and by the Junta de Andalucía.Data were partly collected with the 90 cm telescope at the Observatorio de Sierra Nevada operated by the Instituto de Astrofífica de Andalucía (CSIC).The authors wish to express their sincere thanks to all members of the Calar Alto and Sierra Nevada staff for their expert support of the instrument and telescope operation.We acknowledge financial support from the Agencia Estatal de Investigación (AEI/10.13039/501100011033) of the Ministerio de Ciencia e Innovación and the ERDF "A way of making Europe" through projects PID2021-125627OB-C31, PID2019-109522GB-C5[1:4], PID2019-107061GB-C64, PID2019-110689RB-100, CEX2021-001131-S, and the Centre of Excellence "Severo Ochoa" and "María de Maeztu" awards to the Instituto de Astrofísica de Canarias (CEX2019-000920-S), Instituto de Astrofísica de Andalucía (SEV-2017-0709) and Institut de Ciències de l'Espai (CEX2020-001058-M).This work was also funded by the DFG through grant DR 281/39-1, the Israel Science Foundation through grant No. 1404/22, the Generalitat de Catalunya/CERCA programme, the JSPS KAKENHI grant No.18H05442.We recognise the support from Silvia Sabotta for creating the detection map of GJ 724.We thank the anonymous referee for the comments that improved the quality of this manuscript.

Fig. 1 .
Fig.1.GLS periodograms of the (combined) photometric data and spectroscopic activity indicators from GJ 724 with peaks that reach at least a FAP level of 10% in the period range of 2 d to 200 d.The gray dashed horizontal lines correspond to the FAP levels of 0.1%, 1%, and 10% (from top to bottom, respectively).The mean rotation period of 56 d determined from the photometry is marked by the blue dashed line, while its second harmonic (P rot /2) of 28 d is depicted by the orange dashed line.The period of the 5.1-day planet is highlighted by the red solid line.

Fig. 2 .
Fig.2.Posterior distribution from dQ GP against P GP (top) and Q GP against P GP (bottom) from the GP fit to the photometric data of GJ 724 using a dSHO kernel.The colour bar displays the log-likelihood, while gray points correspond to the samples with ∆ ln L > 10.

Fig. 3 .
Fig.3.Same as Fig.1but for GJ 3988.The mean rotation period of 116 d determined from the photometry is marked by the blue dashed line.The period of the 6.4-day planet is highlighted by the red solid line.

Fig. 4 .
Fig.4.GLS periodograms of the window function (top panel), RVs (middle), and residual for the one planet model (bottom) for GJ 724 using CARMENES and HARPS data sets.The period of the planet, GJ 724 b, is highlighted by the solid red line, whereas its aliases at 1.24 d and 0.83 d are marked with dashed green and magenta lines, respectively.Additionally, we depict the mean rotation period of 57 d determined from the photometry and its harmonic P rot/2 by the dot-dashed blue and orange lines, respectively.

Fig. 5 .
Fig. 5. s-BGLS periodogram of the RV timeseries of GJ 724.The number of data points is plotted against the period.The top colour bar indicates the logarithm of the probability, whereas the side instrument bar shows from which instrument the data originate.CARMENES is shown in orangered and HARPS in bluish-green.The left panel is centred on the 5.1 d signal, middle panel around 56 d, and right panel around 28 d.
d-ecc, dSHO,28d) model.The RVs along this model are shown in Fig. B.1, while the RVs phased-folded to the planet period are shown in Fig. 6.The posteriors and derived parameters are listed in Table 5.A corner plot of the posteriors of the planet parameters is shown in Fig. B.2.

Fig. 6 .
Fig. 6.Phased RVs for GJ 724 b from the best fit model (1P (5 d-ecc) + dSHO-GP 28 d ).The red dots show the CARMENES data while the teal squares depicts the HARPS data.The black lines show the median of 10 000 samples from the posterior.The residuals after subtracting the median model are shown in the lower panel.

Fig. 7 .
Fig. 7. GLS periodograms and window function of the CARMENES and IRD RVs for GJ 3988.The period of the planet, GJ 3988 b, is highlighted by the solid line and the mean rotation period of 116 d determined from the photometry is marked by the blue dashed line.

Fig. 9 .
Fig. 9. Same as Fig. 6 but for the best-model for GJ 3988 b (1P (7 d-circ) + dSHO-GP 116 d ).The red dots show the CARMENES data while the teal squares depicts the IRD data.

Fig. 11 .
Fig. 11.Eccentricity plotted against semi-major axis of the confirmed exoplanets around single stars from the NASA Exoplanet Archive.The red colour represent planets around M-dwarfs, whereas the gray colour illustrates planets around any other stellar mass.Single planets (top panel) are represented by a star symbol, while planets in multi-planetary systems are represented by a square symbol (bottom).The blue star shows the position of GJ 724 b from the solution of the single planet in an eccentric orbit, which is also plotted faintly in the bottom panel (multi-planetary systems) just for visualisation purposes.

Fig. 12 .
Fig. 12.Detection limit map of the RVs of GJ 724.The colour map represents the detection probability of the period-mass combination.The Grid is 50 × 50 in mass and period, with 50 phase samples for each combination.The blue star indicates the planet GJ 724 b.

Fig. 13 .
Fig. 13.Synthetic planetary population as a function of M sin i and eccentricity.The colour shows the orbital period of the planets and the histogram depicts the fraction of theoretically expected observed planets, which is normalised by the number of planets in the same sample.The grey region corresponds to the eccentricity of GJ 724 b within 1σ range.

Fig. 14 .
Fig.14.Planetary mass against orbital period of theoretically computed planetary systems.The four systems are selected from the 1000 synthetic systems ofBurn et al. (2021) by filtering for observable planets with e > 0.4 and maximum displacement of 0.5 times the orbital period and mass of GJ 724 b.Planets with detection probabilities above 50% are plotted with bigger markers and thinner edges.Thin gray lines show the past history of each planet and triangles show planets accreted by the planet closest to GJ 724 b in mass-period space; for simulation 445, this is the second planet.Horizontal black error-bars indicate the distance from apastron to periastron of the planets measured in periods of circular orbits.Two pathways leading to different numbers of observable planets can be distinguished with only the pathway with recent accretion of a massive counterpart and a unobservable exterior chain of planets leading to a single detectable planet.

Fig. 15 .
Fig. 15.(Minimum) mass vs. orbital period plot of confirmed exoplanets from NASA Exoplanet Archive with orbital periods less than 100 d and minimum masses less than 15 M ⊕ .The red stars represents planets around M-dwarfs stars, whereas the gray colour depicts planets around other types of stars.The green star represents the position of GJ 3899 b in the diagram.For completeness, we also include GJ 724 b, represented as a blue star.

Table 2 .
Overview of the RV data.

Table 3 .
Overview of the photometric data.

Table 5 .
; Van Grootel 4 SHERLOCK (Searching for Hints of Exoplanets fRom Lightcurves Of spaCe-based seeKers) code is fully available on GitHub: https: //github.com/franpoz/SHERLOCKArticle number, page 10 of 28 P. Gorrini et al.: Planetary companions orbiting the M dwarfs GJ 724 and GJ 3988 Best-fit posterior parameters and determined planet parameters for the two planetary systems.