© The Authors 2025

1 Introduction

M dwarfs represent 80% of the solar neighbourhood (Henry et al. 1994; Chabrier & Baraffe 2000), with more than 1500 stars within 26 pc, making them ideal targets for statistical analysis. As the radial velocity (hereafter RV) signal of a planet strongly depends on the stellar-to-planet mass ratio, large surveys of M dwarfs have been monitored with an RV planet searcher (e.g. HARPS, Mayor et al. 2003). Since the discovery of the first exoplanet orbiting an M dwarf (Delfosse et al. 1998; Marcy et al. 1998), which was a massive planet, a total of 135 planets have been discovered around M dwarfs using the RV technique¹. They span a wide range of masses, from 0.75 M_⊕ (GJ 54.1 b Astudillo-Defru et al. 2017a) to 1570 M_⊕ (GJ 676A b, 4.9 M_J Forveille et al. 2011; Anglada-Escudé & Tuomi 2012). The study of all these detections as a statistical set provides properties that individual detections and characterisations of planetary systems cannot offer, such as planetary and system occurrence rates. These properties are derived from large samples, which can only be obtained with two detection methods today for short and intermediate periods: photometric transits and RVs. The two Kepler transit missions yielded robust and detailed occurrence statistics of planets around FGK and M dwarfs based on orbital period and radii, especially for periods of under 100 days (e.g. Dressing & Charbonneau 2013; Mulders et al. 2015; Dressing & Charbonneau 2015; Hsu et al. 2020). They predict an increase in occurrence frequency with decreasing star radius. Concerning less massive planets (below 1 R_⊕), the occurrence rates published vary between 50% and 100% for periods of under 10 days and even more so between 10 and 100 days. Morton & Swift (2014) and Dressing & Charbonneau (2013) computed the number of planets per star of radii smaller than 4 R_⊕ orbiting massive M dwarfs and found, respectively, 2 ± 0.45 telluric planets and 0.9${}_{-0.03}^{+0.04}$ below, respectively, periods of 150 and 50 days.

More recently, Ment & Charbonneau (2023) established new occurrence rates of planets around very low-mass M dwarfs using the first 42 TESS sectors. In their study, they compared results with Dressing & Charbonneau (2015) and obtained a higher ratio of rocky-to-sub-Neptune planets for low-mass M dwarfs than for the massive ones.

Regarding the RV method, the extensive dataset of M dwarfs observed with HARPS has previously been the subject of a first occurrence analysis statistics established by Bonfils et al. (2013a), following the method previously proposed by Cumming et al. (2008) and Mayor et al. (2011) fora sample of FGK dwarfs. Based at the time on a sample of 102 M dwarfs observed over more than six years with HARPS, they already announced the preponderance of short-period telluric planets around M dwarfs. Sabotta et al. (2021) established a new planetary occurrence statistics using a sample of 71 M dwarfs observed over five years with the CARMENES visible and infrared spectrometers. Despite the smaller star sample, this new analysis reaffirmed the dominance of short-period terrestrial planets reported. More recently, in their study of 56 massive M dwarfs (M0-M3) observed over six years with HARPS-N, Pinamonti et al. (2022) reported a significantly lower occurrence rate compared to all previously cited studies. They obtained only 0.10 planets per star for periods of between 1 and 10 days.

Theoretical models of planet formation do not yet offer a comprehensive interpretation of the increasing planet occurrence rate around low-mass stars. A lower occurrence of more massive planets around these objects can only partially explain the reported higher abundance of low-mass planets around low-mass stars (Burn et al. 2021). However, detailed comparisons of the same theoretical model to the CARMENES survey have not revealed a statistically significant under-prediction of low-mass planets around M stars (Schlecker et al. 2022).

The statistics that we present in this paper is built on a sample of 197 M dwarfs observed with HARPS over the last two decades. Consequently, we delve into the long-period range that had remained unexplored by previous RV studies due to their limited temporal coverage. This research is a continuation of the analysis of these 197 RV time series (Mignon et al. (2024), hereafter paper I), which allows us to propose a substantial advancement in the methodology employed for constructing the occurrence statistics. This study also aims to investigate the relationship between planetary population and the mass of the central star by comparing the populations obtained around the most and least massive M dwarfs in our sample.

The outline of this paper is the following. In Sect. 2, we present the sample used to build the planetary occurrence statistics. In Sect. 3, we describe the method used to compute the conservative and statistical detection limits. In Sect. 4, we present the implementation of the planetary and system occurrence statistics, and we discuss and compare the results in Sect. 5. We compute the occurrence in the temperate zone and present the results in Sect. 6. Finally, we discuss the limits of this study in Sect.7 and conclude in Sect. 8.

2 Sample

This study focuses on 197² M dwarfs selected from paper I, which were observed for at least ten nights using the HARPS instrument between 2003 and 2021. We discuss the impact of this threshold in Sect. 5.4. This sample was meticulously constructed by cross-referencing several catalogues of cool dwarfs³ with the HARPS public archive of ESO⁴ and only considering targets within a distance of less than 26 pc and below a declination of +20°δ.

2.2 Magnitude and mass distributions

In paper I, we described how we selected VGJHK magnitudes from six different catalogues. The absence of a consistent V-band magnitude catalogue for all M dwarfs in our sample necessitates employing only three out of the five empirical massluminosity relations from Delfosse et al. (2000) (relations in the H, K, and J bands only) for stellar mass calculations. As presented in paper I, the mass distribution in our sample ranges from 0.75 to 0.1 M_⊙ with an average (resp. median) value of 0.44 M_⊙ (resp. 0.46 M_⊙).

2.2 Data

All the programmes whose measurements were used in this study are listed in Table A.1. RV time series are extracted using a specific template-matching code designed for M dwarfs that were previously established and used in multiple exoplanet detection studies (e.g. Astudillo-Defru et al. 2015, 2017b, 2020). Following a sequence of successive selection steps, our sample is made up of 197 series comprising 10 513 RVs that were binned by night and corrected for secular acceleration. For this second stage of analysis, we exclusively used the residual RV time series established in paper I to compute the detection limits, and we did not correct them any further.

2.3 Residual time series of RV

Regarding the long-term variability, 16 RV time series were corrected for a dominant linear model and 29 for a quadratic model. Once the long-term model is identified and corrected for, 36 RV time series are then corrected for one or multiple Keplerian signals attributed to either previously published planetary systems (listed in Table A.2) or arising as candidates (7 cases). Finally, we corrected 23 RV time series for the identified signature of the rotation period, which included not only those at the rotational period, but also harmonics. For example, in the case of GJ 581, we detected signals around 127 days in three chromospheric indices (Calcium, Sodium, Hα), but we detected and subtracted a RV signal found at 67 days; we attributed this to the harmonic of the rotation period. We employed a similar approach for GJ 447 (61.9 days and 123 days) and GJ 514 (15.15 days and 31.4 days).

2.4 Final sample

Ultimately, since the targets GJ 1132 and GJ 1214 were included in the blind search following the discovery of a planetary candidate by Mearth, we decided to exclude them from our sample to avoid biasing the statistics. Additionally, in identifying GJ 4254 as a spectroscopic binary - which prevents us from accurately isolating the signal from the primary companion - we excluded it.

3 Detection limits

We determined detection limits using the methodology outlined in Bonfils et al. (2013a) based on the generalised Lomb-Scargle (hereafter GLS) method proposed by Zechmeister & Kürster (2009). Given that all detectable signals within the RV time series are removed, we consider our residual RV time series a realisation of white noise.

3.1 Computation method

