© The Authors 2025

1 Introduction

The number of confirmed exoplanets has steadily increased since the first discoveries were made nearly three decades ago. Although during the initial era of that research most exoplanets were found by the means of Doppler spectroscopy, carried out in long time-baseline programs (e.g., Vogt et al. 1994; Mayor et al. 2003), over recent years most of the discoveries have been made via space-based transit observations with missions such as CoRoT (Baglin et al. 2006), Kepler (Borucki et al. 2010), and TESS (Ricker et al. 2015).

The observed population of currently almost 6000 confirmed exoplanets and candidates serves as a probe for theoretical models on planet formation. Ideally, the synthetic planet populations from models should converge towards the observed one. Unfortunately, the task is not as trivial as it sounds. In order to reproduce the actual planet population by modeling, one needs to understand the physical process in its entirety. In particular, this also includes differentiation on dynamic time scales and stellar, or rather protoplanetary, disk masses. On the other hand, the outcomes of observational surveys heavily depend on the used methods, instruments’ sensitivity, target selection, as well as number and methodology of observations. Those choices of course lead to some selection biases, which need to be accounted for when determining the planet populations from them.

Naturally, the confidence in conclusions derived from observational population studies is highly affected by the number of considered objects. Therefore, the discovery and characterization of thousands of transiting planet candidates by Kepler was a crucial step forward. Since the mission focused on G-type stars, Kepler transiting planets around M dwarfs were highly underrepresented, and thus occurrence rates of planets orbiting M dwarfs relying on that survey come with high uncertainties (Hardegree-Ullman et al. 2019). Still, it was shown that the occurrence rate of small planets on short orbits up to 50 d around M dwarfs is higher than around solar-like stars (Howard et al. 2012).

Across the whole spectrum of stellar types, M dwarfs are particularly interesting to study. They are not only the most common type of stars (Reylé et al. 2021, and references therein), but also favorable, because due to their low mass and size, small and low-mass planets can be detected around them more easily using Doppler spectroscopy and the transit method, respectively. It was found very early on, and has been confirmed thereafter, that the planet population around M dwarfs differs from those around other types of host stars. As the occurrence rate of giant planets was observed to be correlated with stellar mass (Johnson et al. 2010), in particular hot Jupiters were believed to be rare around M dwarfs compared to their occurrence rate around hotter and more massive stars (Endl et al. 2006; Johnson et al. 2012; Hartman et al. 2015). While the actual percentage of cool stars hosting hot Jupiters is still under debate due to high uncertainty levels in occurrence rate analyses, as determined by, for example, Obermeier et al. (2016), it appears that the frequency of hot Jupiters peaks at G-type dwarfs and decreases in both directions, towards hotter and cooler stars (Gan et al. 2023).

Differential studies of planet occurrence across the full M-dwarf spectral type are particularly valuable because they span a wide range of masses, temperatures, and stellar environments, all of which can significantly influence planet formation and evolution. By examining the occurrence rates of planets around early-, mid-, and late-type M dwarfs, one can explore how factors such as stellar mass, disk properties, and stellar activity shape the frequency and characteristics of formed planets. These comparative studies help to identify important trends, such as whether low-mass planets are more prevalent around specific subtypes of M dwarfs and how stellar properties influence the potential habitability of orbiting planets. Such analyses do not only refine planet formation theories, but also improve our understanding of where potentially habitable planets are most likely to be found. Estimations of occurrence rates of planets around M dwarfs have only been reported for the earlier spectral types, of about M3.5–4.0 V, which corresponds to stellar masses down to around 0.33 M_⊙ (e.g., Bonfils et al. 2013a; Pinamonti et al. 2022). For M dwarfs of later spectral types, however, the picture becomes comparatively incomplete. Because of their faintness, it is more challenging to search for planets around late, low-mass M dwarfs.

Still, among M dwarfs, very low-mass stars (M ≲ 0.16 M_⊙) are of particular interest. Despite their smaller numbers in brightness-limited samples, their lower masses and cooler temperatures create unique conditions for planet formation, offering crucial insights into planetary systems that form in the least massive and faintest stellar environments. Understanding planet formation around these stars is essential for developing a complete picture of how planets form across the full range of stellar masses. Moreover, the low luminosity of very low-mass stars shifts the habitable zone to much closer orbits, facilitating not only the RV detection of Earth-like planets, but also demographic studies on climatic conditions within and outside of the habitable zone (e.g., Checlair et al. 2019; Turbet et al. 2019; Schlecker et al. 2024). These planets also orbit in environments where stellar activity and flaring may influence atmospheric retention and habitability (Tarter et al. 2007). Additionally, very low-mass stars provide advantageous targets for atmospheric characterization. The large planet-to-star size ratios and the proximity of habitable zones result in deeper transit signals and stronger atmospheric features in transmission spectroscopy, facilitating the study of the atmospheres of small planets, including those that could potentially harbor life (Trifonov et al. 2021; Kuzuhara et al. 2024).

The CARMENES¹ spectrograph, installed at the 3.5 m telescope of the Calar Alto Observatory in Almería, Spain, is ideally suited for investigating planetary systems around M dwarfs, including very low-mass stars. Specifically designed to conduct an RV survey of around 350 M dwarfs, the instrument has been operational since January 2016, covering both the visual (VIS) and near-infrared (NIR) wavelength ranges between 520 nm and 1710 nm with spectral resolutions of R = 94 600 in the VIS and R = 80 400 in the NIR (Quirrenbach et al. 2014). Soon after the first planet discovery findings, Sabotta et al. (2021) conducted a foundational study based on a sample of 71 M dwarfs observed with CARMENES, and reported an abundance of short-period planets, particularly around the latest-type M dwarfs. They highlighted the need for further investigation into how occurrence rates vary across the M dwarf spectral type. Based on these results, Ribas et al. (2023) expanded this analysis using a larger sample of 238 M dwarfs from the CARMENES guaranteed time observations (GTO). Their refined study confirmed the overabundance of short-period planets around late-type M dwarfs, further emphasizing the importance of differentiating planet occurrence rates by stellar mass. With this refinement, Ribas et al. (2023) provided an overall occurrence rate of 1.44 ± 0.20 planets per star, illustrating that nearly every M dwarf hosts at least one planet, while also revealing significant trends related to stellar mass.

Our present study serves as a specific follow-up to the studies by both Sabotta et al. (2021) and Ribas et al. (2023), focusing exclusively on the least massive stars (M ≲ 0.16 M_⊙, spectral type about M5.5 V and later) from the CARMENES survey. With the focus on Earth-like planets that are detectable by current RV surveys, we studied only companions with orbital periods of up to 100 d and planetary masses below 10 M_⊕. By refining the occurrence rates of Earth-like planets in this low- mass regime, we took an important first step toward completing the full picture of planetary distribution around M dwarfs. Although the current stellar sample is small and subject to statistical uncertainties, this work provides a primary basis for a more comprehensive analysis when the full CARMENES survey is complete. Our ultimate goal is to provide refined occurrence rates of Earth-like planets specifically for very low-mass stars, contributing to a broader understanding of planetary system formation across the entire M-dwarf spectrum.

The selection criteria for the targets of this study, as well as the final stellar sample itself, are summarized in Sect. 2. Four newly discovered exoplanets around three stars of the sample are presented in Sect. 3. In Sect. 4 the methods and steps for estimating the occurrence rates are documented and explained. Thereafter, in Sect. 5, the results are discussed with respect to their uncertainties and robustness of the methods, and are put into context by comparing them to previous analyses on the subject. In addition to that, the distribution of detected planets around the targets in our stellar sample are compared to predictions from planet formation theory. Finally, we summarize and conclude our work in Sect. 6.

2 Stellar sample

For this study, our aim was to intensively observe stars for which we could detect Earth-mass planets. For this reason, we set several constraints on our sample. As a baseline catalog, we took the CARMENES input catalog, namely Carmencita (version 106; Caballero et al. 2016a). Carmencita contains several dozens parameters for about 2200 nearby, bright M dwarfs, from which the GTO targets were selected. An updated summary of Carmencita was provided by Cortés-Contreras et al. (2024). From there, we selected all targets with masses M ≤ 0.1617 M_⊙. We chose this exact value because an Earth-mass planet in an orbit of 10 d around a star with this stellar mass would induce an RV amplitude of 1 m s⁻¹ , which, based on our experience, is the minimum amplitude that can be detected by CARMENES with a reasonable number of observations (Zechmeister et al. 2019; Luque et al. 2022; Kemmer et al. 2022; Kossakowski et al. 2023; Suárez Mascareño et al. 2023). Stellar masses in Carmencita are computed following the methodology of Schweitzer et al. (2019), which is determining bolometric luminosities from the integration of the spectral energy distribution and precise Gaia parallaxes (Cifuentes et al. 2020), effective temperatures from spectral synthesis (e.g., Passegger et al. 2022, and references therein), stellar radii from the two parameters above and Stefan-Boltzmann law, and stellar masses from a mass-radius relationship calibrated with double-lined, detached, eclipsing binaries.