Following the approach presented in Bonfils et al. (2013a), we first shuffle the RV values of a given star while keeping the original temporal sampling and then compute the GLS periodogram. This is performed 1000 times and allows us to determine the distribution of power for each period that is compatible with no planet (i.e. due to noise). We then considered each period P in the periodogram. We inject a planetary signal with a period of P and no eccentricity that is of a given semi-amplitude K in the RV time series: 12 realisations of equi-spaced phases are performed. We start with K corresponding to the fitted sinusoid at Pin the original RV time series. The power p in the periodogram at period P is computed. K is then progressively increased until it reaches the fap level of the power of the peak down to 1%. It corresponds to the power of the 990th values in the power distribution determined above, after ranking them in increasing order. The amplitude is then converted into a planetary mass (using the stellar mass estimated in our previous study). We therefore obtain, for each period P, 12 planetary mass estimations above which we can rule out the presence of the planet with a 1% confidence level. The amplitude of the signal is phase-dependent, and thus we can choose either a highly conservative threshold by considering the worst phase and the largest amplitude (referred to as the conservative limit), or a statistical threshold by considering the median value (statistical limit). In this statistical study, we adopted the same assumption as Bonfils et al. (2013a) and exclusively used the statistical detection limits. This assumes that it is highly improbable for all systems to be observed under the most unfavourable phase conditions for detection.

This is done for periods of between 1.6 and 10 000 days and the 197 stars in our sample. As an example, Fig. 1 shows the statistical (dark turquoise) and conservative (yellow) detection limits for GJ 163 as a function of the period.

It is worth noting that this methodology incorporates the injection of sinusoidal functions, not Keplerian ones, into the RV series. According to Cumming & Dragomir (2010), this modelling remains robust up to an eccentricity of 0.5, suggesting that any signal with a higher eccentricity would be underestimated. However, given that 81.4% of the detected planets possess an eccentricity lower than 0.4⁵, we can consider this modelling as sufficiently conservative in this statistical study. Furthermore, this method is correct under the assumption that the noise on the data is white. We discuss the consequences of this in detail in Sect. 6.2.

3.2 Detection map

The linear-frequency sampling of the GLS periodogram makes the detection limits extremely erratic and noisy at short periods (as illustrated in the lower plot of Fig. 1). As detailed by Zechmeister & Kürster (2009), the frequency step corresponds to the inverse of the longest tested period, which in our case is the same for all targets: 10 000 days. This means that very high frequencies (corresponding to very short periods) are tested with inadequate sampling (e.g. 11 nights over a single season), leading to the oversampling of these short periods, especially compared to the uncertainties associated with the fitted planetary periods. To mitigate erratic rates, we smooth these limits over the period range of 1.6-10 days (400 points sampled on a log-uniform scale). We then generate the detection map for our dataset by stacking the statistical detection limits⁶. Figure 2 illustrates the distribution of detection limits of projected mass as a function of the probed orbital period. In this representation, the colour-coded scheme enhances the visualisation of iso-frequency curves (the fraction of stars in our sample around which a given planet could be detected). The grey area at the bottom of the map corresponds to the out-of-reach domain (below the lowest detection limit of 0.5%). This region encompasses planets with projected masses ranging from 0.6 to 1 M_⊕ for periods of fewer than ten days; this extends up to 20 M_⊕ for periods as long as 1000 days. Conversely, the white area in the upper left corner of the map corresponds to the domain where planets cannot be missed (beyond the highest detection limit) with a 99% confidence level. This region encompasses planets with projected masses ranging from 20 to 50 M_⊕ for periods shorter than ten days and up to 1000 M_⊕ for periods shorter than 1000 days. The iso-frequency contours on this map are interpreted as follows: if planets in the green zone exist around M dwarfs in this sample, we should expect to detect them in 30-40% of cases (around 60 to 80 stars among our 197) in this specific region.

Finally, this gives the local detectability metric for the mass-period map, which will is used to weigh each planet detection in the occurrence rate calculations and to characterise the accessibility of any size of domain in the next section.

Fig. 1 Conservative (yellow) and statistical (dark turquoise) detection limits of GJ 163 as function of projected planetary mass and periods, with the respective positions of the detections of GJ 163 b, c, and d. Upper panel shows the limits in a range of periods between 1.6 and 10 000 days, while the lower panel shows a zoomed-in view of these limits for periods between 1.6 and ten days.

Fig. 2 Iso-frequency map obtained by stacking the 197 statistical detection limits in projected masses as a function of the period. Iso-contours are highlighted for several percentage levels of the sample and describe the distribution of maximum mass missed. The grey part represents the unknown domain, below the lowest statistical detection limit (corresponding to 0.5%), and the white part is the domain beyond 100%.

4 Construction of the occurrence statistic

Planetary occurrence rates provide the expected count of planets around stars, considering both the detected planets and the detection limits. Consequently, a rate of 100% implies, for example, an expectation of 100 planets orbiting 100 stars, without specific knowledge of their distribution within the systems. We partition the map of projected planetary mass and orbital period (hereafter referred to as the mass-period map) into sub-domains where we compute these occurrence rates. In the subsequent sections of this paper, a domain jk therefore refers to a region on this map delimited by projected mass values between m_k and m_k+1 and period values between p_j and p_j+1.

4.1 Detections

In Bonfils et al. (2013a), for a given domain the occurrence rate was the ratio of the number of detections, ndet, to the number of accessible targets, N_eff. To take into account the variation in the detectability of a planet within a domain, in this work we adopted another method to compute the occurrence rate of a given domain: the ratio of the number of predicted planets in the domain to the total number of targets in the sample. In the jk domain of the mass-period map, we computed the local occurrence frequency f_jk as $${\mathrm{f}}_{\mathrm{j}\mathrm{k}}=\frac{{\tilde{\mathrm{n}}}_{\mathrm{j}\mathrm{k}}}{\mathrm{N}}=\frac{1}{\mathrm{N}}\times \sum _{\mathrm{i}=1}^{{\mathrm{n}}_{\mathrm{j}\mathrm{k}}}{{\alpha }_{({\mathrm{m}}_{\mathrm{i}};{\mathrm{p}}_{\mathrm{i}})}}^{-1}$$(1)

where α_{(mi; pi)} is the local detectability rate of the planet i, n_jk is the total count of detections in the jk domain and N is the total number of stars of the sample. ${\tilde{n}}_{jk}$ is therefore the weighted number of detections in the domain. As an illustration, a unique detection with a local detectability rate of about 50% on a sample of 100 stars gives an occurrence rate of 2%.

The two ratios give identical results under two circumstances: when detections are uniformly distributed within the domain, or when the domain’s detectability rate remains constant (which only occurs in the top left corner and the bottom right corner). Otherwise, rates obtained with the ratio used in Bonfils et al. (2013a) can be overestimated if detections are in the top of the domain, whereas the number of accessible targets is computed for the average mass and period in the range in the domain. Conversely, these rates can also be underestimated if detections are located in the least-massive-planet sections of the domain. As a result, occurrence rates are very sensitive to the position of the detected planet in the mass-period map. We also note that they can be over 100%, which implies local multiplicity.

4.2 Coefficient of multiplicity

In a general case, we compute the statistical bounds of the occurrence rates using a binomial law: we consider a detection as a success of the event. This law admits a maximum frequency of 100% for N planets detected around N stars. This is why multiplicity plays an important role in this work. The local coefficient multiplicity mult_jk is the average multiplicity of all the stars accessible (with and without detection) in the jk domain. This definition of local multiplicity in a domain combines observed multiplicity and probability of additional planets. We compute the multiplicity mult_{i jk} of each star in the jk domain as $$\mathrm{m}\mathrm{u}\mathrm{l}{\mathrm{t}}_{\mathrm{i}\mathrm{j}\mathrm{k}}=\mathrm{m}\mathrm{u}\mathrm{l}{\mathrm{t}}_{\mathrm{o}\mathrm{b}\mathrm{s}\mathrm{i}\mathrm{j}\mathrm{k}}+\overline{{\alpha }_{\mathit{\text{jk}}}}=\mathrm{m}\mathrm{u}\mathrm{l}{\mathrm{t}}_{\mathrm{o}\mathrm{b}\mathrm{s}\mathrm{i}\mathrm{j}\mathrm{k}}+1-<\alpha {>}_{\mathit{\text{jk}}}$$(2)

where mult_{obs i jk} is the observed multiplicity of the system i in the jk domain; < α >_jk is the average detectability rate in the domain, which is a probability of detection; and, therefore, $\overline{{\alpha }_{jk}}$ is the probability of a missed planet in a system in the considered jk domain. This coefficient gives a more complex multiplicity than the observed one (average number of planets in each system in the domain) because it takes into account the domain’s detectability limit. For a given system of one planet, if the local detectability rate is around 100% (as in the region corresponding to massive planets at a short period, for example), the weighted multiplicity is only one, because we know we cannot miss a detection with a confidence level beyond 99%. Conversely, if the local detectability rate of the domain is 60%, for example, there is a non-negligible probability of missing a second planet; therefore, the weighted multiplicity is higher than one for this system.

4.3 Statistical bounds

We compute statistical bounds using a binomial law of N pickings with replacements because events are not linked to each other and can only have two possible outcomes; i.e. detection or no detection. As the occurrence rate can reach 100% and beyond, we use the multiplicity coefficient to compute the statistical bounds. Returning to a given jk domain, with a mult multiplicity, an f occurrence frequency, and n detections, we compute the 1σ statistical bounds (lower f_low and upper f_up) as $${\int }_0^{f_{up}}\left({}_{N_{eff}\times mult}^n\right)\times f^n(1-f{)}^{(N_{eff}\times mult-n)}df=0.84,$$(3a) $${\int }_{f_{low}}^1\left({}_{N_{eff}\times mult}^n\right)\times f^n(1-f{)}^{(N_{eff}\times mult-n)}df=0.84,$$(3b)

where N_eff is the number of accessible stars in the domain used in Bonfils et al. (2013a). We obtain N_eff by multiplying the total number of stars within the sample by the average probability of detection in the domain (obtained using a log-uniform Monte Carlo sampling of 15 000 trials in the domain). This computation assumes a log-uniform probability distribution in both mass and period, as in Bonfils et al. (2013a). Finally, the greyed-out part of Fig. 2 is not included in the subsequent calculations.

Through this computation, we establish the 1σ confidence interval for the binomial probability distribution associated with achieving n successes in N_eff.mult trials with a f probability of success. In domains where multiplicity computed is below one, we consider only one trial per accessible target; i.e. mult parameter set to 1 in Eq. (3). Differently from previous studies, we consider multiplicity, enabling the possibility of drawing a star anew, irrespectively of the outcome of the previous draw. This adaptation stems from the fact that the binomial distribution implies independent random draws. However, it is important to acknowledge a methodological limitation here: within a given domain, the detection of a subsequent planet is not entirely independent of the initial planet detection. This limitation is discussed in Sect. 7. For domains without detection, we exclusively compute an upper limit through the construction of a binomial distribution assuming zero successes within a number of trials equivalent to the number of accessible stars in the domain N_eff.

Fig. 3 Planetary occurrence statistic. Colours in the background refer to the detection limit map, from the non-reachable region (grey) to the unmissable one (white). The grey points represent the considered detections, while the red points represent the candidates. The grey horizontal and vertical lines delineate the probed domains. The computed occurrence rates (along with their corresponding statistical bounds) are indicated at the centre of each domain.

4.4 Planetary system occurrence statistic

Beyond the number of expected planets within a sample of 100 stars, this study also examines the occurrence of planetary systems, namely the number of expected planetary systems per 100 examined stars. It is crucial to ascertain whether an occurrence rate of 100 planets among 100 stars corresponds to 100 single systems or 50 multiple systems, for example. For this computation, we count the number of systems rather than the number of individual planets. Within each domain, a multiple system is counted only once, and the local weighting attributed to the system is the individual weighting of the first planet (easiest to detect) composing the system. The occurrence rate of planetary systems within a domain is therefore ratio similar to the planetary occurrence. Consequently, the statistical bounds are determined analogously to the previous case, but with a number of detections corresponding to the number of the systems and a multiplicity value of one. This computation is a first estimation of the system occurrence rates, according to the method of computation we chose.

5 Results

In this section, we present the planetary and system occurrence rates obtained through the method described above. Additionally, we provide the rates derived from sub-samples, which represent distinct populations aimed at better understanding the relationship between planetary populations and the physical parameters of the central star.

5.1 Planetary occurrence rates

To initiate the presentation of results, Fig. 3 presents the rates from the entire sample of 197 M dwarfs as a function of projected planetary mass and orbital period. Within this visual representation, we overlay the detectability map and the demarcations of the previously discussed domains, in addition to the occurrence rates and their associated statistical bounds. For this first representation, we categorise the planetary masses into four ranges: terrestrial planets of very low mass (m sin i < 3.16 M_⊕), intermediate-mass planets ranging between super-Earths and Neptunes (3.2 < m sin i < 31.6 M_⊕), intermediate-mass volatile planets (31.6 < m sin i < 1 M_J), and Jovian planets (m sin i > 1 M_J). A projected mass of 3162 Earth masses corresponds statistically to an actual mass of approximately 4700 Earth masses, which exceeds the 13-Jupiter-mass threshold considered as the boundary between massive planets and brown dwarfs by the IAU (deuterium fusion critical mass). This is why we chose not to explore the higher mass range and propose this segmentation.

In the lower mass range (m sin i < 3.16 M_⊕), we obtain an occurrence frequency of $120{\mathrm{\%}}_{-23.7}^{+24.8}$ for periods shorter than ten days. Considering that the longer period domains are poorly covered by the detection map, we do not compute the frequencies for these planets beyond ten days. In the higher mass range, including super-Earths and Neptunes, we obtain a rate of $14.3{\mathrm{\%}}_{-3.1}^{+4.8}$ for periods shorter than ten days, $27.9{\mathrm{\%}}_{-5.8}^{+8.9}$ for those between ten and 100 days, $17.7{\mathrm{\%}}_{-4.1}^{+6.1}$ for those between 100 and 1000 days, and, finally, an upper limit of 7.2% beyond 1000 days, as we have no detections in this less explored domain.

Table 1 presents both planetary occurrence and system frequencies within these two specified mass ranges at short orbital periods. Within the lower mass range, for periods shorter than ten days the expected $120{\mathrm{\%}}_{-24}^{+25}$ planetary occurrence frequency is concentrated around $80{\mathrm{\%}}_{-18}^{+16}$ of the targets accessible in the domain. These frequencies could indicate that 20% of the systems investigated with HARPS do not host planets of extremely low masses with short orbital periods. Conversely, the remaining systems would occasionally exhibit the presence of multiple planets, which necessarily implies the existence of distinct and varied planetary systems within the surveyed sample. This hypothesis is reinforced by the multiplicity of 1.6 computed in this domain. Within the higher mass ranges, the expected quantity of planets is fully compatible with the expected number of systems, thus suggesting no multiplicity in the other domains.

In the intermediate volatile-planet mass range (31.6 < m sin i < 1 M_J), we obtain an upper limit of 0.5% for periods shorter than ten days because we have no detections in this domain, even though it is well constrained by detection limits. This absence of planets in this mass range over short periods is not a bias of the study or an observational bias, but rather an absence because such planets could not have been missed during such a period. Between ten and 100 days, we obtain a rate of $0.5{\mathrm{\%}}_{-0.1}^{+1.2}$. The rate is $1.4{\mathrm{\%}}_{-0.4}^{+1.9}$ between 100 and 1000 days and $1.5{\mathrm{\%}}_{-1.5}^{+1.8}$ between 1000 and 10 000 days. In the mass range of Jovian planets (m sin i > 1 M_J), we again obtain an upper limit for periods shorter than ten days because we have no detection in this domain, despite it being well constrained. We obtain a rate of 0.5%${}_{-0.1}^{+1.2}$ between ten and 100 days, 0.6%${}_{-0.2}^{+1.2}$ between 100 and 1000 days, and 3.0%${}_{-0.8}^{+2.0}$ between 1000 and 10 000 days.

Table 2 lists all planetary occurrence rates f_p computed and presented in Fig. 3. It also provides the number of detections n_d and the number of accessible targets in the domain N_eff. In addition, the parameter N_eff′ provides a clearer representation of the difference between the calculation based on the average weighting of the domain and the approach used in this study, which applies local weighting to each detection. The two calculations are equivalent if detections are homogeneously distributed in the domain when N_eff is equal to N_eff′. When N_eff > N_eff′, the detections are concentrated in the upper part of the surveyed domain, and, conversely, when N_eff < N_eff′ the detections are located in the lower part of the domain.