Additionally, we applied a brightness threshold of J ≤ 10 mag, because for fainter stars a precision of 1 m s⁻¹ typically cannot be achieved unless an extraordinarily large number of RVs is collected. We also included only inactive (or, rather, very weakly active) targets, as activity can reduce the detection efficiency, and it is moreover difficult to characterize proper detection sensitivities for those targets. We set three constraints on the activity and kept stars with pseudo-equivalent of the Hα line pEW(Hα) > −1.5 Å, which removed around half of the targets within the mass limit, P_rot > 10 d, and υ sin i < 2 km s⁻¹ . These thresholds made sure that we were not including young active stars with high chromospheric emission, which are expected to also be fast rotators. Less conservative limits on the three parameters were applied to studies of M dwarfs in the CARMENES sample by Schöfer et al. (2019), Cortés-Contreras et al. (2024), and Kemmer et al. (2025).

The 15 remaining targets defined our sample and are listed in Table 1, together with their fundamental parameters and number of valid RVs collected by CARMENES. The stellar rotation periods were collected from the literature, as referenced in the Carmencita catalog (P_rot,ref column in Table 1). They were used in the following to identify activity-induced RV signals and to distinguish them from planetary companions in our analysis. Only for one of our targets, G 109−35 (Karmn J06594+193), the rotation period was not known. Using available photometric data and spectroscopic activity indicators, we determined the missing rotation periods and arrived at $P_{\text{rot }}={110}_{-13}^{+16} \text{d}.$ This analysis is presented in Appendix A. One more star with spectral type M5.0 V, namely GJ 3250 (Karmn J03473+086), also met our defined requirements but was left out of the subsequent analysis, as it had been observed infrequently with CARMENES and only 12 RV measurements were available at the time of the analysis.

Fig. 1 Window function (upper panel) and GLS periodograms of the CARMENES RVs for G 268−110 before (‘0P’, middle) and after subtracting the one-planet model (‘1P’, bottom). The two panels on the left and right represent the same GLS periodograms but plot different regions to better represent the occurring signals. The period of the 1.43-day planet is highlighted by a red solid line. Its first-order aliases at 3.28 d and 0.59 d are marked by red dashed lines, and the second-order aliases at 0.76 d and 0.43 d by red dotted lines. The rotation period of 143 d determined by Newton et al. (2018) is indicated by a solid green line.

3 Planet discoveries

3.1 RV data

The RV measurements used for this work were collected from the M-dwarf survey carried out with the CARMENES spectrograph. It was specifically designed to deliver highly accurate RVs with a precision of the order of 1 m s⁻¹ to search for temperate rocky planets around nearby cool stars (Ribas et al. 2023). Although the instrument provides RVs from two separate spectral channels, we utilized only data from the VIS channel, covering the spectral region up to 960 nm. All acquired spectra were reduced by the extraction pipeline caracal (CARMENES Reduction And CALibration; Caballero et al. 2016b), and the RVs were determined by means of template matching with serval (SpEctrum Radial Velocity AnaLyser; Zechmeister et al. 2018). In order to correct for uncalibrated systematic effects, shared by RVs from the same night, nightly zero-point offsets (NZPs), as described by, for example, Trifonov et al. (2018) and Tal-Or et al. (2019), were computed and applied. All observations of our stellar targets that were used for this publication were obtained between 12 January 2016 and 6 February 2024.

3.2 Discovery of G 268−110 b

From six years of CARMENES observations, we report the discovery of the Earth-mass planet G 268−110 b (GJ 1028, J01048−181). We first derived a tentative orbital period of P = 1.432 630(76) d and a minimum planetary mass of M_pl sin i = 1.52(25) M_⊕ . However, the RV data suffer from strong aliasing, as is evident in the Generalized Lomb-Scargle (GLS) periodogram (Zechmeister et al. 2009) depicted in Fig. 1. It shows a set of very significant signals with low false alarm probabilities (FAPs) at periods of 1.43 d (FAP < 6.1 × 10⁻⁵), 0.59 d (FAP < 7.8 × 10⁻⁵), and 3.28 d (FAP < 1.6× 10⁻⁴), which are related by aliasing due to a sampling frequency of ∼1 d⁻¹ that is dominant in our RV data. This aliasing is further evident through the also significant second-order aliases (with respect to the 1.43- day signal) at periods of 0.76 d and 0.43 d. After subtracting the 1.43-day signal, which has the lowest FAP of all, there are no other significant signals present.

Unfortunately, the signals are all of comparable low FAP in the GLS periodogram, which makes the determination of the true orbital period very difficult. We applied the AliasFinder code, which allows comparing synthetically generated periodograms for different alias periods with the observed one (Dawson & Fabrycky 2010; Stock & Kemmer 2020; Stock et al. 2020). The best match can be an indication of the likely underlying period of the signal. However, as can be seen in Fig. B.1, the periodograms resulting from the different first-order alias periods remain indistinguishable in our case. Therefore, we performed a model comparison for the three different periods, where we also allowed for non-zero eccentricity in the planetary orbit. We tested all of them against the base model (= “0P”), which only includes a jitter and an offset of the CARMENES data.

We used juliet (Espinoza et al. 2019) to perform the fits to the RV data. Thus, Keplerian orbits were parameterized by their period P, RV semi-amplitude K, time of inferior conjunction t₀ , and the ${\mathcal{S}}_1=\sqrt{e}\text{\hspace{0.17em}sin\hspace{0.17em}}\omega$ and ${\mathcal{S}}_2=\sqrt{e}\text{\hspace{0.17em}cos\hspace{0.17em}}\omega$ parameters, which depend on the eccentricity e and argument of periastron ω. For each of the three tested periods we tried both circular and eccentric orbits and constrained the priors for each signal to a narrow range around the highest local peak in order to also avoid issues from yearly aliasing, which is visible in Fig. B.1. An overview of the priors used is given in Table C.1. The results of this test are tabulated in Table 2.

We found that irrespective of the chosen period, all 1P- models are favored over the base model (∆ In 𝒵 > 5; Trotta 2008), while eccentric models are indistinguishable (∆ In 𝒵 < 3) from circular models. We therefore assumed a circular orbit because e is not well determined by our data and thus could lead to an overestimation of its value (Hara et al. 2019, and references therein).

Regarding the different possible periods, ln 𝒵 is only a valid metric for models with the same period priors. Therefore, for the comparison between the different periods we used the maximum log-likelihood (ln ℒ) for the circular or eccentric models individually. However, it turned out that there are no significant differences between the individual maximum likelihoods. In combination with the alias test, this means that it is not possible to identify the true period of the signal from the current data set. A concentrated campaign with high-cadence observations, as would be necessary for the resolution of the 0.5-day alias, is unfortunately difficult to implement with CARMENES due to the low declination of G 268−110.

Still, since the 1.43-day signal shows the lowest FAP in the GLS and slightly outperforms the models to the other periods in terms of ln ℒ, we considered that period as the most probable one. The corresponding phase-folded RVs are plotted in Fig. 2, the complete RV time series can be found in Fig. B.2, and the model parameters, as well as the derived planetary parameters, are listed in Table 3. However, for the sake of completeness, the parameters determined from the sampling of the alternative periods are listed in Table B.1.

For the rotation period, we relied on an earlier measurement of 143 d by Newton et al. (2018) from photometric data from the MEarth project (Irwin et al. 2009). Although this is far from the signals detected in the RVs, we nevertheless explored the CARMENES activity indicators to exclude other forms of activity that are not related to the stellar rotation as the origin of the RV signals. The list of activity indicators studied included all of those regularly computed in the CARMENES data flow, such as the chromatic index (CRX), differential linewidth (dLW), CCF bisector (BIS), and contrast (CON), as well as indicators related to the pEWs of specific lines such as Hα, TiO, and more (see Zechmeister et al. 2018; Schöfer et al. 2019; Lafarga et al. 2020, 2021, for the full list). In order to identify common periodicities appearing in the set of indicators, we first scanned the GLS periodograms of the activity indicators for common periods applying the DBSCAN clustering algorithm (Density-Based Spatial Clustering of Applications with Noise; Ester et al. 1996) implemented in scikit-learn (Pedregosa et al. 2011), as described by Kemmer et al. (2025). For this, we calculated the GLS periodogram for each activity indicator and identified the ten highest peaks. The periodograms start at a period of 2 d to avoid aliases that occur due to the typical sampling of daily observations. The DBSCAN algorithm was then run on the combined sample of all peaks with FAPs below 80 % and periods shorter than the baseline of the observations. In the process, we considered a minimum number of three samples in a neighborhood to be a cluster if the maximum frequency distance between two detected peaks is less than half of the width of the peaks in the GLS (i.e., δf = (t_max − t_min)⁻¹). This clustering analysis is illustrated in Fig. B.3. We did not detect any significant clusters of periods that would hint at a strong influence from stellar activity onto our data. All found clusters are related to the harmonics of one year or the Moon cycle, and none overlap with the planetary signal or its aliases in the RVs. This is also reflected in the GLS periodograms of the activity indicators with signals less than 10% FAP (see Fig. B.4). We therefore reasonably assumed that stellar activity does not need to be considered in our modeling.