Fig. 4 Planetary occurrence statistics obtained from sub-sample of M-dwarfs with M* > 0.35 M⊙.

Fig. 5 Planetary occurrence statistics obtained on sub-sample of M-dwarfs with M* < 0.35 M⊙.

5.2 Dependence on the central mass

In this section, we examine the correlation between the planetary population and the mass of the central star. For this purpose, we split the sample according to the mass of the central star: below and above 0.35 solar masses, limit between partially and fully convective (Chabrier & Baraffe 1997). Therefore, we analyse and compare two distinct samples: the first comprises 145 very massive M dwarfs with a median stellar mass of 0.48 M_⊙; while the second consists of 52 less massive M dwarfs with a median stellar mass of 0.26 M_⊙.

To compare the two planetary populations, we calculate the rates within the same four mass ranges: terrestrial planets of very low masses (m sin i < 3.16 M_⊕), intermediate-mass planets ranging between super-Earths and Neptunes (3.2 < m sin i < 31.6 M_⊕), intermediate-mass volatile planets (31.6 < m sin i < 1 M_J), and massive Jovian planets (m sin i > 1 M_J). Figure 4 displays the occurrence rates obtained for the subsample of more massive M dwarfs (M_* > 0.35 M_⊙), and Fig. 5 shows those obtained for the sub-sample consisting of the less massive ones (M_* < 0.35 M_⊙) and all the values are presented in Table 3 for ease of comparison.

Regarding the smallest planetary mass range, the majority of detected planets orbit less massive M dwarfs, resulting in an occurrence rate of $116{\mathrm{\%}}_{-23}^{+20}$ for periods shorter than ten days for those stars. This rate is based on a total of 11 detections distributed among eight different systems in a domain where only 11 stars are accessible. Concerning the more massive M dwarfs, we count a unique detection (GJ 393 b) only; this is in a domain poorly covered by the detectability map, and consequently we did not calculate the occurrence rate in this mass range domain for the sample of more massive M dwarfs. However, this unique detection occurs in a domain where about 10% of the systems are accessible, suggesting a lower occurrence rate compared to less massive stars. This reinforces the hypothesis that the occurrence frequency of these small terrestrial planets increases as stellar mass decreases (e.g. Ment & Charbonneau 2023).

In the intermediate planetary-mass range (Super-Earth to Neptunes), the occurrence rates are fairly homogeneous and compatible for periods of up to 1000 days for the two subsamples of M dwarfs. Beyond 1000 days, no planets of this mass range have been detected in the sub-samples, resulting in two upper limits of 4.8% and 7.4% for the more massive and less massive M dwarfs, respectively; however, this domain is poorly covered for both stellar population.

Finally, for massive planets (from Neptune to the brown dwarfs limit), there is a hint of a dichotomy between the two planetary populations. Specifically, the only two giant planets detected during short periods (<100 days) orbit low-mass M dwarfs (both orbiting the same system GJ 876), while all giant planets during long periods orbit the more massive M dwarfs (upper right corner of Fig. 4). This difference in the population of giant planets at large separations is discussed in Sect. 7.

5.3 Comparison and evolution

We now compare the results obtained using this new calculation method and new HARPS sample with key studies. To facilitate a meaningful comparison of these rates with our study, we adopt the same mass segmentation and represent the occurrence rates obtained by the studies in Fig. B.1. We report them in Table A.4 for each probed domain. First, we compare with the occurrence rates obtained by Bonfils et al. (2013a) ten years ago on a sample of M dwarfs observed with HARPS. Next, we compare our results with the occurrence rates obtained more recently by Sabotta et al. (2021) on a sample of M dwarfs observed with CARMENES. Finally, we compare the results obtained on the sub-sample of the most massive M dwarfs to the recent study of Pinamonti et al. (2022).

5.3.1 From 2013 to 2023

The detection rates reported by Bonfils et al. (2013a) are based on more than 2000 RV measurements from 102 targets and resulted in 20 detections, including 11 new ones. None of these detections are massive planets (m sin i > 100 M_⊕), leading to an upper limit of 1% of hot Jupiters (period < 10 d). As this study is based on a maximum temporal coverage of seven years, the occurrence rates of massive objects with very long periods were not extensively explored (constrained by detection limits in this regime).

Concerning planets in the smallest range of mass that could be probed (1 < m sin i < 10 M_⊕), Bonfils et al. (2013a) predicts occurrence rates of 36%${}_{-10}^{+24}$ for periods of between 1 and 10 days and 52%${}_{-16}^{+50}$ for periods of between ten and 100 days. In this smallest mass range, we computed a frequency of 89%${}_{-15}^{+18}$ incompatible at a 1-σ level for very short periods (p < 10 days). In this domain, there were five detections from 14 accessible stars (N_eff) in the previous study; while in our study, we have 18 detections out of 48 accessible stars. Within this domain, the detectability rate ranges from 0.5% to 70% in our study (see Fig. 2) and 1% and 50% in that of Bonfils et al. (2013a). Consequently, the incompatibility is due to the new computation method we use in this study, and it highlights the consequences of the assumptions made in both methods, as described in Sect. 4.1.

This situation is also observed for the mass range between 10 and 100 M_⊕ and periods of between 100 and 1000 days. The upper bound obtained in Bonfils et al. (2013a) is incompatible with the rate obtained in this study, based on eight detections inhomogeneously distributed in a domain where detectability rates vary significantly. In the other domains, the rates obtained by both studies are in agreement (within a 1-σ confidence interval), but enlarging the sample size improves uncertainties. Importantly, this new study is also built upon a larger number of detections, ultimately allowing for precise rates to be obtained in five domains where only an upper-bound estimate was previously possible without detections. Finally, this mass segmentation allows us to probe between 10³ and 10⁴ M_⊕; however, it should be noted that we have not included objects larger than 13 M_J (i.e. approximately 4100 M_⊕) in our analysis. Although this mass range extends beyond this limit, any potential detections therein are not included, and thus the rates obtained are underestimated.

5.3.2 CARMENES occurrence rates obtained in RV

Sabotta et al. (2021) (hereafter S21) built their study on the 71 M dwarfs observed over at least 50 nights from the 329 M dwarfs followed in the CARMENES survey. The occurrence rates are based on 27 detections corresponding to 21 planetary systems. As a main result, they computed more than one small planet (1 < m sin i < 10 M_⊕) for periods shorter than 100 days ($132{\mathrm{\%}}_{-32}^{+33}$). As can be seen in Table A.4, the rates obtained from our sample are all consistent (within 1 - σ) with those from the CARMENES study, as well as with the upper bounds in the detection-lacking domains. These occurrence rates were later revised by Ribas et al. (2023) using a larger sample of 238 M dwarfs observed with CARMENES. The updated rates remain consistent with our results across all but two domains (lowest planetary-mass range and longest periods). The origin of this discrepancy is identical to the one discussed in detail in the following paragraph.

Regarding the stellar-mass dependence, we did indeed conduct our calculations in the same domains to facilitate comparison, although we did not include this in the paper. S21 chose to set the stellar mass cut-off at 0.34 M_⊙ - a value very close to our choice of 0.35 M_⊙ - considering the measurement uncertainties for the stellar masses in our samples.

For low-mass M dwarfs (< 0.34 M_⊙ in S21, < 0.35 M_⊙ in our study), our occurrence rates are consistent within 1 σ, except for the rate obtained for long periods (100 < P < 1000 days) with high-mass planets (100 < m sin i < 1000 M_⊕). In this range, S21 reports a rate of ${16}_{-9}^{+15}\mathrm{\%}$, with one detection among the 23 stars in their sub-sample, while we provide an upper limit of 2% since no planets were detected in this range among the 56 stars in our sample. This discrepancy stems from the fact that the method used in S21 to compute the occurrence differs from the method we propose. The computation of the occurrences depends on the detectability rates in the domain probed. However, such long periods are poorly constrained due to CARMENES’ temporal coverage. The local detectability rate of their detection in this range (around 60%) does not correspond to the effective number of stars (Neff) used to compute the rate (around 6.25%). In their computation, only the number of detections in the domain is used to validate the frequencies, and not their location and local rate in the domain. Consequently, the occurrence rate reported by S21 can be overestimated as it is dominated by the very low detectability rates for long periods. By applying the local detectability rate to their detection, we obtain a rate consistent with our upper limit.