Fig. 2 Phased RV plot for G 268–110 b based on the best fit model (1P(1.43 d-circ)). The black line depicts the model based on the parameters listed in Table 3. The shaded area illustrates the 1σ confidence interval.

Fig. 3 Same as Fig. 2 but for G 261−6 b based on the best fit model (1P(5.45 d-circ)).

3.3 Discovery of G 261−6 b

Our intensive RV monitoring campaign of G 261−6 (GJ 1238, J19242+755) over a time span of two years revealed a small temperate planet, which orbits its host star with a period of P = 5.4536(32) d and a minimum mass of M_p1 sin i = 1.37(23) M_⊕. The complete set of its model parameters is given in Table 3, the phase-folded plot of the RVs to the single Keplerian model is illustrated in Fig. 3, and the RVs over time are plotted in Fig. B.5. As depicted in Fig. 4, its signal is highly significant in the GLS periodogram of the system, with FAP < 2.8 × 10⁻⁶. Also prominent are further signals that we attributed to its daily aliases at periods of 1.22 d and 0.84 d with FAP < 0.17% and FAP < 0.28%, respectively. These peaks are no longer notable in the residuals of the model after subtracting the planet signal. However, the residuals suggest three additional signals at periods of 0.95 d, 22.49 d, and 1.83 d. The most prominent of them, at 0.95 d, has a FAP below 1.4% and is related to the 22.49 d signal by daily sampling, while the 1.83 d signal with the lowest significance of FAP < 23% appears to be independent of the others.

We further investigated the consistency of these additional additional signals using stacked-Bayesian GLS (s-BGLS) periodograms (Mortier & Collier Cameron 2017). They are depicted in Fig. B.6 and all of them show significant variability over time. Due to this incoherence and not being sufficiently significant in the GLS, we did not consider them in our model, but since they could turn out to be promising candidates, the system should be further monitored. More data will be helpful to reveal their true nature.

For G 261−6, we also explored the possibility that the additional signals could be spurious and induced by activity. As the analysis of possible period clustering (see Fig. B.7) and the GLSs of CARMENES activity indicators with signals of less than 10 % FAP (see Fig. B.8) show, apparently none of the discussed periodic signals are related to the 114 d rotation period of G 261–6 (Irwin et al. 2011), nor other forms of activity that would imprint periodic signals onto our RVs.

Fig. 4 Same as Fig. 1 but for G 261−6. We indicate periods of the planet (5.45 d), its aliases (1.22 d and 0.84 d), and stellar rotation (114 d) determined by Irwin et al. (2011). Besides, the additional signals at periods of 10.95 d, 1.83 d, and 22.49 d are marked by blue dashed lines.

3.4 Discovery of G 192–15 b and G 192–15 c

The CARMENES RV measurements of G 192−15 (GJ 3380, J06024+498) collected over more than seven years reveal two substellar companions in the system. The GLS periodogram (see Fig. 5) shows a highly significant long periodic peak with an FAP < 9.0 × 10⁻¹⁰ at around 1289 d together with its daily alias at 0.997 d and FAP < 7.2 × 10⁻⁶. After modeling this signal with an eccentric orbit, a second significant signal in the residuals becomes apparent at a period of 2.27 d with an FAP < 2.2 × 10⁻⁴. The one-day alias of this secondary peak at 1.78 d is also evident, but of less significance. We tested different models with varying number of planets, as well as allowing for non-zero eccentricities, and concluded that there are two planets orbiting the stellar host G 192−15. The evidences for the various models are tabulated in Table 4.

The long period planet exhibits a high minimum mass of M_p1 sin i = 14.3 M_⊕ and revolves around its host star on an eccentric orbit with e ≈ 0.68. For the short-period planet we derived a minimum mass of M_pl sin i = 1.03(18) M_⊕ and our model comparison favors a circular Keplerian over an eccentric one. The complete set of model parameters are listed in Table 3, and the phase-folded RVs together with the entire RV time series can be found in Fig. 6.

After subtracting our two-planet solution from the RVs, we found another signal in the residuals at a period of 69.94 d. While it does not reach a high level of significance (FAP around 15 %), a three-planet solution appears at least fairly reasonable with Δ In 𝒵 = 1.6 against our preferred model. For further investigation, we examined all the signals’ stability with increasing number of observation. As depicted in the s-BGLSs in Fig. B.9, the two planetary signals appear coherent and increase in significance over time, whereas the third candidate signal shows some variability. In addition to that, we used Exostriker (Trifonov 2019) to run an orbital evolution for such a three-planet configuration and found that, while the system is stable for the two planets alone, they would undergo significant changes in their orbital eccentricities in the presence of the third companion. For those reasons, we did not consider the third signal any further, but again additional monitoring of the system would be helpful for its understanding.

We also analyzed the activity indices with the most significant signals, as we did for the previous systems. The corresponding GLS periodograms (Fig. B.10), however, do not reveal any significant power at the periods of the discovered companions, nor at the rotation period of the star at 105 d (Díez Alonso et al. 2019).

Fig. 5 Same as Fig. 1 but for G 192–15. Extra panels at the bottom are for the GLS after subtracting the two-planet model. We indicate the periods of the two planets at 1219 d and 2.27 d (red solid lines), their aliases (red dashed lines), and stellar rotation (105 d) determined by Díez Alonso et al. (2019) (green solid line). Besides, the additional signal at a period of 69.9 d is marked by the blue dashed line.

4 Planet occurrences

Planet occurrence rates are estimated from the number of detected planets within particular mass or period bins around a sample of investigated stars. The known planets around the stars from the investigated stellar sample, together with the ones presented in this work, are listed in Table 5. Wherever applicable, the listed planetary masses are rescaled by stellar masses updated since the original publications, as provided in Table 1.

Since the planet candidates proposed by Suárez Mascareño et al. (2023) around G 158–27 (GJ 1002, J00067–075) can not be identified from the CARMENES RVs alone, they were not taken into account in the further statistical analysis regarding planet occurrences. Beside that, as in this work we studied only companions with orbital periods of up to 100 d and planetary masses below 10 M_⊕ , any known planets above these limits were also omitted, which left us with 11 planets in seven systems (i.e. all in Table 1 except for J00067–075 b and c, J06024+498 c, and J08413+594 b and c).

The numbers of planet discoveries are generally affected by the capabilities of the instruments used for the observations of the targets, as well as the methods applied for data reduction and eventually signal and planet detection. Naturally, unknown planets can be missed in RV surveys, which if uncorrected, can lead to a bias and an underestimation of the final occurrence rates. To overcome this issue, we followed the procedures described by Sabotta et al. (2021) and estimated the planet detection completeness within our stellar sample by means of an injection- and-retrieval analysis similar to those previously used (Cumming et al. 1999; Zechmeister et al. 2009; Meunier et al. 2012; Bonfils et al. 2013b; Wolthoff et al. 2022).

Fig. 6 Upper panel: Same as Fig. 2 but for the simultaneous best fit model (= 2P(1218 d-ecc, 2.27 d-circ)) of the two planetary signals around G 192–15. Lower panel: RVs over time for the same model. The black lines show the model based on the parameters listed in Table 3, and the gray shaded areas denote the 68%, 95%, and 99% confidence intervals, respectively. The instrument offset of CARMENES was subtracted from the measurements and the model.

4.1 Preprocessing RV data

As it is preferable to perform the injection-and-retrieval analysis for the estimation of the detection completeness (see Sect. 4.2) on clean RV data, we applied a prewhitening procedure on the time series contained in the analyzed sample. To do so, as in Sabotta et al. (2021), we first performed 3σ outlier clipping and computed GLS periodograms for all targets. We then repeatedly fitted Keplerian orbits to signals with a FAP of 1% or lower, until no such signals were left. In this manner we retrieved the signals of the known planets, but also of others that we attributed either to activity, stellar rotation, or telluric contamination. Activity signals were identified by using the time series of the activity indicators evaluated by serval. On the other hand, if the signals’ periods coincide with known stellar rotation periods or their first high harmonics, we attributed them to stellar rotation. Finally, systematics due to telluric contamination were identified by comparing the RV data with those determined from the same spectra after the correction for telluric absorption following Nagel et al. (2023). However, for some of the found signals, their origin cannot be settled with the current data available to us. All signals found and identified following the described approach are listed in Table 6.

4.2 Detection completeness

In order to estimate the detection sensitivity within our data sets, we injected single artificial planets on circular orbits into our RV data and tested whether we were able to retrieve them. In this step, we neglected eccentricity, since it introduces an additional complexity to the problem, but, as implied by Cumming (2004) and verified by Sabotta et al. (2021), does not significantly change the outcome. We counted individual tests as successful recoveries if the signal in question appears as the highest peak in the resulting GLS periodogram with a FAP below 1 %. We did this 200 times for all combinations of orbital periods P and planet minimum masses M_pl sin i spanned by a log-uniform grid of 30 points in the period range 1–100 d and 30 grid points in the planet minimum mass range 0.5–10 M_⊕, i.e., 180 000 different simulations in total. For all these planetary signals, the phase ϕ was chosen at random, and the semi-amplitude was determined by the common approximation $$K=28.435 \text{m} {\text{s}}^{-1}{\left(\frac{P}{1\text{\hspace{0.17em}yr}}\right)}^{-1/3}\left(\frac{M_{\text{pl}}\sin i}{M_{\text{Jup}}}\right){\left(\frac{M_{\star }}{M_{\odot }}\right)}^{-2/3}.$$(1)

Consequently, the RV signals of the artificial planets $$RV(t)=K\sin \left(\frac{2\pi t}{P}+\varphi \right)$$(2)

were evaluated at the given time stamps of the observations and added up onto the RV time series after they had been prewhitened, as described in Sect. 4.1. The individual detection maps on the chosen grid of parameters for each tested star were next determined by the fraction of successfully recovered planetary signals. After averaging those detection maps over all of the 15 stellar targets, we obtained the final detection sensitivity map of the survey’s subsample, which is plotted in Fig. 7.

Fig. 7 CARMENES detection sensitivity map, averaged over the individual maps of the 15 stellar targets of this study. The light blue markers indicate the 11 planets included for this occurrence rate analysis, and the color map illustrates the detection probabilities of the respective periodmass grid points. The solid lines represent the masses associated to the RV semi-amplitude values equal to the RVs’ root mean square averaged over the time series (blue), the mean RV uncertainties (magenta), and the median of the RV uncertainties.

4.3 Planet occurrence rates

We determined the final occurrence rates within the period range 1–100 d and minimum mass range 0.5–10 M_⊕ by running a Monte Carlo simulation on a grid of planet frequencies ${\overline{n}}_{\text{pl}}$, in terms of number of planets per star, to test how consistent it is with the actual number of detected planets N_pl,det. The simulation was run 1000 times iteratively for each of the test frequencies and consisted of the following steps:

• (a)

  Given the test frequency ${\overline{n}}_{\text{pl}}$, a test planet sample is created, for which the number of test planets N_pl,in is determined by the Poisson distribution $\lambda ={\overline{n}}_{\text{pl}}{N}_{\star }$, where N_✶ is the number of stars within our underlying stellar sample.

• (b)

  A period and minimum mass from the mass-period grid of the detection sensitivity map are assigned to each of the test planets. We applied a log-uniform distribution in the period, whereas for the mass we used a power-law distribution of the form N_pl = a(M_pl sin i)^α, with a = 319.25 and α = −1.06 (Ribas et al. 2023). To test the robustness of the occurrence rate estimations to the assumption of the underlying planet mass distribution, we alternatively applied a log-uniform distribution in mass as well.

• (c)

  The number of test planet retrievals N_pl,out is determined by the count of detections within the sample of the test planets. Here, each of the test planets was accepted as a successful detection based on a random draw with a binary chance of success according to its corresponding detection probability.

For each of the test frequencies ${\overline{n}}_{\text{pl}}$ we counted how often, out of the 1000 simulation runs, the number of retrieved test planets equals the actual number of known detected planets N_pl,out = N_pl,det. Based on these counts, the resulting probability density over the given grid of tested frequencies was finally normalized, and from the locations of its median, as well as its 16th and 84th percentiles, we derived the most probable frequency (occurrence rate) and its corresponding uncertainties, as illustrated by Fig. 8.

Fig. 8 Occurrence rate distributions for different combinations of period and planetary mass ranges. The medians as well as the 16th and 84th percentiles are indicated by vertical solid and dashed lines, respectively. The distributions in blue are based on an underlying power-law distribution for the planet masses, while the orange ones are estimated using a log-uniform distribution.

5 Results and discussion

5.1 Late M dwarfs host Earth-like planets and very few super-Earths

We determined the occurrence rates for four distinct, not overlapping bins in the mass-period plane with orbital periods P between 1 d and 100 d and minimum planetary masses M_pl sin i between 0.5 M_⊕ and 10 M_⊕. They are listed in Table 7, and the underlying distributions are plotted in Fig. 8. Our results indicate a significant dependency on planetary masses. Whereas the occurrence rates for low planetary masses below 3 M_⊕ , namely ${\overline{n}}_{\text{pl}}={0.88}_{-0.28}^{+0.36}$ for short orbital periods below P = 10 d and ${\overline{n}}_{\text{pl}}={0.92}_{-0.39}^{+0.56}$ for periods above P = 10 d, agree well with one planet per star, they decrease substantially for planetary masses between 3 M_⊕ and 10 M_⊕ to ${\overline{n}}_{\text{pl}}={0.11}_{-0.06}^{+0.11}\text{ and }{\overline{n}}_{\text{pl}}={0.06}_{-0.04}^{+0.09}$, respectively. This tendency is illustrated in Fig. 9. For planets below 3 M_⊕ the occurrence rates that we obtained appear to be constant with respect to orbital periods. Therefore, we cannot confirm the implications from Sabotta et al. (2021), who found increased occurrence rates for orbital periods below 10 d for stellar hosts of masses below 0.34 M_⊙, which is twice as high as the mass threshold for the stellar sample that we investigated. While our occurrence rates for the planetary mass range between 3 M_⊕ and 10 M_⊕ may indicate such a trend, with only one known planet within this mass range, the statistics in this regime are not particularly reliable. Consequently, the corresponding uncertainties are high, and the values obtained are mutually consistent within their errors.

Although the occurrence rates are consistently lower when using a log-uniform distribution for the planetary masses, they agree within their uncertainties to those estimated from the power-law distribution, and the relations with periods and masses prevail. However, the rates for planets with masses below 3 M_⊕ using the power-law distribution appear to be significantly increased with respect to those based on the log-uniform distribution when wider orbits are considered. This is due to an increased number of low-mass planets on wide orbits in the simulated planet samples. As those low-mass planets have small RV semi-amplitudes, they naturally fall into regions of low detection probabilities. This in turn leads to higher, possibly overestimated occurrence rates. A more realistic and applicable distribution of planetary masses with a dependence on the orbital period would be helpful to compensate for that possible bias.

All in all, our determined occurrence rates are increased for smaller planet masses for any of the investigated planetary orbit regimes. In fact, the most massive planet in our sample, with M_pl sin i = 3.95 ± 0.43 M_⊕, revolves around the most massive stellar host within our sample, namely CD Cet (J03133+047), with a mass of 0.161 M_⊙. In addition, the comparison to the occurrence rates from Ribas et al. (2023), which were derived from the complete CARMENES stellar sample at the time, indicates a significant dependence of planetary rates on spectral type, even within the M-dwarf regime alone. For planetary masses between 1 M_⊕ and 10 M_⊕ , the authors reported significantly lower occurrence rates of ${\overline{n}}_{\text{pl}}={0.39}_{-0.07}^{+0.10}$ for short orbits below 10 d than for wider orbits up to 100 d, with ${\overline{n}}_{\text{pl}}={0.67}_{-0.15}^{+0.18}$. In contrast to that, the results from the present work, for which only the least massive stars were considered, indicate a significantly increased number of planets on short orbits up to 10 d and for planetary masses from 0.5 M_⊕ to 10 M_⊕. For that same shortperiod regime we consequently arrive at an occurrence rate of ${\overline{n}}_{\text{pl}}$ = 0.99 planets per star, which is at least twice as high as that reported by Ribas et al. (2023). This number can be derived from the combination of the rates for low-mass planets below 3 M_⊕ ${\overline{n}}_{\text{pl}}=0.88$, and ${\overline{n}}_{\text{pl}}=0.11$ for masses between 3 M_⊕ and 10M_⊕, as they represent independent bins of planet masses. This dependence on stellar masses is also evident when our rates are compared with those in other studies. For low-mass planets up to 10 M_⊕ and on orbits between P = 1 d and P = 100 d, Bonfils et al. (2013a) found a rate of ${\overline{n}}_{\text{pl}}={0.36}_{-0.10}^{+0.24}$ for their HARPS M- dwarf sample with a median stellar mass of 0.33 M_⊙, which is in line with Ribas et al. (2023). An even lower rate of ${\overline{n}}_{\text{pl}}={0.10}_{-0.03}^{+0.10}$ was derived by Pinamonti et al. (2022) for an even earlier stellar sample of M dwarfs with types between M0 and M3 and planets on orbits shorter than P = 10 d. Those trends again suggest that the number of small planets on short orbits is increased for late spectral types.