For more massive M dwarfs, our occurrence frequencies are consistent within 1 σ for both short periods (1 < P < 10 days) and long periods (100 < P < 1000 days) across all three planetary mass ranges. However, for intermediate periods (10 < P < 100 days), our rates are inconsistent for various reasons. In the low-mass regime (1 < m sin i < 10 M_⊕), S21 reports ${210}_{-81}^{+113}\mathrm{\%}$ with four detections, while we report ${25}_{-7}^{+12}\mathrm{\%}$ with six detections. Once again, the four detections from S21 are located in the upper left corner of the probed domain. The discrepancy between our rates is therefore due to the computation method once again. The same applies to the intermediate-mass range (10 < m sin i < 100 M_⊕), where the detections are again not uniformly distributed. Finally, in the high-mass regime (100 < m sin i < 1000 M_⊕), S21 reports two detections, yielding an occurrence rate of $6_{-1}^{+4}\mathrm{\%}$, while we provide an upper limit of 1% due to the absence of detection in this range. The two detections from S21 correspond to GJ 876 b and c (Karmn J22532-142), for which the stellar mass reported in Table A1 of S21 is 0.327 M_⊙. This mass is below 0.34 M_⊙, meaning the star was placed in the wrong sub-sample in S21. Without this error, our rates are fully consistent. GJ 876 is also part of our sample, and both planets are indeed accounted for in our low-mass-star sub-sample. In conclusion, aside from this error, the discrepancies between the occurrence rates in our two studies are due to differences in the calculation methods.

5.3.3 HARPS-N occurrence rates obtained in RV

Pinamonti et al. (2022) (hereafter P22) established occurrence frequencies from a sample of 56 massive M dwarfs (M0-M3) observed with HARPS-N. In the lowest mass domain (1 < m sin i < 10 M_⊕), for the shortest periods (1 < P < 10 d) they obtain ${10.3}_{-3.3}^{+8.4}\mathrm{\%}$. This is in agreement with the frequency of ${11.8}_{-3.3}^{+7.3}\mathrm{\%}$ we obtained with the sample of massive M dwarfs for the same mass-period ranges. Concerning the longest periods (10 < P < 100 d), they obtained ${85}_{-19}^{+5}\mathrm{\%},$, and we obtained ${25}_{-7}^{+12}\mathrm{\%}$ % in this study. Once again, our frequencies are not in agreement because the frequencies computed by P22 depends on the average value of the detectability in the domain. Their detections are spread inhomogenously in the domain: all in the upper left corner, where the local rate varies between 0.3 and 0.75; this is far from the average value taken to compute the frequency in that domain. Consequently, with this method, and as explained previously, the frequencies obtained can be overestimated. In conclusion, the incompatibility between our results in that domain of periods is not due to the sample of targets, but to the method of computation.

5.4 Dependence on the number of measurements

This study relies on the sample established in paper I, wherein we decided to consider only M dwarfs observed for at least ten nights with HARPS. Comparing the rates obtained by Sabotta et al. (2021) on the sample of M dwarfs observed for more than 50 nights with CARMENES raises the question of the impact of the initial sample choice. Here, we present the occurrence rates obtained for different sub-samples of our initial sample by selecting targets based on the number of observations. Figure 6 displays the rates obtained in the previously probed mass and period domains for four samples (from left to right): the 197 M dwarfs observed for at least ten nights, the 87 observed for at least 30 nights, the 57 observed for at least 50 nights, and the 36 observed for at least 100 nights.

This selection reduces the total number of targets in the probed sample; thereby, it alters the weighting of each detection. Additionally, since the highest detection limits often come from the less extensive RV series, increasing the minimum number of nights in the series purges the detectability map from the upper end. Specifically, the 90% detection limit for the sample of extensively observed M dwarfs (at least 100 nights) corresponds to a much lower mass than that of the M dwarfs in the original sample (observed for at least 10 nights). These two effects more or less compensate each other, allowing us to distinguish three cases in Fig. 6. When the total number of stars dominates the ratio, the occurrence rates increase with the minimum number of nights (most of the probed domains). When the increase in weighting compensates for the decrease in the total number, the occurrence rates plateau, for example in the domain of very low masses at short periods (bottom left corner). Finally, when the number of detections decreases, the rates fall sharply, as in the domains of long-period Jovian-mass planets (top right corner). In these domains, the number of detections decreases when we only consider the most intensively observed M dwarfs, indicating that systems with such massive cold planets have not been intensively monitored, but rather excluded from the HARPS programs. This is a major bias, as it prevents us from concluding whether there is a correlation, anti-correlation, or non-correlation of the internal and external planetary populations, even though this question is now crucial and has not yet been answered (i.e. Bonomo et al. 2023; Bryan & Lee 2024).

Fig. 6 Comparison of occurrence rates obtained as function of minimal number of measurements to consider in the mass-period domains probed.

6 Occurrence statistics in temperate zone

The range of orbital distances, where the conditions permit the presence of liquid water on a planetary surface, is commonly referred to as the habitable zone (HZ; Huang 1959; Hart 1978). Indeed, this study does not address habitability conditions, which are acknowledged to be significantly more complex than the concept discussed herein. To ensure a clear and unbiased interpretation of this concept, we use a more neutral term, such as temperate zone (TZ). As the TZ mostly depends on the effective temperature of each star, it does not correspond to a fixed domain in a period. To obtain an occurrence rate in the TZs of 197 targets, we first convert the detection limits from period-dependent to insolation-dependent, as in Pinamonti et al. (2022). We then normalise all the detection limits on the inner and outer limits of the TZ. Finally, we compute the new occurrence rates in the validity domain of the model.

6.1 From periods to insolation

First, we compute the detection limits from period to distance using Kepler’s third law, and then we convert it into insolation: $$\mathrm{S}\mathrm{e}\mathrm{f}\mathrm{f}=(\frac{\mathrm{L}}{{\mathrm{d}}^2})$$(4)

where the distance d is converted to effective insolation, S_eff, with the luminosity of the star, L (in L). The luminosity is primarily derived from the TESS catalogue⁷. For the few missing cases, luminosity is computed using the Boltzmann law, with a radius linearly derived from the mass estimated in paper I. Subsequently, we can transform the detection limits into insolation, as illustrated in Fig. 7, which represents the statistical limit of the projected planetary mass as a function of received insolation of GJ 163.

Subsequently, we calculate the values of the inner and outer boundaries (Seff_i) of the TZ as defined by Kopparapu et al. (2014) for 1 M_⊕ (represented by the vertical bars in Fig. 7). To achieve this, we use the i-th set of parameters corresponding to the i-th boundary in the following relation: $$\mathrm{S}\mathrm{e}\mathrm{f}{\mathrm{f}}_{\mathrm{i}}=\mathrm{S}\mathrm{e}\mathrm{f}{\mathrm{f}}_{\mathrm{i}\odot }+{\mathrm{a}}_{\mathrm{i}}\mathrm{T}+{\mathrm{b}}_{\mathrm{i}}{\mathrm{T}}^2+{\mathrm{c}}_{\mathrm{i}}{\mathrm{T}}^3+{\mathrm{d}}_{\mathrm{i}}{\mathrm{T}}^4$$(5)

where T represents the difference between the effective temperature of the star and that of the Sun. The four sets of parameters used are compiled in Table A.5. Following Kopparapu et al. (2014), we opt for the most conservative boundaries of TZ: runaway greenhouse (inner) and maximum greenhouse (outer). After performing this calculation for the 197 targets in our sample, we normalsze them by aligning the inner and outer boundaries. Consequently, it becomes possible to overlay the 197 detection limits and construct a detectability map similar to the one calculated previously, but this time based on insolation rather than the period. However, it is essential to note that this model is only valid for terrestrial planets, and thus the planetary mass to consider falls between 0.1 and 10 M_⊕. Therefore, in the following, we calculate the occurrence rate by only counting and weighting the detected planets in the probed TZ domain. Finally, we obtain the statistical bounds using the same computation as for the planetary occurrence.