To illustrate the difference of the analyzed samples, in Fig. 10 the distribution of stellar masses of our sample is plotted and put into context to the subsample from Ribas et al. (2023) and to the entire CARMENES sample. While the subsample from Ribas et al. (2023) represents the entire GTO sample quite accurately, the stars used for the current work are significantly less massive. We presume that with a sufficiently large pool of targets and observations, it will be possible to further resolve trends in planet demographics and their dependence on host star mass, clarifying the role of host star properties on the planet formation process.

Still, one should take note of the high uncertainties of the estimated rates, which are due to the small number statistics. The analyzed sample contains only 15 stars and 11 already detected planets that lie within our range of period of interest and that could be identified with CARMENES RVs alone. Also, although we tried to compensate by our methods for the incompleteness of the sample, it cannot be ruled out that the number of low-mass planets on wide orbits is still underestimated, as the detection sensitivity for Earth-like planets is still fairly low, even for orbits at around P = 20 d (see the sensitivity map in Fig. 7).

For the sake of completeness, we note that there may be a small number of false-negative non-detections in RV surveys for planets on orbits with periods that are close to the stellar rotation periods of their hosts. Signals of such candidates tend to be falsely discarded, as the data often show concurrent stellar activity. From transit surveys the percentage of such cases can be estimated to be around 5% (Lechuga et al. in prep.) and therefore we conclude that the planet occurrence rates will be underestimated only insignificantly.

Fig. 9 Occurrence rates based on the power-law distribution in planetary masses with respect to orbital periods. The different bins in planetary masses are color and symbol coded.

Fig. 10 Distributions of stellar masses for our stellar sample (green), the subsample used by Ribas et al. (2023) for their occurrence rate analysis (orange), and the whole CARMENES GTO sample (blue). The median stellar masses, indicated by the vertical lines, are 0.137 M⊙, 0.372 M⊙, and 0.348 M⊙ , respectively.

5.2 Detection probabilities from likelihood maps

As the accurate estimation of the occurrence rates is heavily dependent on a proper characterization of the detection limits, we tested another method to derive the sensitivity map for the sample. For this we followed a variation of the approach introduced by Tuomi et al. (2014) and later applied also by Pinamonti et al. (2022). Whereas they used posterior sampling for the estimation of the detection probabilities, we relied on log-likelihood maps in order to identify the area on the plane spanned by planetary masses and orbital periods, where a planetary signal with the corresponding semi-amplitude would likely be detected. To achieve that, for each of our cleaned and prewhitened RV time series (see Sect. 4.1), we estimated the likelihoods of the best fits of models with an additional postulated circular Keplerian signal. This was done on a two-dimensional grid for different orbital periods and planetary masses (RV semi-amplitudes), while the phase was treated as a free parameter. In this way, we arrived at a likelihood map for each target, where for each period the model likelihood decreases with increasing planetary mass. We evaluated the differences in ln ℒ across the whole grids using the models with semi-amplitudes compatible with zero as null hypotheses and, consequently, the average of their corresponding likelihoods as the baseline. From these differences, we derived a binary detection probability map using a threshold of Δ In ℒ = 5 with respect to the null hypothesis. Given the implication that in regions with differences in log-likelihood below that threshold an additional signal cannot be ruled out and, therefore, could be hidden in the data, the corresponding planets could not be detected. Consequently, analogously to Tuomi et al. (2014), we assumed that signals from planets within the complementary area of the mass-orbit plane, with differences in log-likelihood above the threshold, would likely be detected if they were apparent in the given time series. Therefore, we applied a detection probability of 100% for Δ In ℒ > 5 and 0% for Δ In ℒ < 5. This threshold corresponds to a probability ratio between two distinct models of about 150, and has been used in RV data analyses in the past (e.g., Feroz et al. 2011; Gregory 2011; Tuomi 2012). The individual binary maps were then averaged across the entire grid and the stellar sample in order to arrive at the final sensitivity map, which is illustrated in Fig. 11. The alternative occurrence rates based on these detection probabilities were calculated thereafter following the same procedure as outlined in Sect. 4.3, and the results are listed also in Table 7.

Although the absolute numbers differ slightly from those derived using the sensitivity map based on injection-and- retrieval, they agree within their uncertainties and, moreover, they also show the same dependencies with orbital periods and planetary masses. All of this is plausible, as the trends over different period and mass intervals are mostly governed by our actual sample of detected planets, while the magnitude of the absolute numbers depends highly on the sensitivity map used to account for the detection bias. The differences in sensitivity between the two approaches are depicted as differences of detection probabilities P_inj.-retr. – 𝒫_lik. in Fig. 11.

From this comparison, a general trend is apparent for orbital periods above around 5 d. For a given period, the detection probabilities based on the likelihood maps tend to be a bit lower for smaller planetary masses but exceed those based on injection- and-retrieval for higher masses. For shorter orbital periods (P_pl < 5 d) the alternatively derived probabilities are typically higher. The mean of the differences over the entire grid is −2.9%. The comparably lower detection probabilities from the likelihood maps at low planetary masses but wider orbits lead to increased occurrence rates in this region of the grid, whereas the rates are slightly lower than those from the injection-and-retrieval maps anywhere else. This trend is also confirmed by comparison of the probabilities averaged along the paths of the error limits, as can be seen in Figs. 7 and 11.

We stress the fact that, although the final numbers differ slightly, the results are still compatible, as both methods allow for some margin. Therefore, it is difficult to determine which method is superior. Both methods involve choosing a somewhat arbitrary threshold as detection criteria. While for the injection- and-retrieval method one needs to choose a threshold for the FAP of signals being counted as retrieved, one is left with an arbitrary choice of a difference in the Δ ln ℒ distribution at which one assumes a model not compatible with the null hypothesis any more. The limit of Δ ln ℒ = 5 that we chose translates into a probability ratio of 1/150% = 0.67% (see above), which is comparable to the 1% that we used in the injection-and- retrieval approach. Given these considerations, we believe that the compatibility of the resulting occurrence rates from the two methods underlines their robustness, in general. Still, one should keep in mind the margins given by the choice of the different applied methods when comparing the occurrence rates of different surveys and studies.

Fig. 11 Left: same as Fig. 7 but based on log-likelihood estimation and averaged over the individual binary maps of the 15 targets. Right: delta 𝒫inj.-retr. − 𝒫lik. between the CARMENES detection sensitivity maps based on injection-and-retrieval vs. the one based on log-likelihood estimation.

5.3 Comparison to planet formation theory

A comprehensive synthesis of planet formation around M dwarfs in the standard core accretion scenario was presented by Burn et al. (2021). They used the model of Emsenhuber et al. (2021), which integrates the growth by solid (in the form of planetesimals) and gas accretion, orbital migration, as well as N-body interactions of 50 concurrently growing protoplanets. At the start of the simulations, small planetesimals and the largest seed protoplanets are assumed to have formed. The initial disk properties are informed by Class I disk measurements and scaled to lower stellar masses. Based on this model, Burn et al. (2021) found occurrence rates of Earth-like planets of order unity in agreement with our results presented here. Giant planets, however, as reported by Morales et al. (2019), challenge planet formation models, as such ones cannot form around late M dwarfs under the standard assumptions (Schlecker et al. 2022): the required planetesimal densities are not reached at orbital distances with sufficiently short growth timescales (Schlecker et al. 2021). This conundrum might be solved by introducing disk structures acting as traps for migration (e.g., Hasegawa & Pudritz 2011).

Under the standard assumption that accretion of planetesimals and N-body interactions between low-mass planets leading to giant impacts dominate in the inner disk region, we compared the model predictions with our determined occurrence rates. As in Schlecker et al. (2022), we randomly drew inclinations for the systems synthesized by Burn et al. (2021) and created 1000 mock observations of 11 synthetic planets per draw. The chance of drawing each planet was weighted by the M_pl sin i and orbital period-dependent detection sensitivities mentioned in Sect. 4.2 and shown in Fig. 7.

The resulting distribution of synthetic planetary M_pl sin i and orbital periods was then compared to the known small planets detected using CARMENES data only, which are listed in Table 5. The results shown in Fig. 12 indicate some significant differences: The hypothesis that P_pl and M_pl sin i are individually drawn from the same distribution is excluded using a 1D Kolmogorov-Smirnov test at moderate p values of 0.068 and 0.008, respectively. The observed planets orbit their host stars at orbital periods shorter than predicted, and their masses are smaller. In fact, due to the dependency of the detection sensitivity on orbital period and mass, both distributions can no longer be statistically distinguished if synthetic planetary masses are reduced by a factor of two (p values of 0.5 and 0.3 for P_pl and M_pl sin i, respectively).

As planet mass is the most fundamental and constraining property for their formation, a disagreement about the population of small planets demands an explanation. While a thorough model iteration with different parameters is out of the scope of this work, we used analytical mass scales derived by Emsenhuber et al. (2023) to assess what changes in initial conditions would lead to consistent planetary masses. Assuming that the planets did not migrate over large separations, the relevant mass scale is the Goldreich mass $$M_{\text{Gold}}=16a^3{\text{Σ}}_{\text{P}}^{3/2}{\left(\frac{2{\pi }^7{a}^3{\rho }_{\text{p}}}{3M_{\star }^3}\right)}^{1/4},$$(3)

where a is the semi-major axis, ∑_P is the surface density of planetesimals, and ρ_P is the bulk density of the rocky planets and planetary embryos. A reduction of the amount of planetesimals by 37% is sufficient to decrease the expected mass of the planets by 50% and, therefore, to the range of the observed planet mass distribution. For solar-mass stars, a comparison between the synthetic data from the same planet formation model and the HARPS survey results (Mayor et al. 2011) shows no evidence that such a reduction is needed (Emsenhuber et al. 2024). Moreover, the planetary mass range can be probed to lower values than what HARPS can achieve around solar-type stars. Combined, those results hint at a steeper-than-linear dependency of the available mass in planet-building material with stellar mass.

Future studies should further explore whether a possible lower planetesimal surface density for the lowest-mass stars can be explained by planetesimal formation models. The results from Lenz et al. (2019) indicate even higher surface densities than those assumed by Burn et al. (2021) due to the predicted steeper slope of the radial surface density of about −2, in contrast to the adopted −1.5. However, Burn et al. (2021) assumed globally complete conversion of the dust mass to planetesimals, which is perhaps overly optimistic. Moreover, the used dust masses were observed by Tychoniec et al. (2018) in the Perseus association, which seems to be one of the star-forming regions with more massive disks as opposed to, for example, Ophiuchus (Williams et al. 2019, see also the discussion by Tobin et al. 2020). Furthermore, if preferential locations exist for planetesimal formation, such as the water ice line (Drążkowska & Alibert 2017), a different planet formation channel is required for close-in planets.

An alternative scenario to explain the differences between the observed and predicted planetary mass and orbit distributions is migration from the exterior of the snowline, where different mass scales can apply, namely the typical mass at which type-I migration dominates over the relevant solid accretion mechanism (M_pl ≈ 5 M_⊕.). If this is the major origin of the planets, they should contain large amounts of H₂O and other volatile species, which are expected to be retained during photoevaporation (Lopez 2017; Burn et al. 2024). So far, transiting planets in the mass range of the low-mass planets in Table 3 are consistent with being rocky (Luque & Pallé 2022).

Another possibility is that the growth of planets is dominated by pebble accretion. In this case, the pebble isolation mass, which depends on the uncertain temperature and viscosity of the disk, can leave a distinct imprint on the mass distribution (Brügger et al. 2020) similar to what is observed in Fig. 12. However, the work of Liu et al. (2019) focusing on late M dwarfs found lower planetary masses than found here (while the recent work of Venturini et al. 2024 found rocky planets with masses up to 3 M_⊕). Therefore, subsequent giant impacts, excluded by Liu et al. (2019), would be required between the planets at the pebble isolation mass, which would re-introduce scatter in the planetary masses.

Finally, an additional effect omitted in the simulations is the loss of material to debris during giant impacts. Although the simulations by Burn et al. (2021) assumed perfect merging, about a factor of two in mass is lost due to fragmentation in realistic giant impacts (Emsenhuber et al. 2020). When this mass loss is taken into account, the outcome agrees remarkably well with our results and provides an additional avenue to explain the population of small planets around very late M dwarfs, in addition to reduced planetesimal surface densities and dry pebble accretion with subsequent migration.

Fig. 12 Planet detection statistics of observed planets compared with mock detections based on the synthetic planet population by Burn et al. (2021). Included are all planets with orbital periods of 1–100 d and with minimum masses of 0.5–10 M⊕. The left panel shows the Mpl sin i of the observed (orange squares) and synthetic (circles) planets against orbital periods. The transparency and color of the synthetic data are scaled by detection probability. The cumulative distribution of the observed and synthetic Ppl (middle) and Mpl sin i (right) is shown. For the synthetic data, we show with transparent lines 1000 random draws of 11 planets to visualize the expected spread.

6 Summary

We studied planet occurrence rates in a sample of 15 late-type M dwarfs observed with the CARMENES spectrograph. All the stars are low-mass (M ≲ 0.16 M_⊙), relatively bright (J < 10 mag), slowly-rotating (υ sin i < 2 km s⁻¹), and weakly-active (pEW(Hα) > −1.5 Å). We used available photometric data combined with time series of different spectral activity indicators from CARMENES to determine the rotation period for one of the targets without a previous measurement.

In order to minimize possible detection biases in our final occurrence rates, we performed an injection-and-retrieval analysis on the RV data of our targets to determine the detection sensitivities with respect to all data sets at hand. In the process of reevaluating the RV time series and identifying all significant periodic signals, we confirmed ten known planets around five M dwarfs and discovered four new planets around three stars: G 268-110b, with a minimum mass of M = 1.52(25) M_⊕ on a short orbit of P ≈ 1.43 d around its host, G 261–6 b, a 1.37(23) M_⊕ planet with a period of P ≈ 5.45 d, and two companions around G 192–15. One of them is an Earth-mass planet with 1.03 ± 0.18 M⊕ at P ≈ 2.27 d and the other is a 14.3 ± 1.6 M_⊕ planet on a wide and eccentric orbit with P ≈ 1220 d and e = 0.68 ± 0.07.

The final planet occurrences were determined by Monte Carlo simulations within predefined orbital period and minimum mass ranges. We took into account the overall detection probabilities averaged over the individual targets’ detection sensitivity maps in order to account for false statistics from missing detections. We found an apparent trend in planet occurrence as a function of minimum planet masses. The rates decrease substantially from ${\overline{n}}_{\text{pl}}={0.88}_{-0.28}^{+0.36}$ for planetary masses between 0.5 M_⊕ and 3 M_⊕ to ${\overline{n}}_{\text{pl}}={0.11}_{-0.06}^{+0.11}$ for planetary masses between 3 M_⊕ and 10 M_⊕ in the orbital period regime from P = 1 d to P = 10 d, and from ${\overline{n}}_{\text{pl}}={0.92}_{-0.39}^{+0.56}$ to ${\overline{n}}_{\text{pl}}={0.06}_{-0.04}^{+0.09}$ respectively for orbital periods between P = 10 d and P = 100 d. As our results do not show any significant dependency on orbital periods, we could not confirm the general trend of increasing occurrence rates when moving from wider to shorter orbits, as indicated by Sabotta et al. (2021) in a larger mass regime of M dwarfs (M = 0.095–0.34 M_⊙). However, in comparison to previous studies (e.g., Bonfils et al. 2013a; Hardegree-Ullman et al. 2019; Pinamonti et al. 2022), we showed that stars of later spectral types tend to host a larger number of small planets on shorter orbits. This trend is apparent among M dwarfs, as for our sample of targets with the lowest masses we found at least twice as many planets on shorter orbits than indicated by the occurrence rates from Ribas et al. (2023), who analyzed the entire CARMENES sample, back then covering the mass range M = 0.095–0.677 M_⊙.

Since occurrence rates are highly affected by possible detection biases, a proper estimation of the apparent detection limits is crucial. Although for our analysis we relied on an injection- and-retrieval analysis, other methods, such as posterior sampling have been applied for this purpose in the past (Tuomi et al. 2014; Pinamonti et al. 2022). Therefore, we tested another method, based purely on likelihood maps. While the resulting sensitivity maps differ from our initial approach, the final occurrence rates are robust against the choice of method and agree within their uncertainties. The rather high uncertainties in our rate estimations presumably stem mainly from low-number statistics, as our final stellar sample consists of only 15 targets with 11 known planets that satisfy our selection criteria. All except one of these planets have minimum masses around or below 1.5 M_⊕, and in particular below the minimum mass of water worlds, so that their composition is probably rocky (Luque & Pallé 2022). However, with only one planet in the highest mass bin, the rates for the mass bin above 3 M_⊕ should be handled with caution.