Fig. 7 GJ 163 statistical mass limit as function of insolation. The vertical bars represent the inner (orange: recent Venus; yellow: runaway greenhouse) and outer (green: maximum greenhouse; blue: early Mars) boundaries of the HZ defined in Kopparapu et al. (2013).

6.2 Results

Taking the most conservative boundaries of the TZ, we identify four planets within this zone: GJ 551 b (known as Proxima b), GJ 3053 b, GJ 3293 d, and GJ 3323 c. This resulted in an occurrence rate of 45.3%${}_{-16}^{+20}$, as illustrated in Fig. 8. This rate is in agreement with other estimations of the planetary occurrence in HZ obtained using the RV method (e.g. 41% Bonfils et al. 2013a). Our occurrence rate is also consistent to 1σ with those obtained around M dwarfs using the transit method - with a rate of 15.82%${}_{-6.54}^{+16.60}$ obtained by Dressing & Charbonneau (2015) -but inconsistent with the 8.58%${}_{-8.22}^{+17.94}$ obtained more recently by Bergsten et al. (2023) on Kepler surveys using models also based on insolation.

It is crucial to emphasise that our rate is derived from a homogeneous analysis of the 197 RV time series, which did not recover 100% of the planets detected and published to date; therefore, some planets might be missing from our reference sample of detections. However, it should be noted that if a planet has not been detected, its signal has not been subtracted from the RV time series; consequently, this absence is considered in the calculation of the detection limits of this study, thereby affecting the computation of the rate and statistical bounds. In conclusion, the discrepancy between our rate obtained by RV and those previously obtained by transit is mostly due to the non-determination of the radius. Indeed, working on minimum projected mass, the most massive planets determined in the TZ could be sub-Neptunes rather than tellurics, rendering the applied model invalid and our rate overestimated.

Fig. 8 Conservative occurrence rate and planets in TZ.

6.3 Planets in TZ

6.3.1 GJ 551 b: Proxima b

Although Proxima b is currently considered to be within the HZ (e.g. Faria et al. 2022), it is crucial to note that this determination relies heavily on the effective temperature value. Without a measurement from TESS, we selected an effective temperature from the five available values on the CDS. When adopting the most recent value of 3306 K (Maldonado et al. 2020), Proxima b falls within the inner zone, where it receives excessive insolation. The other four references suggest values between 2800 and 3000 K, which are more consistent with the star’s spectral type (M5.5V). Consequently, we selected the second most recent value of 2810 K (Kuznetsov et al. 2019) for this analysis. This case clearly demonstrates the high sensitivity of our study to the effective temperature.

6.3.2 GJ 3053 b and GJ 3293 d

These two planets, despite being within the TZ, possess projected masses near the limit of the model’s applicability (Kopparapu et al. 2014). In this mass range, the determination of the planet’s radius is crucial to distinguishing whether they are super-Earths or sub-Neptunes.

6.3.3 GJ 3323 c

This planet is on the conservative outer boundary of the TZ (maximum greenhouse). Consequently, its position is highly sensitive to the effective-temperature value adopted from the literature (18 references ranging from 3000 to 3300 K), potentially placing it outside the TZ.

6.3.4 GJ 667C c

Although this star is considered out of the TZ in this study, it is located very close to the conservative inner limit (runaway greenhouse). The selection of the effective temperature is again crucial: based on the TESS catalogue, we adopted a value of 3650 K. However, it is noteworthy that the planet falls on the inner boundary if using the temperature from Gaidos & Mann (2014) 3472 K, and within the TZ when using that of Santos et al. (2013) 3351 K. In addition, recent studies using planetary climate models confirm the presence of GJ 667C c in the HZ or extremely close to the internal limit (see e.g. Fig 8 of Turbet et al. 2023). In the latter case, the occurrence rate of planets in the TZ is 49.6%${}_{-19}^{+13}$.

7 Discussion

7.1 Intrinsic limitation of the study framework

In this section, we discuss the impact of assumptions chosen to compute this new occurrence rate.

7.1.1 Limits of correction of activity

As a first limitation of this study, it is important to point out that 23 series have been corrected for sinusoidal signals attributed to the impact of stellar activity, especially stellar rotation in Mignon et al. (2024). Subtracting a Keplerian signal does not clearly remove the contribution of stellar activity, especially in the case of long time series; in this case, a phase change in the signal over time and a variation in its amplitude have not been taken into consideration. Given that we considered all series of at least ten nights of observation, it was not feasible to implement a more complex method, such as Gaussian processes. However, in our sample, 19 stars have more than 150 measurements, and among them only three (GJ 229, GJ 667C, and GJ 393) exhibited RV signals at varying periods across different seasons. For these three stars, we corrected for two sinusoidal signals with close periods, and we may have introduced additional signal into the RV series. In addition, among the 19 stars, some detectionlimit masses might also be overestimated because RV time series have not been completely cleaned from stellar signals. However, since only three out of 197 stars are affected by this inappropriate correction, we believe that the results concerning the detection limits of the 197 stars are not significantly compromised.

7.1.2 Limits of planetary detection method

Out of the 31 planetary systems detected, eight are incomplete when compared to the NASA Exoplanet Archive. As an example, GJ 180 c and d and GJ433 d were detected with complementary data (Tuomi et al. (2014) using UVES data and Feng et al. (2020) using VLT2 data). We verified that all these planets lie below our detection limits, and we list them in Table A.3 the reference paper and provide the method and the data of the discovery.

Regarding the impact of these missing planets on the occurrence rates, we cannot include them in our detectability map as their signals were not removed from the RV series used to determine detection limits. These planets could potentially be detected with more RV data or by using more advanced methods. With well-sampled RV series, we could compute much lower detection limits, which would likely increase the number of detected planets. However, in the context of our study, this would mean focusing only on the most observed targets. As shown in Fig. 6, such targets are those for which we have detected planets. As a result, this would create a highly biased sample in the current state of the data.

In conclusion, while the method we chose does not recover all previously published planets, it allows us to analyse and process a large number of RV series. We consider this to be the least biassed method possible given the current dataset. The missed planets are all below our detection limits and are therefore accounted for in our statistical boundaries.

7.1.3 Limits of the binomial model

In this study, we compute statistical boundaries using a binomial distribution constructed within the surveyed domain, where each draw results in either detection or non-detection. As explained in Sect. 4, the distribution reaches its limit in domains with abundant multiple systems. Computing an average multiplicity in the proposed domain allows us to artificially increase the number of draws to account for this multiplicity. Specifically, the conducted experiment involves drawing each star mult times in a draw with replacement, based on the assumption of independent draws. However, this assumption can be questioned on two distinct aspects. Firstly, at the level of multiple system architecture, the presence of a planet with specific mass and period does not necessarily allow the presence of another planet with the same parameter set. During the random draw, it may be preferable to narrow down the surveyed zone of the domain by excluding parameter sets prohibited by the first planet, for example.

Secondly, the discovery of one planet has an impact on the discovery of a second. Indeed, observations are not conducted blindly, and the RV time series that are better sampled are those selected and strengthened after an initial discovery. Our set of RV series is therefore biased by an iterative process of searching for planets in the same system. While the assumption of independent draws is probably not the optimal approach, manipulating the size of the surveyed domains can render these effects more or less marginal.

7.1.4 Limits of the new computation of frequencies