To put our results in context, we compared the currently known planets around the stars within our stellar sample with predictions from state-of-the-art planet formation models. We simulated mock observations of planets drawn from a synthetic planet population around M dwarfs based on the standard core accretion scenario (Burn et al. 2021), and weighted them by our estimated detection probabilities. The observed orbital periods are shorter and the minimum masses are smaller than predicted. In the framework of standard core accretion, such a deviation could in principle be explained by a reduced planetesimal surface density in models. While such a correction was not needed based on previous occurrence rate studies, we probed here the very low-mass star regime. Therefore, the dependency of available mass in planet building blocks on stellar mass needs to be further investigated. Apart from that, the observed lower masses could also be explained by alternative formation scenarios such as accretion of dry pebbles within the snowline or mass loss during giant impacts.

Altogether, our results and the discussion show the importance of differentiation in stellar masses when discussing planet occurrences and, ultimately, planet formation processes. Once the CARMENES survey of M dwarfs is completed and at least 50 RV epochs have been obtained for all targets, we will conclude our analysis of occurrence rates on the entire sample.

Acknowledgements

We thank the anonymous referee for a very quick and constructive report. 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 Min- isterio 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, Insti- tut 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 Forschungsge- meinschaft (DFG) through the Major Research Instrumentation Programme and Research Unit FOR2544 “Blue Planets around Red Stars” (RE 2694/8-1), the Klaus Tschira Stiftung, the states of Baden-Württemberg and Niedersachsen, and by the Junta de Andalucía. We used data from the CARMENES data archive at CAB (CSIC-INTA). We acknowledge financial support from the Agencia Estatal de Investigación (AEI/10.13039/501100011033) of the Ministe- rio de Ciencia e Innovación and the ERDF “A way of making Europe” through projects PID2022-137241NB-C4[1:4], PID2021-125627OB-C31, 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 (CEX2021-001131-S) and Institut de Ciències de l’Espai (CEX2020- 001058-M). This work was also funded by the Generalitat de Catalunya/CERCA programme, the DFG under Germany’s Excellence Strategy EXC 2181/1390900948 Exploratory project EP 8.4 (Heidelberg STRUCTURES Excellence Cluster), and the Bulgarian National Science Fund (FNI) program “VIHREN- 2021” project No. KP-06-DV/5. The results reported herein benefitted from collaborations and/or information exchange within NASA’s Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA’s Science Mission Directorate under Agreement No. 80NSSC21K0593 for the program “Alien Earths”.

Appendix A Rotational period of G 109−35

For the determination of the stellar rotational period of G 109−35 (Karmn J06594+193) we used all available photometric data and spectroscopic activity indicators. The activity indicators are measured from CARMENES spectra and are available as a time series of 30 data points corresponding to our RV measurements collected over four years. As two observations did not provide information on drift corrections during the night, the corresponding precise RVs were discarded, while the determination of the activity indicators is still valid.

To search for common periods in the GLS periodograms of the indicators, we applied the DBSCAN clustering algorithm as described in Sect. 3.2. The result of this analysis is illustrated in Fig. A.1. It revealed two clusters of signals at periods of around 95 d and 125 d, but with only three activity signals being significant at FAPs of 4.0 × 10⁻⁷−2.0 × 10⁻⁴ and 8.8 × 10⁻⁷, respectively. We used the positions of those signals as prior knowledge in our further analysis.

We then analyzed photometric data from the All-Sky Automated Survey (ASAS²; Pojmanski 1997) and the Northern Sky Variability Survey (NSVS³; Woźniak et al. 2004) as compiled by Díez Alonso et al. (2019). The ASAS data available to us cover around seven years between December 2002 and November 2009, and were taken by the V-band wide-field camera at the survey's station located in Las Campanas Observatory, Chile (ASAS-3), while the NSVS photometric measurements that we used were collected in Los Álamos, USA, between April 1999 and April 2000. An inspection of the corresponding GLS periodograms did not reveal any significant signals at the periods of interest, as already reported by Díez Alonso et al. (2019). We therefore continued with Gaussian process (GP) regression fits to the combined photometric data using a double harmonic oscillator (dSHO) kernel. A description of the kernel and its application to photometric data using Juliet was provided by Kossakowski et al. (2021). The default priors for the parameters of the kernel used in this modeling are listed in Table C.2. During a first step, we used a wide unconstrained prior for the period to look for promising signals, of U(10,200) d. The posterior distribution shows evidence for periods between 90 d and 130 d, as well as some at a period of around half a year. Next, using our prior findings from the analysis of the activity indicators, we constrained the range of possible rotation periods to further sample the regions around the two clusters in Fig. A.1. The shape of the posterior of this model is bimodal with peaks at periods of 101 d and 119 d. Although the shorter period appears to be slightly more favored in this distribution, we still continued to test both periods with models, for which the priors in periods were further narrowed, in order to sample those peaks individually. The results of these models yielded rotational periods of P_rot,ı = 101.9 ± 5.0 d and P_rot,2 = 119.0 ± 7.2 d, respectively. Unfortunately, the difference between the likelihoods of the two models is marginal and, therefore, does not indicate which of them is favored by the given data. For this reason, we declared both periods as genuine candidates for the true rotation period but noted that activity appears to be slightly more significant at the shorter one, as well as that it is slightly more favored during sampling with priors allowing both of them. Eventually, we report a single value $P_{\text{rot}}={110}_{-13}^{+16}$ d in Table 1, which is determined by averaging the two candidate periods. Their original errors were used to derive the final uncertainties.

Fig. A.1 Diagram of peaks appearing in the GLS periodograms of the activity indicators accessible with CARMENES for G 109−35 (J06594+193). Big open circles at 94.7 d (in blue) and 124.8 d (brown) mark the center periods of the two main DBSCAN clusters. Vertical dashed lines indicate one year and its higher harmonics, while the horizontal dotted line indicates FAP = 0.1 %.

Appendix B Analyses on planet discoveries

B.1 Discovery of G 268−110 b

Fig. B.1 Alias test for G 268−110. Each row represents the result of simulating a different underlying period: first row: 1.43 d, second row: 3.28 d, third row: 0.59 d (each marked by the blue dashed line, respectively). The solid black line shows the median periodogram from 1000 simulations, with the interquartile ranges and the ranges of 90 % and 99 % denoted by the different gray-shaded areas. The observed periodogram is depicted by the red solid line. In addition to the periodograms, the resulting peak phases can also be compared in the circles above (the same colors as for the periodograms, but with the gray shade showing the standard deviation of the sampled peak phases).

Fig. B.2 RVs over time for the best model (= 1P(ı.4 d-circ)) fitted to the CARMENES RVs of G 268−110. The black lines show the model based on the parameters listed in Table 3, and the gray shaded areas denote the 68 %, 95 %, and 99 % confidence intervals, respectively. The instrumental offset of CARMENES was subtracted from the measurements and the model.

Fig. B.3 Diagram of peaks appearing in the GLS periodograms of the activity indicators accessible with CARMENES for G 268−110. For each activity indicator, the GLS periodogram was created and the ten highest-occurring peaks determined. If the FAP of the GLS peak is below 80 %, it is written to a table. This list of peaks is used to run a DBSCAN clustering algorithm. A cluster is a group of peaks with at least 3 members, where the distance to the nearest neighbor is less than the resolution of the GLS.

Fig. B.4 GLS periodograms of the activity indicators with signals of less then 10 % FAP for G 268−110. The period of the 1.4-day planet and the 3.3-day alias are highlighted by the red solid and dashed lines, respectively. The rotation period of 143 d determined by Newton et al. (2018) is marked by the green solid line.

B.2 Discovery of G 261-6 b

Fig. B.5 Same as Fig. B.2 but for the best model (= 1P(5d-circ)) fitted to the CARMENES RVs of G 261−6.

Fig. B.6 Stacked Bayesian GLS periodograms of the three noteworthy signals with periods around 0.95 d, 1.83 d and 22.49 d, detected in the residuals of the RVs after the planetary signal is removed.

Fig. B.7 Diagram of peaks appearing in the GLS periodograms of the activity indicators accessible with CARMENES for G 261−6. For each activity indicator, the GLS periodogram was created and the ten highest-occurring peaks determined. If the FAP of the GLS peak is below 80 %, it is written to a table. This list of peaks is used to run a DBSCAN clustering algorithm. A cluster is a group of peaks with at least 3 members, where the distance to the nearest neighbor is less than the resolution of the GLS.

Fig. B.8 Same as Fig. B.4 but for the activity indicators with significant signals for G 261−6. The period of the 5.5-day planet is highlighted by the red solid line, while the locations of the possible candidate signals in the residuals at 0.95 d, 1.83 d and 22.49 d are indicated by the blue dashed lines, respectively. The rotation period of 114 d determined by Irwin et al. (2011) is marked by the green solid line.

B.3 Discovery of G 192−15 b and G 192−15 c

Fig. B.9 Stacked Bayesian GLS periodograms of the signals detected in the RVs. Left panel: Zoom in on the long-period planetary signal at 1218.5 d. Middle panel: Zoom in on the 2.27 d period signal. Right panel: Zoom in on the 69.94 d variable signal, which appears in the RV residuals after the 2-planet model is subtracted.

Fig. B.10 Same as Fig. B.4 but for the activity indicators with significant signals for G 192−15. The periods 2.27 d and 1218.5 d of the planet signals are highlighted by the red solid lines, and the signal found in the residuals at 69.94 d is indicated by the blue dashed line, respectively. The rotation period of 105 d determined by Díez Alonso et al. (2019) is marked by the green solid line.

Appendix C Model priors

References

All Tables

All Figures

Fig. 1 Window function (upper panel) and GLS periodograms of the CARMENES RVs for G 268−110 before (‘0P’, middle) and after subtracting the one-planet model (‘1P’, bottom). The two panels on the left and right represent the same GLS periodograms but plot different regions to better represent the occurring signals. The period of the 1.43-day planet is highlighted by a red solid line. Its first-order aliases at 3.28 d and 0.59 d are marked by red dashed lines, and the second-order aliases at 0.76 d and 0.43 d by red dotted lines. The rotation period of 143 d determined by Newton et al. (2018) is indicated by a solid green line.

In the text

	Fig. 2 Phased RV plot for G 268–110 b based on the best fit model (1P(1.43 d-circ)). The black line depicts the model based on the parameters listed in Table 3. The shaded area illustrates the 1σ confidence interval.
In the text

	Fig. 3 Same as Fig. 2 but for G 261−6 b based on the best fit model (1P(5.45 d-circ)).
In the text

	Fig. 4 Same as Fig. 1 but for G 261−6. We indicate periods of the planet (5.45 d), its aliases (1.22 d and 0.84 d), and stellar rotation (114 d) determined by Irwin et al. (2011). Besides, the additional signals at periods of 10.95 d, 1.83 d, and 22.49 d are marked by blue dashed lines.
In the text

Fig. 5 Same as Fig. 1 but for G 192–15. Extra panels at the bottom are for the GLS after subtracting the two-planet model. We indicate the periods of the two planets at 1219 d and 2.27 d (red solid lines), their aliases (red dashed lines), and stellar rotation (105 d) determined by Díez Alonso et al. (2019) (green solid line). Besides, the additional signal at a period of 69.9 d is marked by the blue dashed line.

In the text

Fig. 6 Upper panel: Same as Fig. 2 but for the simultaneous best fit model (= 2P(1218 d-ecc, 2.27 d-circ)) of the two planetary signals around G 192–15. Lower panel: RVs over time for the same model. The black lines show the model based on the parameters listed in Table 3, and the gray shaded areas denote the 68%, 95%, and 99% confidence intervals, respectively. The instrument offset of CARMENES was subtracted from the measurements and the model.

In the text

Fig. 7 CARMENES detection sensitivity map, averaged over the individual maps of the 15 stellar targets of this study. The light blue markers indicate the 11 planets included for this occurrence rate analysis, and the color map illustrates the detection probabilities of the respective periodmass grid points. The solid lines represent the masses associated to the RV semi-amplitude values equal to the RVs’ root mean square averaged over the time series (blue), the mean RV uncertainties (magenta), and the median of the RV uncertainties.

In the text

	Fig. 8 Occurrence rate distributions for different combinations of period and planetary mass ranges. The medians as well as the 16th and 84th percentiles are indicated by vertical solid and dashed lines, respectively. The distributions in blue are based on an underlying power-law distribution for the planet masses, while the orange ones are estimated using a log-uniform distribution.
In the text

	Fig. 9 Occurrence rates based on the power-law distribution in planetary masses with respect to orbital periods. The different bins in planetary masses are color and symbol coded.
In the text

	Fig. 10 Distributions of stellar masses for our stellar sample (green), the subsample used by Ribas et al. (2023) for their occurrence rate analysis (orange), and the whole CARMENES GTO sample (blue). The median stellar masses, indicated by the vertical lines, are 0.137 M⊙, 0.372 M⊙, and 0.348 M⊙ , respectively.
In the text

	Fig. 11 Left: same as Fig. 7 but based on log-likelihood estimation and averaged over the individual binary maps of the 15 targets. Right: delta 𝒫inj.-retr. − 𝒫lik. between the CARMENES detection sensitivity maps based on injection-and-retrieval vs. the one based on log-likelihood estimation.
In the text

Fig. 12 Planet detection statistics of observed planets compared with mock detections based on the synthetic planet population by Burn et al. (2021). Included are all planets with orbital periods of 1–100 d and with minimum masses of 0.5–10 M⊕. The left panel shows the Mpl sin i of the observed (orange squares) and synthetic (circles) planets against orbital periods. The transparency and color of the synthetic data are scaled by detection probability. The cumulative distribution of the observed and synthetic Ppl (middle) and Mpl sin i (right) is shown. For the synthetic data, we show with transparent lines 1000 random draws of 11 planets to visualize the expected spread.

In the text

	Fig. A.1 Diagram of peaks appearing in the GLS periodograms of the activity indicators accessible with CARMENES for G 109−35 (J06594+193). Big open circles at 94.7 d (in blue) and 124.8 d (brown) mark the center periods of the two main DBSCAN clusters. Vertical dashed lines indicate one year and its higher harmonics, while the horizontal dotted line indicates FAP = 0.1 %.
In the text

Fig. B.1 Alias test for G 268−110. Each row represents the result of simulating a different underlying period: first row: 1.43 d, second row: 3.28 d, third row: 0.59 d (each marked by the blue dashed line, respectively). The solid black line shows the median periodogram from 1000 simulations, with the interquartile ranges and the ranges of 90 % and 99 % denoted by the different gray-shaded areas. The observed periodogram is depicted by the red solid line. In addition to the periodograms, the resulting peak phases can also be compared in the circles above (the same colors as for the periodograms, but with the gray shade showing the standard deviation of the sampled peak phases).

In the text

	Fig. B.2 RVs over time for the best model (= 1P(ı.4 d-circ)) fitted to the CARMENES RVs of G 268−110. The black lines show the model based on the parameters listed in Table 3, and the gray shaded areas denote the 68 %, 95 %, and 99 % confidence intervals, respectively. The instrumental offset of CARMENES was subtracted from the measurements and the model.
In the text

Fig. B.3 Diagram of peaks appearing in the GLS periodograms of the activity indicators accessible with CARMENES for G 268−110. For each activity indicator, the GLS periodogram was created and the ten highest-occurring peaks determined. If the FAP of the GLS peak is below 80 %, it is written to a table. This list of peaks is used to run a DBSCAN clustering algorithm. A cluster is a group of peaks with at least 3 members, where the distance to the nearest neighbor is less than the resolution of the GLS.

In the text

	Fig. B.4 GLS periodograms of the activity indicators with signals of less then 10 % FAP for G 268−110. The period of the 1.4-day planet and the 3.3-day alias are highlighted by the red solid and dashed lines, respectively. The rotation period of 143 d determined by Newton et al. (2018) is marked by the green solid line.
In the text

	Fig. B.5 Same as Fig. B.2 but for the best model (= 1P(5d-circ)) fitted to the CARMENES RVs of G 261−6.
In the text

	Fig. B.6 Stacked Bayesian GLS periodograms of the three noteworthy signals with periods around 0.95 d, 1.83 d and 22.49 d, detected in the residuals of the RVs after the planetary signal is removed.
In the text

Fig. B.7 Diagram of peaks appearing in the GLS periodograms of the activity indicators accessible with CARMENES for G 261−6. For each activity indicator, the GLS periodogram was created and the ten highest-occurring peaks determined. If the FAP of the GLS peak is below 80 %, it is written to a table. This list of peaks is used to run a DBSCAN clustering algorithm. A cluster is a group of peaks with at least 3 members, where the distance to the nearest neighbor is less than the resolution of the GLS.

In the text

Fig. B.8 Same as Fig. B.4 but for the activity indicators with significant signals for G 261−6. The period of the 5.5-day planet is highlighted by the red solid line, while the locations of the possible candidate signals in the residuals at 0.95 d, 1.83 d and 22.49 d are indicated by the blue dashed lines, respectively. The rotation period of 114 d determined by Irwin et al. (2011) is marked by the green solid line.

In the text

	Fig. B.9 Stacked Bayesian GLS periodograms of the signals detected in the RVs. Left panel: Zoom in on the long-period planetary signal at 1218.5 d. Middle panel: Zoom in on the 2.27 d period signal. Right panel: Zoom in on the 69.94 d variable signal, which appears in the RV residuals after the 2-planet model is subtracted.
In the text

	Fig. B.10 Same as Fig. B.4 but for the activity indicators with significant signals for G 192−15. The periods 2.27 d and 1218.5 d of the planet signals are highlighted by the red solid lines, and the signal found in the residuals at 69.94 d is indicated by the blue dashed line, respectively. The rotation period of 105 d determined by Díez Alonso et al. (2019) is marked by the green solid line.
In the text