By choosing to weigh each detection by its local detectability rate, rather than the average rate over the entire probed domain, we assign significant importance to the planetary parameters (the best fits presented in Table A.2. Regarding period uncertainties, the smoothing we applied ensures a uniform rate, as the precision of the map is of the same order of magnitude as the typical error (which is very small in RV measurements).

For uncertainties of the fitted masses, we performed the following calculation to quantify their impact. We calculated occurrence rates based on 1000 different detection sets, with each set being a random realisation of masses assuming Gaussian errors. We then computed the average value and standard deviation of the resulting occurrence-rate distribution. In the vast majority of the probed domains, the standard deviations are much smaller than the statistical bounds we calculated (from 0.01 to a few percent), demonstrating that the uncertainties in the masses are not the dominant factor in the occurrence-rate calculations; rather, the statistical bounds are. However, it should be noted that for lower mass planets, which are closer to the detection limits, the impact of the uncertainty becomes more significant (22%).

Finally, it is important to recall that the previous calculation, based on the average detectability rate over the entire probed domain, would imply error bars on the masses of the order of the domain size, which is not realistic for such large parameter spaces. In conclusion, by using the local detectability rate for each detection, we more accurately capture the distribution of detections within the domain.

7.1.5 Limits on the occurrence rates at long-term periods

This study is based on the sample of RV series compiled in paper I where long-term linear and quadratic trends were identified and subtracted from the series. Some series were corrected for potential signatures of massive companions whose orbits are still too poorly constrained to be characterised, leading to the determination of a minimum period and a minimum mass. However, the detection limits were calculated based on the residuals of these series, resulting in lowered limits for very long periods (beyond 1000 days). Consequently, the number of detections during these very long periods and the obtained occurrence rate are possibly underestimated. It is noteworthy that the vast majority of the massive candidates resulting from the previous study orbit the most massive M dwarfs in our sample. In conclusion, potentially underestimating the occurrence rate of massive planets at large separations does not fundamentally challenge the observed dichotomy between the least massive and most massive stars concerning this planetary population, it may even tend to reinforce it.

7.1.6 Limits of the TZ domain model used

The definition of the TZ we used has been computed for small terrestrial planets orbiting various stellar types and considering three different planetary masses: 0.1 M_⊕, 1M_⊕, and 10 M_⊕ (Kopparapu et al. 2014). The stellar type is defined through the effective temperature of the star, from 2600 K to 7200 K. Outside this range, the inner edge of the TZ is no longer valid. As we do not have access to the real masses of the planets or the composition of their atmosphere in our sample, we used the TZ estimated for one Earth mass assuming terrestrial composition. This assumption is commonly used in the literature to compensate the lack of data. Kopparapu et al. (2014) proposed different inner limits: the ’recent Venus’, ’runaway greenhouse’, and ’moist greenhouse’. In this work, we used the runaway greenhouse, which is the most conservative definition. Moreover, it is interesting to notice that an Earth-like planet (in terms of mass and atmospheric composition) is supposed to be tidally locked by assuming a 4.5 Gyr tidal-locking timescale for stars cooler than 4400 K (Kopparapu et al. 2014). In such a scenario, the global circulation in the atmosphere will be drastically different, and the inner limit of the TZ it will be closer to the star (e.g. Yang et al. 2013; Kopparapu et al. 2017). This configuration requires 3D modeling to accurately estimate the inner edge of the TZ, but no generic equation has been proposed so far to cover this regime.

7.2 Interpretation of the planetary occurrence around M dwarfs

7.2.1 Lack of hot Jupiters

With no detection of hot Jupiters among the 197 targets, we established an upper limit of 0.5% on their occurrence rate around M dwarfs. This updated upper limit represents the lowest within RV studies around M dwarfs, which is in agreement with prior studies: <1.27% in Endl et al. (2006), <1% in Bonfils et al. (2013a), and <3% in Sabotta et al. (2021) (see Fig. B.1). This new upper limit accentuates the apparent lack of hot Jupiters around M dwarfs compared to solar-type stars predicted in Endl et al. (2006). Indeed, occurrence rates of hot Jupiters detected via RV around FGK-type stars obtained so far are slightly higher: 1.5% in Cumming et al. (2008), 1.2 ± 0.4 in Wright et al. (2012), and more recently 0.84% in Wittenmyer et al. (2020). Without any detection of a hot Jupiter around an M dwarf in our probed sample of 197 M dwarfs, this upper limit only depends on the total volume of the sample. Thanks to the extensive volumes surveyed and the high sensitivity to hot Jupiters, the transit technique emerges as the most adept in refining the population contrast between M dwarfs and FGK stars. Notably, a recent study by Gan et al. (2023) based on the TESS survey of M dwarfs reported a rate of 0.27 ± 0.09%; this is slightly lower than the rates of ${0.43}_{-0.06}^{+0.07}\mathrm{\%}$ by Masuda & Winn (2017) and ${0.57}_{-0.12}^{+0.14}\mathrm{\%}$ by Petigura et al. (2018), which were both based on FGK-type samples using the transit method. Although still consistent within 2 σ, these frequencies obtained on transit surveys will need to be refined to account for the observed lack of hot Jupiters around M dwarfs as the RV method alone will not provide a large enough sample of M dwarfs to establish statistically significant conclusions on this matter.

7.2.2 Dependence of planetary population on the central mass

It is of particular interest to probe the dependency of planetary occurrence with stellar mass as this dimension can reveal new constraints on planet formation. The most accessible, first-order analysis involves the occurrence of planets as a function of stellar mass. From transit missions, a decreasing occurrence with increasing stellar mass was found (Dressing & Charbonneau 2013; Mulders et al. 2015; Dressing & Charbonneau 2015; Hsu et al. 2020). Recently, based on updated stellar parameters, the conclusion was questioned again (Bergsten et al. 2023). Although based on a smaller sample of stars, our results based on the radial velocity method can help to complement these findings.

There is no significant difference in planetary occurrence for the planets with m sin i from 3.16 to 31.6 M_⊕ at each respective period range. The lowest mass bin (Figs. 3, 4, and 5) with planets below 3.16 M_⊕ is inaccessible for the more massive stars. At larger planetary masses (<100 M_⊕), there is evidence that more planets exist around low-mass stars, but the number of planets is still limited.

The occurrence for FGK stars based on the RV method can be estimated from the Mayor et al. (2011) sample using EPOS (Mulders et al. 2018) and a broken power law in orbital period and mass space. In order not to be influenced by a potentially different giant planet population, we constrained the fit to the planets with m sin i < 100 M_⊕ and investigated the super-Earth-to-sub-Neptune mass bin ((3.2 < m sin i < 31.6 M_⊕)) shown in Fig. 9. Since this is a regime that is challenging to assess around solar-mass stars, useful constraints can only be obtained for planets with P < 100 days. In the innermost period bin (1.5 < P < 10 days), we find an occurrence of ${1.9}_{-1.5}^{+4.9}\mathrm{\%},$ which is significantly smaller than around M stars as reported here. This observation was also made by Sabotta et al. (2021) and was discussed in Schlecker et al. (2022). They hypothesised that a migration trap that depends on the disc temperature (and therefore also on stellar properties) could produce this outcome.

For the bin extending from 10 < P < 100 days, the EPOS analysis results in an occurrence of ${14}_{-10.7}^{+29.9}$, which is, given the large uncertainty, in agreement with the same mass-period region around M stars. Given that transit observations mainly probe the innermost orbital periods, the decrease in planet occurrence with increasing stellar mass might mainly be a close-orbit effect, and we cannot provide evidence with RV surveys that it also holds at larger orbital periods.

Fig. 9 Occurrence of low-mass planets around FGK stars based on Mayor et al. (2011) and obtained using EPOS (Mulders et al. 2018). The marked bins are the ones also used in Figs. 3, 4, and 5, and occurrence is indicated in percent and colour-coded in logarithmic mass against period space.

7.2.3 Comparison with planetary-population synthesis

Theoretical models have so far failed to comprehensively explain the correct stellar-mass-dependent planet distribution. For the aforementioned overall occurrence trend with stellar mass, Burn et al. (2021) predicted a slight increase of occurrence from FGK to early M before a decrease due to solid-mass limitation. The reason for the lower occurrence around FGK stars is the more common presence of massive giant planets and their gravitational perturbations of lower mass planets. Further effects are growth to larger mass bins and differences in minimum spacing between two planets if they become more massive due to their more extended Hill sphere. Another effect could be the halting of the pebble flux to the inner system, which is not included in that work (Lambrechts et al. 2019; Mulders et al. 2021). However, if the decrease in planetary occurrence with stellar mass is as significant as reported by transit observations, the observed occurrence of giants around every fourth FGK star is insufficient to explain the full effect. Furthermore, if the effect acts, it should also produce a strong super-Earth-to-cold-Jupiter anti-correlation, which has not yet been concretely observed in studies (Barbato et al. 2018; Bryan et al. 2019) and could not be proven by this study due to the observational bias highlighted in Sect. 5.4.

Despite the theoretical uncertainties, to test planet formation models, it is interesting to compare synthetic planetary populations against the here-derived occurrence rates. Figure 10 shows the occurrence of synthetic planets based on the modelling by Burn et al. (2021). Their model includes comprehensive disc evolution, growth of planets by planetesimal accretion, N-body interactions, orbital migration, and an evolving planetary envelope and interior resolved in 1D.

The general conclusion from the comparison is that there are too few giant planets forming in the synthetic models around M dwarfs, while small planets are over-produced. The over-production is more significant for the simulation around 0.5 M_⊙ stars compared to the late M dwarfs. However, we note that a similar over-production of model planets is present for solar-mass stars.

Potential model components that should be adjusted are to include pebble accretion, which might promote more growth to giant planets while also only including initial seeds at realistic locations and masses (Voelkel et al. 2021). Furthermore, the radial planetesimal distribution could be less steep or show different features compared to their assumed steep profile (exponent of −1.5). As discussed in detail in Schlecker et al. (2022), the innermost population of planets should be revised around more massive stars.

Fig. 10 Occurrence rates of theoretically calculated planets by Burn et al. (2021). The same bins as in Figs. 3, 4, and 5 are used. Values of sin i are randomly assigned to the synthetic systems following the expected distribution for randomly oriented systems (Schlecker et al. 2022). The top panel shows the synthetic occurrence for the combined discrete stellar masses of 0.3 M⊙ and 0.5 M⊙. The bottom panels show the deviation from observations measured in standard deviations (fp,synthetic – fp)/Sfp for the indicated stellar masses compared to the observational data from Figs. 5 (left) and 4 (right). Where observational upper limits are in agreement with theory, the box is left blank.

7.2.4 1D versus 3D atmospheric models

The definition of the TZ given by Kopparapu et al. (2014) is largely used thanks to the useful and generic equation they provide. However, they used a 1D climate model (Kopparapu et al. 2013) for their calculation. Some studies (e.g. Yang et al. 2013; Kopparapu et al. 2016, 2017; Turbet et al. 2023) showed that intrinsically 3D processes such as the dynamics and the clouds strongly affect the position of the inner edge of the TZ.

It could be interesting to use a more accurate definition of the TZ, based on 3D climate modeling, to revisit the occurrences proposed in this work. In the same way, considering atmospheric compositions consistent with the mass of the planet (e.g. Earthlike for masses around one Earth mass, and H₂-dominated for masses higher than 5 Earth masses) could greatly improve our estimations. This is the aim of a future work.

8 Conclusion

Built upon a prior homogeneous analysis of 197 M dwarfs, our study enables us to present the most precise occurrence rates in RV studies to date. We confirm the prevalence of low-mass, short-period planets, with a rate of 120%${}_{-24}^{+25}$ within a ten-day orbital period for planets under 3 M_⊕, and we will now begin to explore planets with projected masses below 1 M_⊕. The obtained rates are consistent with those from other RV studies, and we propose a new rate for planets at very large separations facilitated by the unique longevity of HARPS. By dividing our sample of M dwarfs by stellar mass, we reveal the expected dichotomy between populations of massive planets at large separations, none of them being detected around the less massive M dwarfs. In addition, extending this study to the insolation calculation yields an occurrence rate of 45.3%${}_{-16}^{+20}$ in the TZ of M dwarfs; this suggests a large number of potential targets for this domain of research given the number of M dwarfs in the solar neighbourhood. Beyond these main results, this study highlight the limitations faced when computing occurrence statistics, even when considering a sample that is analysed homogeneously. Finally, this study also shows the impact of sampling and the choices made over the last 20 years in M-dwarf monitoring programmes.

Acknowledgements

This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France. The original description of the VizieR service was published in A&AS 143, 23 This work has been supported by a grant from LabEx OSUG@2020 (Investissements d’avenir - ANR10LABX56) This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This research has made use of the Exoplanet Follow-up Observation Program (ExoFOP; DOI: 10.26134/ExoFOP5) website, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. NCS acknowledges funding by the European Union (ERC, FIERCE, 101052347). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. This work was supported by FCT -Fundação para a Ciência e a Tecnologia through national funds and by FEDER through COMPETE2020 - Programa Operacional Competitividade e Interna-cionalização by these grants: UIDB/04434/2020; UIDP/04434/2020. N.A.-D. acknowledges the support of FONDECYT project 1240916. X.D acknowledge from support by the French National Research Agency in the framework of the Investissements d’Avenir program (ANR-15-IDEX-02), through the funding of the “Origin of Life” project of the Univ. Grenoble-Alpes. GC acknowledge the financial support of the SNSF (grant number: 200021_197176 and 200020_215760).

Appendix A Tables

Appendix B Graphs

Fig. B.1 Planetary occurrence rates from Sabotta et al. (2021) (orange), those derived from this study (green) and those from Bonfils et al. (2013a) (blue), compared for each range of projected planetary mass and period.

References

All Tables

All Figures

	Fig. 1 Conservative (yellow) and statistical (dark turquoise) detection limits of GJ 163 as function of projected planetary mass and periods, with the respective positions of the detections of GJ 163 b, c, and d. Upper panel shows the limits in a range of periods between 1.6 and 10 000 days, while the lower panel shows a zoomed-in view of these limits for periods between 1.6 and ten days.
In the text

Fig. 2 Iso-frequency map obtained by stacking the 197 statistical detection limits in projected masses as a function of the period. Iso-contours are highlighted for several percentage levels of the sample and describe the distribution of maximum mass missed. The grey part represents the unknown domain, below the lowest statistical detection limit (corresponding to 0.5%), and the white part is the domain beyond 100%.

In the text

Fig. 3 Planetary occurrence statistic. Colours in the background refer to the detection limit map, from the non-reachable region (grey) to the unmissable one (white). The grey points represent the considered detections, while the red points represent the candidates. The grey horizontal and vertical lines delineate the probed domains. The computed occurrence rates (along with their corresponding statistical bounds) are indicated at the centre of each domain.

In the text

	Fig. 4 Planetary occurrence statistics obtained from sub-sample of M-dwarfs with M* > 0.35 M⊙.
In the text

	Fig. 5 Planetary occurrence statistics obtained on sub-sample of M-dwarfs with M* < 0.35 M⊙.
In the text

	Fig. 6 Comparison of occurrence rates obtained as function of minimal number of measurements to consider in the mass-period domains probed.
In the text

	Fig. 7 GJ 163 statistical mass limit as function of insolation. The vertical bars represent the inner (orange: recent Venus; yellow: runaway greenhouse) and outer (green: maximum greenhouse; blue: early Mars) boundaries of the HZ defined in Kopparapu et al. (2013).
In the text

	Fig. 8 Conservative occurrence rate and planets in TZ.
In the text

	Fig. 9 Occurrence of low-mass planets around FGK stars based on Mayor et al. (2011) and obtained using EPOS (Mulders et al. 2018). The marked bins are the ones also used in Figs. 3, 4, and 5, and occurrence is indicated in percent and colour-coded in logarithmic mass against period space.
In the text

Fig. 10 Occurrence rates of theoretically calculated planets by Burn et al. (2021). The same bins as in Figs. 3, 4, and 5 are used. Values of sin i are randomly assigned to the synthetic systems following the expected distribution for randomly oriented systems (Schlecker et al. 2022). The top panel shows the synthetic occurrence for the combined discrete stellar masses of 0.3 M⊙ and 0.5 M⊙. The bottom panels show the deviation from observations measured in standard deviations (fp,synthetic – fp)/Sfp for the indicated stellar masses compared to the observational data from Figs. 5 (left) and 4 (right). Where observational upper limits are in agreement with theory, the box is left blank.

In the text

	Fig. B.1 Planetary occurrence rates from Sabotta et al. (2021) (orange), those derived from this study (green) and those from Bonfils et al. (2013a) (blue), compared for each range of projected planetary mass and period.
In the text