Issue 
A&A
Volume 664, August 2022



Article Number  A65  
Number of page(s)  21  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202142828  
Published online  10 August 2022 
HADES RV Programme with HARPSN at TNG
XV. Planetary occurrence rates around earlyM dwarfs^{★,}^{★★}
^{1}
INAF – Osservatorio Astrofisico di Torino,
Via Osservatorio 20,
10025
Pino Torinese, Italy
email: m.pinamonti.astro@gmail.com
^{2}
INAF – Osservatorio Astronomico di Palermo,
piazza del Parlamento 1,
90134
Palermo, Italy
^{3}
INAF – Osservatorio Astrofisico di Catania,
Via S. Sofia 78,
95123
Catania, Italy
^{4}
Institut de Ciències de l’Espai (ICE, CSIC),
Campus UAB, C/ de Can Magrans s/n,
08193
Bellaterra, Spain
^{5}
Institut d’Estudis Espacials de Catalunya (IEEC),
C/ Gran Capità 24,
08034
Barcelona, Spain
^{6}
Instituto de Astrofísica de Canarias (IAC),
38205
La Laguna, Tenerife, Spain
^{7}
Universidad de La Laguna,
Dpto. Astrofísica,
38206
La Laguna, Tenerife, Spain
^{8}
INAF – Osservatorio Astronomico di Trieste,
via G. B. Tiepolo 11,
34143
Trieste, Italy
^{9}
INAF – Osservatorio Astronomico di Padova,
vicolo dell’Osservatorio 5,
35122
Padova, Italy
^{10}
INAF – Osservatorio Astronomico di Capodimonte,
Salita Moiariello 16,
80131
Napoli, Italy
^{11}
Centro de Astrobiología (CSICINTA),
Carretera de Ajalvir km 4,
28850
Torrejón de Ardoz, Madrid, Spain
^{12}
INAF – Osservatorio Astronomico di Cagliari & REM,
Via della Scienza, 5,
09047
Selargius CA, Italy
^{13}
Dipartimento di Fisica e Astronomia, Università di Padova,
via Marzolo 8,
35131
Padova, Italy
^{14}
Fundación Galileo Galilei  INAF,
Ramble José Ana Fernandez Pérez 7,
38712
Breña Baja, TF, Spain
^{15}
INAF – Osservatorio Astronomico di Brera,
via E. Bianchi 46,
23807
Merate, Italy
Received:
3
December
2021
Accepted:
28
February
2022
Aims. We present the complete Bayesian statistical analysis of the HArpsn red Dwarf Exoplanet Survey (HADES), which monitored the radial velocities of a large sample of M dwarfs with HARPSN at TNG over the last 6 yr.
Methods. The targets were selected in a narrow range of spectral types from M0 to M3, 0.3 M_{⊙} < M_{★} < 0.71 M_{⊙}, in order to study the planetary population around a welldefined class of host stars. We take advantage of Bayesian statistics to derive an accurate estimate of the detectability function of the survey. Our analysis also includes the application of a Gaussian Process approach to take into account stellaractivityinduced radial velocity variations and improve the detection limits around the mostobserved and mostactive targets. The Markov chain Monte Carlo and Gaussian process technique we apply in this analysis has proven very effective in the study of Mdwarf planetary systems, helping the detection of most of the HADES planets.
Results. From the detectability function we can calculate the occurrence rate of smallmass planets around earlyM dwarfs, either taking into account only the 11 already published HADES planets or adding the five new planetary candidates discovered in this analysis, and compare them with the previous estimates of planet occurrence around Mdwarf or solartype stars: considering only the confirmed planets, we find the highest frequency for lowmass planets (1 M_{⊕} < m_{p} sin i < 10 M_{⊕}) with periods 10 d < P < 100 d, , while for shortperiod planets (1 d < P < 10 d) we find a frequency of , significantly lower than for laterM dwarfs; if instead we also take into account the new candidates, we observe the same general behaviours, but with consistently higher frequencies of lowmass planets. We also present new estimates of the occurrence rates of longperiod giant planets and temperate planets inside the habitable zone of earlyM dwarfs: in particular we find that the frequency of habitable planets could be as low as η_{⊕} < 0.23. These results, and their comparison with other surveys focused on different stellar types, confirm the central role that stellar mass plays in the formation and evolution of planetary systems.
Key words: techniques: radial velocities / stars: lowmass / stars: activity / methods: statistical / planets and satellites: detection
All RV and activity data are only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/cat/J/A+A/664/A65
Based on: observations made with the Italian Telescopio Nazionale Galileo (TNG), operated on the island of La Palma by the INAF – Fundación Galileo Galilei at the Roque de Los Muchachos Observatory of the Instituto de Astrofísica de Canarias (IAC); photometric observations made with the APACHE array located at the Astronomical Observatory of the Aosta Valley; photometric observations made with the robotic telescope APT2 (within the EXORAP programme) located at Serra La Nave on Mt. Etna.
© ESO 2022
1 Introduction
Due to the many observational advantages facilitating the detection of rocky planets in close orbits around them, M dwarfs have become increasingly popular as targets for extrasolar planet searches (e.g. Dressing & Charbonneau 2013; Sozzetti et al. 2013; AstudilloDefru et al. 2017). However, these advantages are counterweighted by the difficulties in both detection and characterisation of exoplanetary systems caused by the stellar activity of the host stars, which can produce radial velocity (RV) signals comparable to those from actual planets, leading to possible misinterpretations (e.g. Bonfils et al. 2007; Robertson et al. 2014; AngladaEscudé et al. 2016). Moreover, Mdwarf planetary systems offer an interesting laboratory to test planet formation theories, because they form under different conditions with respect to FGKdwarf systems, with different protoplanetary disk masses, temperature, and density profiles, and gasdissipation timescales (e.g. Ida & Lin 2005).
Several previous studies attempted to compute the occurrence rate of extrasolar planets of different kinds around lowmass stars, taking into account the detection biases of different surveys. Bonfils et al. (2013) analysed the M dwarf sample observed as part of the RV search for southern extrasolar planets with the ESO/HARPS spectrograph (Mayor et al. 2003). Tuomi et al. (2014) studied a similar though smaller sample of M dwarfs observed with both the HARPS and UVES spectrographs (Dekker et al. 2000). The statistical properties of exoplanets in the Kepler M dwarfs sample were analysed by several authors; for example Dressing & Charbonneau (2015) and Gaidos et al. (2016). Their general results pointed towards a high number of lowmass planets orbiting at different distances from their hosts, and confirmed the paucity of giant planets already suggested by earlier surveys (Endl et al. 2003, 2006; Cumming et al. 2008). Recently, Sabotta et al. (2021) computed the occurrence rates from a subsample of 71 stars observed within the CARMENES exoplanet survey. Their sample covered a wide range of hoststar masses between 0.09 M_{⊙} and 0.70 M_{⊙}, and found slightly higher occurrence rates, although compatible with those from Bonfils et al. (2013).
In this framework, the Harpsn red Dwarf Exoplanet Survey (HADES, Affer et al. 2016) programme aims to fully characterise the population of exoplanetary systems on a consistent sample of stars with wellknown properties in order to minimise the effects of varying stellar properties on the planetary frequency measurements. To this end, a catalogue of earlyM dwarfs in a narrow range of stellar masses was selected. HADES is a collaboration between the Italian Global Architecture of Planetary Systems (GAPS, Covino et al. 2013; Desidera et al. 2013; Poretti et al. 2016) Consortium, the Institut de Ciències de l’Espai de Catalunya (ICE), and the Instituto de Astrofísica de Canarias (IAC). Perger et al. (2017a) carried out a performance study of the survey based on the first 4 yr of HADES observations: these authors performed a simulation analysis based on the expected population of earlyM planets, and predicted a yield of significant detections of 2.4 ± 1.5 planets over the whole sample. However, the results of Perger et al. (2017a) were based on the number of observations collected at that time, which more than doubled by the end of the survey (see Sect. 2). Following the same assumptions, in their Fig. 13 these latter authors derived the expected yields of planets for different numbers of targets and observations per target. We can therefore derive the expected number of detected planets corresponding to the final number of HADES observations, that is, 3.8 ± 1.9 planets. To date, 11 planets have been detected as part of the survey (Affer et al. 2016, 2019; Suárez Mascareño et al. 2017; Perger et al. 2017b, 2019; Pinamonti et al. 2018, 2019; ToledoPadrón et al. 2021; GonzálezÁlvarez et al. 2021; Maldonado et al. 2021), which is already exceeding more than 3σ from the adjusted prediction. This shows how the previous knowledge of earlyM dwarf planetary populations was incomplete. Moreover, during recent years, advanced techniques to mitigate stellar noise have been adopted, further improving the detection limits of RV surveys, and increasing the yield of detected planets. The published HADES planets are shown in Fig. 1. Moreover, in Pinamonti et al. (2019) we studied the population of extrasolar planets orbiting M dwarfs detected using the RV technique. We found moderate evidence of a correlation between planetary mass and stellar metallicity, with different behaviours for small and large planets, which is expected from theory (e.g. Mordasini et al. 2012). In Maldonado et al. (2020), we analysed a large sample of M dwarfs in a homogeneous way, and confirmed that giant planet frequency exhibits a correlation with metallicity, while the same does not seem to be true for lowmass planets; these trends are similar to those previously observed for FGK stars (e.g. Sousa et al. 2008; Adibekyan et al. 2012). Pinamonti et al. (2019) also confirmed the effect pointed out by Luque et al. (2018) of different mass distributions between single and multiple planetary systems. However, these studies were conducted over a heterogeneous sample of planets detected from different surveys, which could not be corrected for detection biases. Therefore, their results are to be considered preliminary, because observational biases may have a strong influence on the observed distributions.
In this work, we present a thorough Bayesian characterisation of the HADES planetary population statistics and global detectability of the survey. We model our analysis on the Bayesian approach proposed by Tuomi et al. (2014). We implement this technique with the Monte Carlo Markov chain (MCMC) and Gaussian Process (GP) regression framework, which we successfully applied in the detection of several of the planetary systems of the survey (e.g. Affer et al. 2016; Pinamonti et al. 2018, 2019).
In Sect. 2, we describe the HADES programme and summarise the characteristics of its targets. In Sect. 3, we detail the MCMC analysis technique and main results in the identification of planetary signals, while in Sect. 4 we discuss the detection limits of the survey derived from our Bayesian analysis. In Sect. 5, we discuss the planetary occurrence rates we derive for early M dwarfs, and in Sect. 6, we compare them with previous results of similar surveys and also with the statistics of FGKdwarf systems. We summarise our conclusions in Sect. 7.
Fig. 1 Overview of the HADES detected planetary systems. The published planets of the sample are shown as red circles: the symbol size is proportional to the minimum planetary mass. The conservative and optimistic limits of the habitable zone of each system (see Sect. 5.1) are shown as thick dark green and light green bands, respectively. 
2 The survey
As discussed in Affer et al. (2016), the initial HADES sample was constructed by selecting 106 targets from the Lépine & Gaidos (2013) and PMSU (Palomar/Michigan State University, Reid et al. 1995) catalogues: the stars were required to have spectral types between M0 and M3, with masses between 0.3 M_{⊙} and 0.7 M_{⊙}, a visual magnitude V < 12, and both to have a high number of Gaia mission (Gaia Collaboration 2016) visits and to be part of the target catalogue of the photometric survey APACHE (A PAthway towards the CHaracterization of Habitable Earths, Sozzetti et al. 2013). Several targets were successively discarded, as they were discovered to be illsuited to planet searches because of close binary companions, fast rotation (v sin i > 4 km s^{−1}), very high chromospheric activity levels (log R’_{HK} > −4.3), incorrect spectral type (T_{eff} > 4000 K), or for having been classified as subgiants/giants (log ɡ < 4). All the stellar parameters of the targets were consistently rederived by Maldonado et al. (2017) from the HADES HARPSN measurements by applying the techniques from Maldonado et al. (2015), which allowed precise estimates of spectral types, metallicities, and effective temperatures. For the current study, we also selected only targets that were observed at least ten times, because with fewer observations our Bayesian analysis technique discussed in Sect. 3 could not be robustly applied.
Applying the aforementioned selection criteria, the HADES sample was reduced to 56 targets, with spectral types between M0 and M3. The parameters important for this analysis of all the stars in the sample, as computed by Maldonado et al. (2020), are listed in Table 1. The RVs were derived with the TemplateEnhanced Radial velocity Reanalysis Application (TERRA, AngladaEscudé & Butler 2012): this templatematching approach has proven to be more effective than the CCF technique of the standard HARPSN Data Reduction System pipeline (DRS, Lovis & Pepe 2007), in particular in the analysis of earlyM dwarf spectra (Affer et al. 2016; Perger et al. 2017a). Moreover, we computed the TERRA RVs considering only spectral orders redder than the 22nd (λ > 453 nm)^{1} in order to avoid the lowsignaltonoiseratio (S/N) blue part of the M dwarfs spectra. The observations were carried out from August 2012 to December 2018, with an average timespan 〈T_{s}〉 = 1760 d (median 1890 d). The average number of RVs per target is 77 (median 60), which is more than twice the average number of observations collected at the time of the performance study by Perger et al. (2017a). Moreover, it is worth noting that this typical number of observations per target is significantly larger than the average number of observations used by Bonfils et al. (2013) (20), and Tuomi et al. (2014) (62). All RV and activity data used in this analysis are available at the CDS in machinereadable format.
Stellar rotation periods
Stellar surface inhomogeneities often produce RV signals close to the rotation period of the star, which can easily mimic the planetary signals (e.g. Vanderburg et al. 2016). This is particularly significant around M dwarfs, because these stellar signals have periods and amplitudes comparable to those of rocky planets in the habitable zone of their stars (Newton et al. 2016). For these reasons, several studies attempted to obtain a precise characterisation of the stellar rotation periods and activity features of the HADES sample: Maldonado et al. (2017) and Scandariato et al. (2017) studied the activity indicator behaviour in the M dwarf sample monitored by the HADES programme, identifying many useful relationships between activity, rotation, and stellar emission lines; GonzálezÁlvarez et al. (2019) studied the Xray luminosity as a proxy of the coronal activity of the stars in the sample, and its relationship with the stellar rotation periods. Suárez Mascareño et al. (2018) instead used spectral activity indicators and photometry to investigate the presence of signatures of rotation and magnetic cycles in 71 HADES M dwarfs, providing activityindex derived rotation periods for 33 of those stars and 36 rotation periods estimated from activity—rotation relationships. All the 33 targets with activityderived rotation periods are included in the sample of the current analysis, along with 23 of those with periods derived from activity relationships. The relevant rotation periods are listed in Table 1. Moreover, 23 HADES targets have photometric rotation periods detected in the APACHE light curves (Giacobbe et al. 2020), which shows a good correspondence with the spectroscopically derived rotation periods, being either equal or with one close to the first harmonic of the other.
This information over the stellar rotation periods of the sample is very important for our study, because it is paramount to the identification of real Keplerian and activityinduced RV variation during the Bayesian modelling of the time series, as discussed in Sect. 3.4. It is worth noting that, in the analysis of the impact of activity, the active region lifetime is also an important parameter (e.g. Giles et al. 2017; Damasso et al. 2019); however, as not all of our stars have timeseries that are sufficiently well sampled for the application of sophisticated analysis techniques such as GP regression, we focus on the rotation period as the main parameter for stellar activity modelling.
3 MCMC analysis
3.1 Statistical analysis technique
Our analysis technique follows the approach by Tuomi et al. (2014), who used Bayesian statistics to estimate the occurrence rate and detectability function of extrasolar planets around M dwarfs.
For each timeseries of the sample, the observed number of planets in a given periodmass interval Δ_{P,M} can be expressed as: (1)
where the subscript i indicates the ith system: f_{occi} is the occurrence rate of planets around the ith star, and p_{i} is the detectability function, that is, the probability of detecting a planet in the parameter interval Δ_{P,M} in the ith time series. f_{obs,i} is simply computed as the number of planets detected from each timeseries in each parameter interval.
Assuming that k_{i} Keplerian signals were detected in the ith dataset, the Bayesian technique from Tuomi et al. (2014) estimates the detectability function p_{i} by applying a k_{i} + 1signals model to the timeseries by means of a sampling algorithm. This way, the region of the (P, M) parameter space explored by the parameters of the hypothetical k + 1th testplanet, which cannot be significantly detected in the data, represents the region of the parameter space where the detection technique is not able to reveal significant signals, that is, the RV amplitudes are so low that the likelihood function does not rule out the presence of a hypothetical additional planetary signal. More precisely, the areas that were not explored by the Markov chains are the areas where signal detection could have been possible, if significant signals were actually present in the data. The detectability function can therefore be approximated as: (2)
where is equal to one in the regions of parameters space explored by the k + 1th Keplerian, and zero otherwise (as detailed in Sect. 3.2.2). An example of the detection function p_{i} derived with this method is shown in Fig. 2 for the planethosting star GJ 15 A (Pinamonti et al. 2018).
The observed frequency of planets in a given range of orbital period and planetary mass over the whole sample of N datasets can be computed by the sum over i of Eq. (1): (3)
assuming the occurrence rate f_{occ} to be common for all stars in the sample, f_{occ} = f_{occ,i} for all i. This is expected to be the case for a welldefined sample of host stars such as the HADES sample described in Sect. 2.
As f_{obs}(Δ_{P,M}) is known and the square brackets term can be easily calculated, Eq. (3) allows us to calculate the occurrence rate of planets across the (P, M) parameter space. The global detection probability function of the sample, p can instead be computed simply by dividing the square brackets term by N: (4)
Fig. 2 Detection function map of the RV timeseries of GJ 15 A. The white part corresponds to the area in the period—minimum mass space where additional signals could be detected if present in the data, i.e. p_{i} = 1, while the black region corresponds to the area where the detection probability is negligible, i.e. p_{i} = 0. The red circle marks the position in the parameter space of the planet GJ 15 A b. 
Stellar parameters and summary of the observations of the 56 HADES targets analysed in this work.
3.2 Algorithm implementation
This Bayesian technique was applied using an adapted Python version of the publicly available emcee Affine Invariant MCMC Ensemble sampler by ForemanMackey et al. (2013) as a sampling algorithm. We also take advantage of GP regression to treat stellaractivity correlated noise from RV timeseries, which was applied through the GEORGE Python library (Ambikasaran et al. 2015). These are the same algorithms that were previously and successfully applied in the study of several HADES targets (e.g. Affer et al. 2016; Pinamonti et al. 2018, 2019). We used a variable number of random walkers, N_{walk}, to sample the parameter space, depending on the number of parameters, N_{par}, of the final model applied to each timeseries, always selecting N_{walk} ≳ 10 N_{par} to ensure a good sampling of the parameter space. The convergence of the different MCMC analyses was evaluated calculating the integrated correlation time for each of the parameters, stopping the code after a number of steps equal to 100 times the largest autocorrelation times of all the parameters (ForemanMackey et al. 2013), after applying an initial burnin phase (Eastman et al. 2013).
The general form for the RV models that we applied to all RV timeseries in this work is given by the equation: (5)
where N_{pl} is the number of planets in the system, and the RV Baseline Model, ΔRV_{BM}(t), is defined as: (6)
where γ is the systemic velocity, d is the linear acceleration, and is the mean epoch of the timeseries. The Keplerian RV model ΔRV_{Kep,j}(t) is defined as: (7)
The loglikelihood function, which is maximised by our fitting algorithm, ln ℒ, is defined as follows: (8)
where y(t) is the ‘clean’ RV timeseries (see Sect. 3.3), σ(t) is the RV total uncertainty at epoch t, computed as , with σ_{data} the RV HARPSN internal error, and σ_{jit} the additional uncorrelated ‘jitter’ term, which is fitted by the model to take into account additional sources of error such as uncorrelated stellar noise or instrument drifts. For a brief discussion on all the fitted parameters and the adopted priors, see Appendix A.
3.2.1 Model selection
We identified the best model for each timeseries by computing the Bayesian information criterion (BIC, Schwarz et al. 1978), defined as: (9)
and, starting from the RV Baseline Model, ΔRV_{BM}(t), accepting a model with a higher number of parameters when its BIC value was at least 10 points lower than the previous model, which is strong statistical evidence in favour of the lowerBIC model (Kass & Raftery 1995)^{2}.
3.2.2 Test planet model
When the best model, namely that containing N_{pl} planetary signals, has been identified for each timeseries, in order to compute the detectability function p_{i}, we run one additional MCMC fit adding one additional ‘Keplerian’ term to Eq. (5), which corresponds to the testplanet used to sample . The RV model therefore becomes: (10)
where is the testplanet RV model. It is worth noting that, as the N_{pl}planet model was already the ‘bestmodel’ for the corresponding timeseries, the N_{pl} + 1 test model will not be statistically favoured with respect to the former, i.e. ΔBIC < 10. Moreover, for the sake of computational efficiency, we restricted the testplanet RV model to the circular case, imposing e_{test} = ω_{test} = 0, i.e.:
This assumption should not strongly affect our results, because eccentricity up to 0.5 is not expected to have a strong impact on detection limits (Endl et al. 2002; Cumming & Dragomir 2010). Finally, from the posterior distribution of P_{test}, and that of K_{test} converted to minimummass, it is possible to compute as discussed above.
3.3 Stellar noise correction
As previously discussed, an important issue in the analysis ofMdwarf RV timeseries is the presence of stellar noise. It is known that not all stellar RV variations are directly related to the stellar rotation period (e.g. Dumusque et al. 2015, and references therein), and therefore we had to adopt a multistep approach to correct as much as possible of the correlated and uncorrelated stellar noise in order to improve the estimated planetary parameters, and the derived detectability function.
For all the HADES targets selected for this study, we computed the stellar activity indexes based on the Ca II H and K, Hα, Na I D_{1} D_{2}, and He I D_{3} spectral lines, applying the procedure described in Gomes da Silva et al. (2011) to all the available HARPSN spectra. We then applied the following two correction criteria to all timeseries.
Outlier removal. first we used the activity index timeseries to identify any potential outlier caused by transient stellar effects such as flares. This can be very important, because flares have been observed to be very common even around earlytomid M dwarfs (Hawley et al. 2014), and have been observed to manifest even around stars showing little to no periodic stellar signals produced by active regions (Yang et al. 2017). To identify these outliers around each target, we performed a 3σ clipping around the median of each activity timeseries; if any longterm trend was present in the activity timeseries, the clipping was performed over the linearly detrended timeseries in order to prevent the longterm variation from enlarging the standard deviation of the data, which could hide significant outliers. We then removed the epochs corresponding to the identified outliers from all the activity and RV timeseries for that target. The number of outliers identified in each timeseries is listed in Table 1.
Correlation with activity indexes. Another known stellar phenomenon that needs to be accounted for is that of magnetic cycles, which have been observed to produce significant variations in RV timeseries (e.g. Lovis et al. 2011; Dumusque et al. 2011). This effect can be identified by the presence of a longterm correlation between the RV and spectroscopic activity index timeseries. For this reason, after removal of outliers, for each target we computed the Pearson correlation coefficients between the RVs and all activity indicators, and if any strong correlation was found (ρ > 0.5) we detrended the RV timeseries via a linear fit with the corresponding index timeseries^{3}. We chose to take into account only strong correlations, because even if lower correlations might still suggest a magnetic cycle contamination in the RVs, they would be more difficult to correct via a simple linear fit, with the risk of injecting additional noise into the timeseries if the RVs were inappropriately detrended. The timeseries were corrected for significant activity correlation, and the correlated activity indexes are listed in Table C.1.
These corrections were applied to all the analysed systems, producing the ‘clean’ RV timeseries to which the MCMC algorithm was applied to estimate the planetary parameters and compute the detection limits as previously discussed. In addition to these, we adopted additional steps for those systems with correlated periodic stellar signals in the RV timeseries corresponding to either the star’s rotation period or its first harmonics. We used two alternative techniques to model these stellar signals. The first approach was to fit the activity RV signals as a sine wave, that is, to include them in the general model of Eq. (5) (in which case the summation is computed over N_{signals} = N_{pl} + N_{act}, including the number of activity signals, N_{act}).
The second approach was to apply GP regression to those systems that presented strong activity RV signals. We adopted the GP regression only for those systems that had a sufficiently large number of observations in order to avoid overfitting the data due to the adaptive nature of the GP. We chose this threshold to N_{obs} > 70 RVs per target in order to have at least ten times the number of minimum parameters in the model (3BM + 4GP). For those systems, we computed the likelihood function to be maximised by the MCMC sampler via GP regression, that is, the loglikelihood function in Eq. (8) was redefined as: (11)
where r is the residual vector r = y(t) − ΔRV(t), with ΔRV(t) either from Eq. (5) or Eq. (10), and K is the covariance matrix. The covariance matrix is defined by the GP kernel: in this work we adopted the commonly used quasiperiodic (QP) kernel, which is defined as the multiplication of an expsinsquared kernel with a squaredexponential kernel (Ambikasaran et al. 2015). The QP kernel can be expressed as (12)
where K_{i,j} is the ij element of the covariance matrix, and the covariance is described by four hyperparameters: h is the amplitude of the correlation, λ is the timescale of decay of the exponential correlation, θ is the period of the periodic component, and w is the weight of the periodic component. The last term of Eq. (12) describes the white noise component of the covariance matrix, where δ_{i,j} is the Kronecker delta and σ is the RV total uncertainty defined as in Eq. (8).
Finally, it is important to point out that both these approaches to the fitting of correlated stellar signals were applied only when the model including the stellaractivity correction was preferred over the previous model by our model selection criterion discussed in Sect. 3.2.1. The adopted priors for the GP hyperparameters are discussed in Appendix A, while in Table C.1 all the timeseries on which GP regression was applied are listed, as well as all timeseries in which sinusoidal activity signals were fitted.
3.4 Planetary candidates and activity signals
We now briefly discuss the Keplerian RV signals and stellar activity features identified in the bestfit models of our analysis. A summary of all the adopted models for each system in the survey is given in the Appendix in Table C.1.
Published systems. We used the algorithm discussed in the previous sections to reanalyse all the published HADES systems (see Fig. 1) in order to compute the detection function p_{i}. These independent analyses confirmed the previous results, producing a good agreement in the fitted planetary and activity parameters, usually consistent within 1σ with the published bestfit values. Some systems required additional care in the reanalysis, because in the present study we used only the HADES HARPSN RV timeseries collected during the programme, while other studies took advantage of additional RV data from the literature or from other instruments (e.g. Pinamonti et al. 2018; Perger et al. 2019). In the case of GJ 15 A in particular, it was not possible to model the longperiod planet GJ 15 Ac with only the HARPSN data, which have a much shorter timespan than its orbital period (Pinamonti et al. 2018): we therefore decided not include the longperiod Keplerian of planet c in the RV model, but to use the acceleration term d of the Baseline Model to fit its contribution in the HADES dataset^{4}.
Planetary candidates. In addition to the already published HADES planetary systems, there are a few other systems in which our analyses identified promising periodic signals, with no obvious relationship with the stellar rotation period or other activitydriven phenomena, and therefore could be considered planetary candidates. However, to ascertain their true planetary nature would require an indepth analysis which is beyond the scope of this paper. For this reason we do not present detailed analyses or complete orbital parameters for these systems, which will be discussed in future specific publications, two of which are currently in preparation: González Hernández et al. and Affer et al. will discuss GJ 21 and GJ 3822 systems, respectively. However, as these additional planets would affect the occurrence rates derived from our survey, we take into account the presence of these additional candidates, as discussed in Sect. 5. The bestfit orbital periods and minimum masses of these candidates are listed in Table 2, and in Appendix B a brief overview of the corresponding RV signals is presented.
Activity signals. Many additional significant signals have been identified in the modelling of the remaining HADES timeseries, which either had a direct counterpart in the activity indexes or were classified as being of stellar origin after further investigation. All these activity signals were still included in our models, either as quasiperiodic signals in the GP regression (Eq. (12)) or as sinusoidal signals included in the general model in Eq. (5). A general discussion of the RV signals of stellar origins identified in the HADES timeseries can be found in Suárez Mascareño et al. (2018).
Candidate planetary signals detected in the RV modelling.
3.5 Longterm trends
In addition to the periodic RV signal identified during the MCMC analysis, another interesting aspect is the presence of longterm trends in the RV timeseries: these could indicate the presence of longperiod companions, which could not otherwise be detected due to the limited temporal baseline of the survey. To account for this, our RV Baseline Model applied to all systems, ΔRV_{BM}(t), always includes an acceleration term, (see Eq. (6)).
In Table 3 the targets that presented significant (> 3σ) accelerations in the final RV model are listed. It is worth noting that, as discussed in Sect. 3.3, the analysed RV timeseries were checked and corrected for correlations with the activity indexes in order to reduce longterm signals due to stellar activity or magnetic cycles. However, two of the targets that showed significant longterm trends, GJ 119A and GJ 694.2, were previously identified as having longterm activity cycles, expected to produce the RV variation (Suárez Mascareño et al. 2018), even though we did not find significant correlation between RVs and activity indexes. Nevertheless, the RV trend in the timeseries for GJ 694.2 has a peaktopeak amplitude of >100 m s^{−1}, which is extremely large compared to other M dwarfs with detected RV activity cycles, and is therefore probably produced by other sources (see below).
A few systems showed additional longterm evolution, with significant curvature that could not be fitted via a simple linear trend, and therefore we introduced in their analyses a quadratic acceleration term, changing the Baseline Model in Eq. (6) to: (13)
where q is the quadratic acceleration coefficient. These quadratic trends were included in the final RV models whenever their addition to the model was accepted by our model selection criteria (see Sect. 3.2.1). The significant quadratic longterm trends detected in our sample are listed in Table 4. Of these three targets, GJ 119A shows a significant quadratic acceleration, but this is suspected to be related to activity cycle previously identified by Suárez Μascareño et al. (2018). Moreover, GJ 649.2 and GJ 3997 both have massive companions which are most probably the source of the observed trends (see below).
We checked the targets showing longterm trends to identify any stellar companions that could be the origin of the observed RV shift. GJ 793 and GJ 4092 are both single stars with no known stellar companion. GJ119A, GJ412 A, and StKM1650 all have common propermotion companions, but with projected separations of 340 AU, 160 AU, and 980 AU, respectively, the maximum RV acceleration produced by the stellar companions^{5} is on average one order of magnitude lower than the bestfit d values listed in Table 3. Therefore, the longterm trends observed in these systems cannot be ascribed to the nearby stellar companions, but could in fact be due to the presence of asofyet unknown perturbers. Moreover, GJ 3997 and NLTT 21156 both have another stellar object with a small onsky separation, but the lineofsight separation measured from the Gaia EDR3 parallaxes is in both cases >10^{4} AU, which again means that the nearby objects could not produce the observed RV trends.
Finally, we checked for propermotion anomalies (PMAs) between Gaia Early Data Release 3 (Gaia EDR3) and HIPPARCOS (Kervella et al. 2022), and other indications of unresolved massive companions in Gaia EDR3, such as Renormalised Unit Weight Error (RUWE) and excess noise (Gaia Collaboration 2021). Three targets show significant PMA, GJ 15A, GJ 694.2, and GJ 119 A: in GJ 15A the anomaly is probably due to the stellar companion GJ 15B; in GJ694.2 the PMA is caused by a subarcsecond stellar companion, which is also the probable source of the measured RV longterm variation; in GJ 119A the PMA is also probably caused by the stellar companion, which is nevertheless too distant to produce the observed RV acceleration (as discussed above). Two of the other targets, StKM1650 and GJ 3997, while not having PMA measurements, present high RUWE and excess noise values in Gaia EDR3, which are potential evidence of orbital motion, indicating the presence of possible brown dwarf companions (Lindegren et al. 2021). It is worth noticing that these analyses on the presence of longperiod planetary companions in the HADES sample are still preliminary, but a more indepth study to ascertain the nature of the observed longterm RV trends goes beyond the scope of this paper.
Significant linear trends identified in the RV modelling.
Significant quadratic trends identified in the RV modelling.
4 Detection limits
To derive the detectability function, we first mapped the parameter space into a 150 × 150 logarithmic grid: the period range covered was [1,3000] d, obtained from the prior of P_{test} defined in Appendix A; the investigated minimummass range was instead [1,1000] M_{⊕}. We then used this (P, M) grid and the posterior distributions derived from the MCMC analysis of each HADES system to compute , from which the global detection probability function p of the HADES sample could be computed as the average of the detection functions of all the systems, as shown in Eq. (4). The resulting detection map of the survey is shown in Fig. 3.
As expected, the detectability function increases for larger masses and shorter periods. For periods between 1 d and 10 d, the average p = 90% detection level corresponds to m_{p} sin i = 9.3 M_{⊕}, while for the longest considered periods, [1000,3000] d, it corresponds to masses as large as m_{p} sin i ≃ 180 M_{⊕}.
It is worth noting that most of the planets and candidates are found in the region of the parameter space with intermediate detectability, around the p = 50% level, while a few are found in regions with very low p (around the targets with the largest and bestsampled timeseries).
5 Occurrence rates
Given the detectability function, p, the planetary occurrence rate, f_{occ}, can be computed as the average number of planets per star in a given region of the parameter space Δ(P, M) from the Poisson distribution: (14)
where k = k_{Δ(P,M)} is the number of planets detected in the chosen parameter space interval Δ(P, M), and the expected value is computed as the product between n = n_{Δ(P,M)}, the number of targets sensitive to planets in Δ(P, M), and f_{occ}. The number of sensitive targets, n_{Δ(P,M)}, can be computed as the mean value of p over Δ(P, M) multiplied by the number of stars in the survey N^{6}. We can compute the planetary occurrence rate from Eq. (14) by considering it as the (unnormalised) posterior distribution of f_{occ}, for given values of k and N. Therefore, bestfit values of planetary occurrence can be computed as the medians and 68% confidence intervals. In the case of k = 0, that is no detected planet in a given interval, the upper limit of the planetary occurrence rate was computed as the 68th percentile of the distribution.
To compute the occurrence rates f_{occ}, we defined the following parameter space interval Δ(P, M): the m_{p}, sin i was divided in three intervals, [1., 10.] M_{⊕}, [10., 100.] M_{⊕}, and [100., 1000] M_{⊕}; the P axis was divided in four intervals, [1,10] d, [10,10^{2}] d, [10^{2},10^{3}] d, and [10^{3},3 × 10^{3}] d. We chose these large intervals, even if most of our detections are grouped at short periods, because these are the most common period and mass intervals used in the literature for computing the planetary occurrence rates. In particular, the two previous studies on Mdwarf planetary populations from RVs of Bonfils et al. (2013) and Tuomi et al. (2014) use these same definitions, making a comparison between our results and theirs possible. The averaged detection function 〈p〉_{Δ(P,M)} is shown in Fig. 4.
We computed the occurrence rates first using only the confirmed HADES planets, shown in Fig. 1, inserting in Eq. (14) the number of confirmed planet per parameter space interval k_{p} = k_{p,Δ(P,M)}. We then also considered the candidate planets presented in Sect. 3.4, thus computing the total number of both confirmed planets and candidates per parameter space interval k_{c} = k_{c,Δ(P,M)}. The resulting occurrence rates derived from the confirmedonly and confirmed+candidates planets are listed in Tables 5 and 6, respectively.
In addition to the 2D occurrence rates shown in Tables 5 and 6, we computed the 1D occurrence rate distribution for both period and minimum mass. This was done using the same bins defined previously for the 2D map, but integrating over the whole range of the other parameter, thus computing the frequency of planets of all masses as a function of period, and of planets of all periods as a function of minimum mass. The results are shown in Figs. 5a and 6a. Furthermore, we computed the cumulative planet frequency for both the period and minimum mass. This was done to show the finer structure of the planetary occurrence rates. In both cases, we restricted the computation over the intervals [1,100] d and [1,100] M_{⊕} in order to focus on the region of the parameter space in which there were detected planets and candidates. We then divided these intervals into a fine grid of bins in order to reduce the number of planets in each bin by as much as possible. We used 80 bins for the period distribution, and 100 bins for the minimum mass distribution. The cumulative planetary frequency was computed as the occurrence rate, integrating the detection function from the minimum value to the value of each bin, marginalising over the other parameter. The results are shown in Figs. 5b and 6b.
Fig. 3 Global HADES detection map. The color scale expresses the global detection function, p. The red circles mark the position in the parameter space of the confirmed HADES planets (see Fig. 1), while the yellow circles indicate the additional candidates presented in this study (see Table 2). 
Fig. 4 Averaged HADES detection map. The color scale expresses the global detection function, p, averaged over the intervals Δ(P, M) defined for the computation of the occurrence rates. The red circles mark the positions in the parameter space of the confirmed HADES planets (see Fig. 1), while the yellow circles indicate the additional candidates presented in this study (see Table 2). 
5.1 HZ occurrence rates
An aspect of great interest in the studies of occurrence rates of lowmass planets is the parameter η_{⊕}, the frequency of lowmass habitable planets, that is, planets with masses m_{p} sin i < 10 M_{⊕} orbiting at a distance from their host star that allows the presence of liquid water on their surface, in other words, within the socalled habitable zone (HZ, Kasting et al. 1993). We compute the HZ limits following the recipe from Kopparapu et al. (2013b), defining a conservative and an optimistic HZ: the conservative HZ is computed adopting the runaway greenhouse and maximum greenhouse coefficients for its inner and outer limits respectively; the optimistic limits of the HZ, a_{HZ,in} and a_{HZ,out} are instead computed with the recent Venus and early Mars coefficients, respectively (Kopparapu et al. 2013a).
Thus far no planet has been confirmed inside the conservative or optimistic HZs of any HADES target (see Fig. 1). However, we can still compute an upper limit of the frequency of habitable planets. Therefore, we generated an additional detection map of the HADES samples as a function of the minimum mass and of the position inside the HZ: this was done for each target by computing the detectability functions p_{i} as in Sect. 4, but the posterior of the period P_{test} was first converted into the posterior of semimajor axis a_{test} and then the semimajor axis was converted in HZ logarithmic scale, a_{test,HZ}: (15)
this way a_{test,HZ} = 0 corresponds to the inner edge of the optimistic HZ, a_{test,HZ}, and a_{test,HZ} = 1 to the outer edge, a_{test,out} This conversion allow us to compare the HZs of different stars, which are located at different intervals of the semimajor axis parameter space.
For each star, we computed over a 100 × 100 logarithmic grid in the [a_{HZ}, M] parameter space, where the range of the converted semimajor axis, a_{HZ}, was [0,1], and the minimummass range was [1,100] M_{⊕}. We then computed the global detection function inside the HZ p_{HZ} following Eq. (4). The result is shown in Fig. 7.
From the detection function p_{HZ}, taking into account only the confirmed planets none of which were detected inside the HZ (i.e. k_{HZ} = 0), we were able to use Eq. (14) to compute the upper limits of the occurrence rate f_{occ,HZ} = η_{⊕} as before: for lowmass planets, m_{p} sin i < 10 M_{⊕}, we obtained an upper limit f_{occ,HZ} = η_{⊕} < 0.23. Moreover, we also computed the occurrence rate of highermass planets inside the HZ, that is 10 M_{⊕} < m_{p} sin i < 100 M_{⊕}, obtaining f_{occ,HZ} < 0.03.
Moreover, we detected one potential planetary signal located inside the optimistic HZ of its host star: GJ 399 32.9 d candidate is located inside the HZ close to its inner edge. However, this signal is slightly less significant than the threshold we adopted in Sect. 3.2.1 (ΔBIC = 9.5), and is therefore not considered as a planetary candidate in Sect. 3.4 and Table 2. Nonetheless, we chose to take it into account to compute a second estimate of η_{⊕}. Taking this candidate into account, the occurrence rate for m_{p} sin i < 10 M_{⊕} increases to .
Occurrence rates of planets in the HADES sample.
5.2 Longperiod companions occurrence rates
It is more difficult to estimate the occurrence rates of longperiod planets, that is those with P ≳ 3000 d, which we were only able to detect as linear RV trends due to the relatively short timespan of our survey. As we detected some linear trends with no apparent stellar origin, as discussed in Sect. 3.5, we discuss a few approximations and caveats adopted in order to try and estimate the occurrence rates of the associated longperiod planets. It should be noted that the uncertainties on such estimates are intrinsically very large, and very strong approximations are required, but the results could provide interesting information about the littleknown field of longperiod planets around M dwarfs, and therefore we proceeded anyway.
First we needed to estimate the planetary parameters from the information we can obtain from the linear trends: were able to compute an estimate of the minimum orbital period P_{min} in the approximation of a circular orbit as four times the timespan of the observations P_{min} ≳ 4 T_{s}^{7}; similarly, we computed the minimum RV semiamplitude K_{min} as the total RV variation during the observations, that is K_{min} ≳ T_{s} d. From P_{min} and K_{min}, the corresponding ‘minimum’ minimum mass (m_{p} sin i)_{min} can be derived. The resulting minimum planetary parameters are listed in Table 7. In the case of the GJ 3997 quadratic trend, we were able to approximate the minimum RV amplitude K_{min} as the difference between the maximum and minimum RV values of the quadratic curve during the timespan of the observations, K_{min} = RV_{max} − RV_{min}, while we computed the period as a function of the quadratic acceleration, q, as (Kipping et al. 2011).
Another obstacle to overcome is the computation of the detection efficiency for such longperiod signals, which cannot be computed following our standard recipe described in Sect. 3.1 because the emcee testplanet approach cannot be robustly applied to periods much longer than the timespan of the timeseries. We decided to assess the detection efficiency of the HADES timeseries towards longperiod signals, that is P ≳ 3000 d, by computing for each detected longterm trend d the number of timeseries in which that trend could be detected with a significance larger than >3σ, that is in which d > 3 σ_{d,i}, where σ_{d,i} is the measured uncertainty of the acceleration term in the bestfit MCMC solution of the ith timeseries. The resulting numbers of sensitive timeseries for each trend, n_{t}, are listed in the last column of Table 7. The mean value of n_{t} can be then used along with the number of significant trends, k_{t}, to estimate the occurrence rate f_{occ,t} from Eq. (14).
As discussed in Sect. 3.5, we conducted a few tests to ascertain whether or not the measured RV trends could be caused by other phenomena such as magnetic cycles or stellar companions. A more indepth analysis is beyond the scope of the present work, but we can compute some preliminary estimates of the occurrence rates, assuming that all the observed trends are due to actual longperiod planetary companions. The resulting frequency of planets with P ≳ 3000 d and m_{p} sin i ≳ 10 M_{⊕} is therefore As the minimum masses listed in Table 7 are quite diverse, and cover two different intervals of mass defined in the analysis of the shortperiod planet occurrence rates (Fig. 4), we were able to divide the trend sample into two subsamples with minimum masses larger and smaller than 100 M_{⊕}. Computing the occurrence rates for these two subsamples separately, we obtain for longperiod planets with m_{p} sin i > 100 M_{⊕}, and for longperiod planets with 10 M_{⊕} < m_{p} sin i < 100 M_{⊕}.
Fig. 5 Planetary occurrence rate f_{occ} as a function of the orbital period, with the corresponding 1σ uncertainties (top panel). The median cumulative planet frequency is shown in the lower panel. The red and yellow distributions correspond to the occurrence rates derived from the confirmedonly, k_{p}, and confirmed+candidates planets, k_{c}, respectively. 
Fig. 6 Planetary occurrence rate f_{occ} as a function of the minimum mass, with the corresponding 1 σ uncertainties (top panel). The median cumulative planet frequency is shown in the lower panel. The red and yellow distributions correspond to the occurrence rates derived from the confirmedonly, k_{p}, and confirmed + candidates planets, k_{c}, respectively. 
Fig. 7 Habitable zone HADES detection map, derived as described in Sect. 5.1. The colour scale expresses the detection function, as in Fig. 3. 
Longperiod candidate planetary signals assumed from the modelled linear RV trends.
6 Discussion
We now discuss the main features of our results concerning the sensitivity of the HADES survey and the occurrence rates of extrasolar planets around early M dwarfs, comparing them with previous studies on the planetary population of M dwarfs and earliertype stars, both from RV surveys and photometric surveys such as Kepler. It is worth noting that the comparison between RV and photometric surveys is far from trivial, because it strongly depends on the mass—radius relations which are still largely unknown, in particular for small lowmass planets.
6.1 Comparison with other RV surveys
We can see in Fig. 6a that the occurrence rate of planets around earlyM dwarfs increases strongly towards lower masses. If we take into account only the confirmed HADES planets (red distribution), the occurrence rate rises from , for 10 M_{⊕} < m_{p} sin i < 100 M_{⊕}, to , for 1 M_{⊕} < m_{p} sin i < l0 M_{⊕}. This confirms the expected behaviour that was previously observed in other surveys of laterM dwarfs (Bonfils et al. 2013; Tuomi et al. 2014). The occurrence rate as a function of the orbital period, as shown in Fig. 5a, instead appears to peak at intermediate periods, with the highest value , for 10 d< P < 100 d, which is roughly 2σ higher than the value for P < 10 d,. These behaviours are accentuated if we also take into account the additional planetary candidates presented in Sect. 3.4 (yellow distributions in Fig. 5). It is worth noting that this is due to the fact that most of the detected planets and candidates in the HADES samples are gathered around periods of a few tens of days and masses of around 10 M_{⊕} (Fig. 3), even if the detection probability in that region of the parameter space is quite low, p ≃ 50%. This confirms what was previously reported by Tuomi et al. (2014) in their study of UVES and HARPS RV data, and could also correspond to the overabundance of transiting planets with similar periods and radii smaller than 3 R_{⊕} observed around Kepler M dwarfs: Dressing & Charbonneau (2015) observed a higher occurrence rate of planets with 10 d< P < 50 d with respect to those with P < 10 d, albeit not as significant as that observed in our sample. Moreover, it is worth noting that the cumulative planet distribution as a function of the minimum mass (Fig. 6b) suggests that the planetmass distribution might show a valley between 3 and 5 M_{⊕}, where the cumulative frequency decreases steadily due to the increasing sensitivity of the survey. However, due to the limited number of detected planets, the uncertainties on the fine details of the planet distributions are still quite large.
As shown in Table 5, we derive an occurrence rate of lowmass (1 M_{⊕} < m_{p} sin i < 10 M_{⊕}) shortperiod (1 d < P < 10 d) planets around earlyM dwarfs of , which is lower (≃2σ) than the frequency derived by Bonfils et al. (2013) for laterM dwarfs, . This could indicate that earlyand latetype M dwarfs have different populations of lowmass shortperiod planets. Even taking into account all the candidate planets that we discussed in Sect. 3.4, our estimate of f_{occ} is still ≃1.5σ lower than that of Bonfils et al. (2013) (Table 6), which could still be evidence of this phenomenon. On the other hand, the same effect does not appear to be present for periods 10 d < P < 100 d, where our estimate of is perfectly compatible with the value obtained for laterM dwarfs, or even higher if all the candidates planets were confirmed (although still within the 1σ uncertainty of the value obtained by Bonfils et al. (2013) . It is worth noting that we find very large upper limits (>1) for lowmass planets with periods larger than 100 d, and this is in line with previous results (Tuomi et al. 2014): this is due to the low sensitivity of RV surveys to lowmass planets at these periods. However, it is interesting to notice that, due to this lack of observational constraints, lowmass planets at intermediate periods (100 d < P < 1000 d) could be very abundant, even if only four such planets have been detected between both RV and transit surveys to date^{8}.
Sabotta et al. (2021), analysing a subsample of the CARMENES survey, found slightly higher occurrence rates than Bonfils et al. (2013), although compatible within 1σ with their results in all period and minimummass bins in which both surveys had detected planets. It is interesting to notice that Sabotta et al. (2021) found a relatively high frequency of giant planets (100 M_{⊕} < m_{p} sin i < 1000 M_{⊕}) with periods up to 1000 d, , which is somewhat larger than our upper limit of f_{occ} < 0.02. However, as they mention in their work, this could be explained by the fact that the results of Sabotta et al. (2021) could be biased by the selection criteria of the analysed subsample. For intermediate masses (10 M_{⊕} < m_{p} sin i < 100 M_{⊕}), Sabotta et al. (2021) found increasing planetary frequencies towards longer periods up to 1000 d, which appears to be in contrast with our occurrence rates that peaks at intermediate periods (see Table 5); it is nevertheless worth noticing that our estimate of the frequencies of longperiod planets (see Sect. 6.3) is quite large, which could confirm the increasing trend observed by Sabotta et al. (2021). Looking at the frequency of lowmass shortperiod planets, Sabotta et al. (2021) found , which is higher than both our estimates and that from the HARPS Mdwarf survey (Bonfils et al. 2013). This again could be explained by the difference in spectral type between the three samples, as the subsample analysed by Sabotta et al. (2021) includes laterM dwarfs, which are not observed by the other surveys. Moreover, considering stars with masses higher and lower than 0.34 M_{⊕} separately, the CARMENES occurrence rates become f_{occ} < 0.22 and , respectively: while the former value is compatible with our estimate of , the latter is significantly higher (>3σ) than our value, which again suggests that the frequency of smallmass planets strongly increases with later spectral types.
6.2 Comparison with transiting planets
Comparing with the frequency of small planets around Kepler M dwarfs, it is evident that our derived occurrence rates for earlyM dwarfs are much lower than the average per star of 2.2 ± 0.3 planets with 1 R_{⊕} < R_{p} < 4 R_{⊕} and 1.5 d< P < 180 d derived by Gaidos et al. (2016).^{9} However, HardegreeUllman et al. (2019) found that the number of planets per star with orbital periods shorter than 10 d for Kepler midM dwarfs decreases from M5 V to M3 V stars. Assuming this effect also applies to intermediateperiod planets, this could mitigate the difference between our RVbased results and Kepler occurrence rates, as the sample analysed by Gaidos et al. (2016) includes many latertype M dwarfs which are not included in our survey. The results from HardegreeUllman et al. (2019) also confirm the discrepancy we observed between our earlyM sample and the latertype survey of Bonfils et al. (2013).
Moreover, Muirhead et al. (2015) found that of Kepler earlyM dwarfs host compact multiplanets systems with orbital periods lower than 10 d: while this frequency could be consistent with the frequency of shortperiod lowmass planetary systems we derive from our survey, namely , when also taking into account the candidate planets^{10}, we found only one system hosting two shortperiod planets (GJ 3998 Affer et al. 2016), while all other detected systems contained only a single planet. This apparent incompatibility could be mitigated by selection effects, because many multiple systems discovered by Kepler are composed of small 1−2 R_{⊕} planets, which could easily fall below our M_{⊕} detection threshold for planets with periods shorter than 10 d (Marcy et al. 2014). Moreover, many Kepler multiplanet systems host planets with very short periods below 2 d, which are inherently difficult to detect in groundbased surveys because of daily sampling limitations. Another possible explanation of the discrepancy between the number of multiplanet systems found in our RV survey and the statistics of Kepler planets could be found in the metallicity of our sample: recent studies (e.g. Anderson et al. 2021) found evidence that Kepler compact multiple planetary systems are hosted by M dwarfs that are statistically more metalpoor than the general population of systems with no detected transiting planets. This is coherent with the observed behaviour of Sunlike stars, which show increasing frequencies of compact multiple systems for metallicities [Fe/H] < −0.3 (Brewer et al. 2018). The complete HADES sample has a mean [Fe/H] = −0.13, and the planethosting subsample has a similar mean metallicity of [Fe/H] = −0.15 (see Fig. 8): as the M dwarfs in our sample are relatively metalrich compared to the bulk of compactsystemhosting Kepler stars, this could explain the scarcity of multiple planetary systems detected by our survey. This could also be in line with the results from the analysis in Maldonado et al. (2020), which showed that there might be an anticorrelation between the frequency of lowmass planets and stellar metallicity.
Fig. 8 Distribution of stellar metallicities for the HADES samples and for the subsample of planethosting stars (including also planetary candidates from Sect. 3.4). 
6.3 Longperiod and giant planets
From our analysis of longterm trends, in the assumption that all the trends listed in Table 3 are in fact of planetary origin, we can estimate the frequency of longperiod (P > 3 × 10^{3} d) planets. We estimate that, for intermediate masses (10 M_{⊕} < m_{p} sin i < 100 M_{⊕}), such planets could be quite abundant around earlyM dwarfs, , with a much higher frequency than the estimate for the same masses and periods 10^{3} d < P < 3 × 10^{3} d (Table 6). This result could be confirmed by the similar occurrence rate f_{occ} ≃ 0.20 computed by Tuomi et al. (2014) for planets with periods longer than 1000 d, even if their survey suffered from the same limitations discussed in Sect. 5.2 due to the similar timespan of the observations. Moreover, this would confirm the trend observed by Sabotta et al. (2021), who measured increasing intermediatemass planet frequencies at increasing periods. Therefore, our survey confirms the scarcity of giant planets (100 M_{⊕} < m_{p} sin i < 1000 M_{⊕}) at shorttointermediate periods around M dwarfs, setting an upper limit to their occurrence f_{occ} < 0.02 for periods P < 3000 d (Table 5). This is consistent with the behaviour observed for laterM dwarfs by Bonfils et al. (2013), who estimated the frequency of giant planets of periods up to 1000 d to be f_{occ} ≃ 0.02, even if their sample did in fact contain two such planets. Moreover, having detected two verylongperiod giant planets (GJ 849 b and GJ 832b), they estimated the frequency of giant planets with periods 10^{3} d < P < 10^{4} d to be . Even if the timespan of our survey did not allow for the detection of such longperiod planets, this result is perfectly compatible with our estimate of the frequency of longperiod giant planets derived from the observed RV trends (Sect. 5.2). However, the occurrence rate we derive is larger (≃2σ) than the estimated frequency of Jupiterlike planets around M dwarfs computed combining RV and microlensing surveys, (Clanton & Gaudi 2014). Moreover, it is worth noting that the f_{occ,t} we derived is very close to the recent results by Wittenmyer et al. (2020), who derived the frequency of cool Jupiters to be from an RV survey of FGK stars. This could suggest that the frequency of longperiod giant planets does not vary as strongly as the frequency of lowmass planets from Solartype stars to M dwarfs.
Comparing Fig. 1 and Table 2 to Table 3, it is interesting to notice that, apart from GJ 15A, none of our stars hosting shortperiod lowmass planets show a longperiod trend compatible with the presence of an outer planetary companion. This appears to be in contrast with theoretical models, as Izidoro et al. (2015) predicted that systems with only one lowmass planet with orbital period shorter than 100 d should also harbour a Jupiterlike planet on a more distant orbit: otherwise both in situ formation and inward migration models predict that superEarth and warm Neptunes should form in rich systems hosting many closein planets. However, of all of our detected planets and candidates, only one system appears to host more than a single closein planet (GJ 3998), and, as previously mentioned, the other systems do not show any evidence of the expected outer massive companions.
6.4 Habitable planets
The HADES program has not been able to confirm the presence of any planet in the HZ so far, and therefore even if the detection function was generally low for lowmass planets in the HZ (Fig. 7), we derived an upper limit to the frequency of habitable planets around earlyM dwarfs of η_{⊕} < 0.23 (see Sect. 5.1). Bonfils et al. (2013) estimated a higher frequency of habitable planets around M dwarfs . However, this value takes into account the presence of the habitable planet GJ 581 d, which was later discarded as an artifact of stellar activity (Baluev 2013; Robertson et al. 2014). Therefore, recomputing their value of η_{⊕} considering only the other habitable planet detected in their survey, the obtained value is lower η_{⊕} ≃ 0.20, which is much closer to the value derived from our analysis of the HADES survey. Moreover, Tuomi et al. (2014) found a similar value of considering planets with masses 3 M_{⊕} < m_{p} sin i < 10 M_{⊕}, which is compatible within 1σ with our estimate. The upper limit we derived, η_{⊕} < 0.23, could still allow a higher number of lowmass planets in the HZ of earlyM dwarfs compared to Gtype stars, for which the fraction of Earthlike planets is estimated to be η_{⊕} ≃ 0.10 (Perryman 2018). Instead, if we consider the weak candidate planet detected in this analysis orbiting in the outer rim of the HZ of its host star, our estimate of η_{⊕} increases to , which would be compatible with the previous estimates from RV Mdwarf surveys, and confirm the abundance of habitable planets around lowmass stars compared to Solartype stars.
7 Conclusions
We presented the statistical analysis of the spectroscopic observations of a magnitudelimited sample of nearby earlyM dwarfs observed within the HADES programme. The sample is composed of 56 targets, each with an average of 77 highprecision RVs. We analysed all of the timeseries in the sample with a uniform Bayesian technique in order to obtain consistent and unbiased estimates of the detection limits and planetary occurrence rates. Moreover, we applied GP regression to refine the planetary parameters and improve the detection efficiency around the most observed and active targets in the survey.
The sample includes ten published planetary systems discovered and characterised as part of the HADES programme, and in this work we present five new planetary candidates. Moreover, we discuss eight RV longterm trends that could be produced by longperiod planetary companions. The new candidates will be further analysed and discussed in future focused publications, two of which are in preparation and soon to be submitted (González Hernández et al. on GJ 21 and Affer et al. on GJ 3822). The other candidates, and in particular the weak HZ candidate discussed in Sect. 5.1, will need additional RV observations to confirm or disprove their planetary nature. Similarly, an indepth analysis of the longterm trends will be the focus of a future dedicated work, as a survey focused on the characterization of latetype systems hosting longperiod planets is currently being carried out within the GAPS programme (Barbato et al. 2020).
We confirm that giant planets (100 M_{⊕} < m_{p} sin i < 1000 M_{⊕}) are very rare around M dwarfs at shorttointermediate periods (P < 3000 d), f_{occ} < 0.02. On the other hand, lowmass planets (1 M_{⊕} < m_{p} sin i < 10 M_{⊕}) appear to be common, with frequencies varying between for short periods (1 d < P < 10 d) to for longer periods (10 d < P < 100 d). While high compared to Solartype stars, these frequencies appear to be lower than for laterM systems, therefore confirming the strong dependence of planetary occurrence rates on stellar mass. It is worth noting that early and lateM dwarfs have very different internal structures, as the latter are fully convective, and this affects the efficiency of tidal interaction between a star and a closeby planet (Barker 2020), which could explain the significant difference in planetary occurrence rates between early and lateM dwarfs. We also estimated the frequency of habitable planets around earlyM dwarfs, finding an upper limit of η_{⊕} < 0.23, which suggests that the frequency of Earthlike planets around such stars is as high as if not higher than around Solartype stars.
Finally, it is worth noting how the current number of HADES closein planets (11–16) is incompatible with the prediction of 3.8 ± 1.9 detections (Perger et al. 2017a) based on previous Mdwarf planet population models. This highlights the importance of longterm highprecision surveys focused on narrow intervals of stellar masses in improving our knowledge of planetary populations, and therefore formation mechanisms, through different classes of stellar hosts.
Acknowledgements
M.Pi., A.P., A.M., and A.S. acknowledge partial contribution from the agreement ASIINAF n.201816HH.0. M.Pi. acknowledges partial contribution from OB.FU. 1.05.02.85.13 “Planetary Systems At Early Ages (PLATEA)”. J.M., and G.M. acknowledge support from the accordo attuativo ASIINAF n.20215HH.0 “Partecipazione italiana alla fase B2/C della missione Ariel”. E.G.A. acknowledge support from the Spanish Ministery for Science, Innovation, and Universities through projects AYA201679425C31/2/3P, AYA201569350C32P, ESP201787676C52R, ESP201787143R. The Centro de Astrobiología (CAB, CSICINTA) is a Center of Excellence “Maria de Maeztu”. M.D. acknowledges financial support from the FP7SPACE Project ETAEARTH (GA no. 313014). P.G. gratefully acknowledges support from the Italian Space Agency (ASI) under contract 201824HH.0. M.Pe. and I.R. acknowledge support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grant PGC2018098153B C33, as well as the support of the Generalitat de Catalunya/CERCA programme. J.I.G.H. acknowledges financial support from Spanish Ministry of Science and Innovation (MICINN) under the 2013 Ramón y Cajal program RYC2013 14875. B.T.P. acknowledges Fundación La Caixa for the financial support received in the form of a Ph.D. contract. A.S.M. acknowledges financial support from the Spanish MICINN under the 2019 Juan de la Cierva Programme. B.T.P., A.S.M., J.I.G.H. and R.R. acknowledge financial support from the Spanish MICINN project PID2020117493GBI00, and from the Government of the Canary Islands project ProID2020010129. The HARPSN Project is a collaboration between the Astronomical Observatory of the Geneva University (lead), the CfA in Cambridge, the Universities of St. Andrews and Edinburgh, the Queen’ss University of Belfast, and the TNGINAF Observatory. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. We also thank the anonymous referee for the insightful review.
Appendix A MCMC analysis: Planetary parameters, GP hyperparameters, and adopted priors
We adopted uninformative priors for all the parameters and hyperparameters in our MCMC models. This was done for consistency across the statistical analysis, and to avoid adulterating the computed detection function by adopting incorrect informative priors. For most of the RV model parameters and GP hyperparameters we adopted a uniform prior, except for the correlation decay timescale, λ, for which we adopted a uniform prior in logarithmic scale to avoid oversampling the long scales, as λ can be of the order of tens or hundreds of days and more, depending on the activity level of the target. When an eccentric Keplerian signal was included in the model, instead of fitting separately e_{j} and ω_{j}, we use the auxiliary parameters · cos ω_{j} and · sin ω_{j} to reduce the covariance between e_{j} and ω_{j}, especially for low eccentricity value. Moreover, when the eccentricity was consistent with e_{j} = 0 we adopted a circular model, fixing e_{j} = ω_{j} = 0, to reduce the number of parameters and the computational weight of the analysis.
The analysis of GJ 15A required some different priors for the RV model because, as discussed in Pinamonti et al. (2018), the complete orbital solution was obtained by a combination of HIRES + HARPSN data, and was not recoverable from HARPSN data only. For this reason, we applied Gaussian priors to the GP hyperparameters, to the semiamplitude and period of the inner planets GJ 15A b, and to the acceleration coefficient d, which as discussed in Sect. 3.4 corresponds to the RV variations caused by the outer planet GJ 15A c over the timespan of the HARPSN observations.
Finally, particular care was required for the three parameters of the testplanet RV model, whose prior could directly impact the resulting function. For the orbital period P_{test} we adopted a uniform prior in logarithmic scale over the interval log , choosing the upper limit to be roughly twice the average timespan of the observations. The prior of RV amplitude K_{test} required a particularly large prior to avoid constraining the derived minimummass values over all the explored orbital period values: for this reason we adopted a very broad logarithmic uniform prior log , with an upper limit sufficiently larger than the maximum r.m.s. of the HADES timeseries (see Table 1). Lastly, to fit the reference time of the testplanet orbit, T_{0,test}, we defined the auxiliary parameter t_{p,test} = T_{c}/P_{test}, with T_{c} the time of inferior conjunction: t_{p,test} could be thus defined over a uniform prior , avoiding the multimodal distribution of T_{c}, which could greatly impact the efficiency of emcee (ForemanMackey et al. 2013).
Moreover, in Fig. A.1 is shown the detection function map of GJ399, with marked the position of the 32.9 signal identified in the RV timeseries which, as discussed in Sect. 5.1, is just below the acceptance threshold ΔBIC = 10: it is worth noticing that the weakcandidate signal mass, M_{p} sin i = 7.0 M_{⊕}, is just below the M_{⊕} detection threshold for periods in [20,40] d. This confirms the good correspondence between the adopted planetary acceptance criterion (Sect. 3.2.1) and the MCMCderived detectability function (Sect. 3.2).
Fig. A.1 Detection function map of the RV timeseries of GJ 399 (as in Fig. 2). The yellow circle shows the position of the 32.9 d weak candidate signal. 
Appendix B New planetary candidates
In Fig. B.1 to B.5 the relevant plots describing the planetary candidates detected in this analysis, as discussed in Sect. 3.4, are shown. The periodic signals in the RV data are shown via the Generalized Lomb Scargle periodogram analysis (GLS, Zechmeister & Kürster 2009) of the RV timeseries, after all stellar signals and the Base Model have been removed.
Fig. B.1 Planetary candidate in the GJ 21 system. a) GLS periodogram of the RV timeseries after subtracting the Base Model and the GP activity signal. The vertical blue dashed line marks the orbital period of the candidate, while the vertical red dotdashed and dotted lines show the stellar rotation period and its first harmonic, respectively. b) Phasefolded for the RV curve of the planetary candidate. The red solid line represents the bestfit Keplerian model, while the red circles indicates the binned RVs with the corresponding rms. 
Fig. B.2 Planetary candidate in the GJ 1074 system. a) GLS periodogram of the RV timeseries after subtracting the Base Model. The vertical blue dashed line marks the orbital period of the candidate, while the vertical red dotdashed line shows the stellar rotation period as derived by Suárez Μascareño et al. (2018). b) Phasefolded for the RV curve of the planetary candidate. The red solid line represents the bestfit Keplerian model, while the red circles indicates the binned RVs with the corresponding rms. 
Fig. B.3 Planetary candidate in the GJ 9404 system. a) GLS periodogram of the RV timeseries after subtracting the Base Model. The vertical blue dashed line marks the orbital period of the candidate. b) Phasefolded for the RV curve of the planetary candidate. The red solid line represents the bestfit Keplerian model, while the red circles indicates the binned RVs with the corresponding rms. 
Fig. B.4 Planetary candidate in the GJ 548A system. a) GLS periodogram of the RV timeseries after subtracting the Base Model. The vertical blue dashed line marks the orbital period of the candidate. b) Phasefolded for the RV curve of the planetary candidate. The red solid line represents the bestfit Keplerian model, while the red circles indicates the binned RVs with the corresponding rms. 
Fig. B.5 Planetary candidate in the GJ 3822 system. a) GLS periodogram of the RV timeseries after subtracting the Base Model. The vertical blue dashed line marks the orbital period of the candidate. b) Phasefolded for the RV curve of the planetary candidate. The red solid line represents the bestfit Keplerian model, while the red circles indicates the binned RVs with the corresponding rms. 
Appendix C MCMC analysis: bestfit RV models
The bestfit models of all the HADES targets fitted in this analysis are summarised in Table C.1. The bestfit model components are abbreviated as follow: GP bestfit model includes GP regression as described in Sect. 3.3; BM Base Model from Eq. 6; Quad quadratic component added to the Base Model as in Eq. 13; Kpl(ecc) full Keplerian component from Eq. 7; Kpl(circ) circular Keplerian component with e = ω = 0 in Eq. 7; Act stellar activity RV signal, modelled as a sine wave as discussed in Sect. 3.3..
Bestfit RV models of the targets of the survey.
References
 Adibekyan, V. Z., Sousa, S. G., Santos, N. C., et al. 2012, A&A, 545, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Affer, L., Micela, G., Damasso, M., et al. 2016, A&A, 593, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Affer, L., Damasso, M., Micela, G., et al. 2019, A&A, 622, A193 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ambikasaran, S., ForemanMackey, D., Greengard, L., Hogg, D. W., & O'Neil, M. 2015, IEEE Trans. Pattern Anal. Mach. Intell., 38, 252 [Google Scholar]
 Anderson, S. G., Dittmann, J. A., Ballard, S., & Bedell, M. 2021, AJ, 161, 203 [NASA ADS] [CrossRef] [Google Scholar]
 AngladaEscudé, G., & Butler, R. P. 2012, ApJS, 200, 15 [Google Scholar]
 AngladaEscudé, G., Tuomi, M., Gerlach, E., et al. 2013, A&A, 556, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AngladaEscudé, G., Tuomi, M., Arriagada, P., et al. 2016, ApJ, 830, 74 [CrossRef] [Google Scholar]
 AstudilloDefru, N., Forveille, T., Bonfils, X., et al. 2017, A&A, 602, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baluev, R. V. 2013, MNRAS, 429, 2052 [Google Scholar]
 Barbato, D., Pinamonti, M., Sozzetti, A., et al. 2020, A&A, 641, A68 [EDP Sciences] [Google Scholar]
 Barker, A. J. 2020, MNRAS, 498, 2270 [Google Scholar]
 Bonfils, X., Mayor, M., Delfosse, X., et al. 2007, A&A, 474, 293 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bonfils, X., Delfosse, X., Udry, S., et al. 2013, A&A, 549, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brewer, J. M., Wang, S., Fischer, D. A., & ForemanMackey, D. 2018, ApJ, 867, L3 [NASA ADS] [CrossRef] [Google Scholar]
 Clanton, C., & Gaudi, B. S. 2014, ApJ, 791, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Covino, E., Esposito, M., Barbieri, M., et al. 2013, A&A, 554, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cumming, A., & Dragomir, D. 2010, MNRAS, 401, 1029 [NASA ADS] [CrossRef] [Google Scholar]
 Cumming, A., Butler, R. P., Marcy, G. W., et al. 2008, PASP, 120, 531 [CrossRef] [Google Scholar]
 Damasso, M., Pinamonti, M., Scandariato, G., & Sozzetti, A. 2019, MNRAS, 489, 2555 [NASA ADS] [CrossRef] [Google Scholar]
 Dekker, H., D'Odorico, S., Kaufer, A., Delabre, B., & Kotzlowski, H. 2000, in Proc. SPIE, 4008, Optical and IR Telescope Instrumentation and Detectors, eds. M. Iye, & A. F. Moorwood, 534 [NASA ADS] [CrossRef] [Google Scholar]
 Desidera, S., Sozzetti, A., Bonomo, A. S., et al. 2013, A&A, 554, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dressing, C. D., & Charbonneau, D. 2013, ApJ, 767, 95 [Google Scholar]
 Dressing, C. D., & Charbonneau, D. 2015, ApJ, 807, 45 [Google Scholar]
 Dumusque, X., Lovis, C., Udry, S., & Santos, N. C. 2011, in IAU Symposium, 276, The Astrophysics of Planetary Systems: Formation, Structure, and Dynamical Evolution, eds. A. Sozzetti, M. G. Lattanzi, & A. P. Boss, 530 [NASA ADS] [Google Scholar]
 Dumusque, X., Glenday, A., Phillips, D. F., et al. 2015, ApJ, 814, L21 [Google Scholar]
 Eastman, J., Gaudi, B. S., & Agol, E. 2013, PASP, 125, 83 [Google Scholar]
 Endl, M., Kürster, M., Els, S., et al. 2002, A&A, 392, 671 [CrossRef] [EDP Sciences] [Google Scholar]
 Endl, M., Cochran, W. D., Tull, R. G., & MacQueen, P. J. 2003, AJ, 126, 3099 [Google Scholar]
 Endl, M., Cochran, W. D., Kürster, M., et al. 2006, ApJ, 649, 436 [Google Scholar]
 Feng, F., Butler, R. P., Shectman, S. A., et al. 2020, ApJS, 246, 11 [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
 Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaidos, E., Mann, A. W., Kraus, A. L., & Ireland, M. 2016, MNRAS, 457, 2877 [Google Scholar]
 Giacobbe, P., Benedetto, M., Damasso, M., et al. 2020, MNRAS, 491, 5216 [Google Scholar]
 Giles, H. A. C., Collier Cameron, A., & Haywood, R. D. 2017, MNRAS, 472, 1618 [Google Scholar]
 Gomes da Silva, J., Santos, N. C., Bonfils, X., et al. 2011, A&A, 534, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 GonzálezÁlvarez, E., Micela, G., Maldonado, J., et al. 2019, A&A, 624, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 GonzálezÁlvarez, E., Petralia, A., Micela, G., et al. 2021, A&A, 649, A157 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 HardegreeUllman, K. K., Cushing, M. C., Muirhead, P. S., & Christiansen, J. L. 2019, AJ, 158, 75 [Google Scholar]
 Hawley, S. L., Davenport, J. R. A., Kowalski, A. F., et al. 2014, ApJ, 797, 121 [Google Scholar]
 Ida, S., & Lin, D. N. C. 2005, ApJ, 626, 1045 [NASA ADS] [CrossRef] [Google Scholar]
 Izidoro, A., Raymond, S. N., Morbidelli, A., Hersant, F., & Pierens, A. 2015, ApJ, 800, L22 [Google Scholar]
 Kass, R. E., & Raftery, A. E. 1995, J. Am. Stat. Assoc., 90, 773 [Google Scholar]
 Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108 [Google Scholar]
 Kervella, P., Arenou, F., & Thévenin, F. 2022, A&A, 657, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kipping, D. M., Hartman, J., Bakos, G. Ä., et al. 2011, AJ, 142, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013a, ApJ, 770, 82 [NASA ADS] [CrossRef] [Google Scholar]
 Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013b, ApJ, 765, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Lépine, S., & Gaidos, E. 2013, Astron. Nachr., 334, 176 [CrossRef] [Google Scholar]
 Lindegren, L., Klioner, S. A., Hernandez, J., et al. 2021, A&A, 649, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lovis, C., Dumusque, X., Santos, N. C., et al. 2011, ArXiv eprints [arXiv:1187.5325] [Google Scholar]
 Lovis, C., & Pepe, F. 2007, A&A, 468, 1115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Luque, R., Nowak, G., Pallé, E., et al. 2018, A&A, 620, A171 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maldonado, J., Affer, L., Micela, G., et al. 2015, A&A, 577, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maldonado, J., Scandariato, G., Stelzer, B., et al. 2017, A&A, 598, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maldonado, J., Micela, G., Baratella, M., et al. 2020, A&A, 644, A68 [EDP Sciences] [Google Scholar]
 Maldonado, J., Petralia, A., Damasso, M., et al. 2021, A&A, 651, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marcy, G. W., Isaacson, H., Howard, A. W., et al. 2014, ApJS, 210, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20 [NASA ADS] [Google Scholar]
 Mordasini, C., Alibert, Y., Benz, W., Klahr, H., & Henning, T. 2012, A&A, 541, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Muirhead, P. S., Mann, A. W., Vanderburg, A., et al. 2015, ApJ, 801, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Newton, E. R., Irwin, J., Charbonneau, D., BertaThompson, Z. K., & Dittmann, J. A. 2016, ApJ, 821, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Perger, M., GarciaPiquer, A., Ribas, I., et al. 2017a, A&A, 598, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perger, M., Ribas, I., Damasso, M., et al. 2017b, A&A, 608, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perger, M., Scandariato, G., Ribas, I., et al. 2019, A&A, 624, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perryman, M. 2018, The Exoplanet Handbook (Cambridge University Press) [Google Scholar]
 Pinamonti, M., Damasso, M., Marzari, F., et al. 2018, A&A, 617, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pinamonti, M., Sozzetti, A., Giacobbe, P., et al. 2019, A&A, 625, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Poretti, E., Boccato, C., Claudi, R., et al. 2016, Mem. Soc. Astron. Italiana, 87, 141 [NASA ADS] [Google Scholar]
 Reid, I. N., Hawley, S. L., & Gizis, J. E. 1995, AJ, 110, 1838 [Google Scholar]
 Robertson, P., Mahadevan, S., Endl, M., & Roy, A. 2014, Science, 345, 440 [Google Scholar]
 Sabotta, S., Schlecker, M., Chaturvedi, P., et al. 2021, A&A, 653, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Scandariato, G., Maldonado, J., Affer, L., et al. 2017, A&A, 598, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schwarz, G., et al. 1978, Ann. Stat., 6, 461 [Google Scholar]
 Selsis, F., Kasting, J. F., Levrard, B., et al. 2007, A&A, 476, 1373 [CrossRef] [EDP Sciences] [Google Scholar]
 Sousa, S. G., Santos, N. C., Mayor, M., et al. 2008, A&A, 487, 373 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sozzetti, A., Bernagozzi, A., Bertolini, E., et al. 2013, in European Physical Journal Web of Conferences, 47, 03006 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Suárez Mascareño, A., Gonzalez Hernandez, J. I., Rebolo, R., et al. 2017, A&A, 605, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Suárez Mascareño, A., Rebolo, R., Gonzalez Hernandez, J. I., et al. 2018, A&A, 612, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 ToledoPadrón, B., Suarez Mascareno, A., Gonzalez Hernandez, J. I., et al. 2021, A&A, 648, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tuomi, M., Jones, H. R. A., Barnes, J. R., AngladaEscudé, G., & Jenkins, J. S. 2014, MNRAS, 441, 1545 [Google Scholar]
 Vanderburg, A., Plavchan, P., Johnson, J. A., et al. 2016, MNRAS, 459, 3565 [NASA ADS] [CrossRef] [Google Scholar]
 Wittenmyer, R. A., Wang, S., Horner, J., et al. 2020, MNRAS, 492, 377 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, H., Liu, J., Gao, Q., et al. 2017, ApJ, 849, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577 [CrossRef] [EDP Sciences] [Google Scholar]
The MCMCderived detectability function we adopt in this analysis shows good correspondence with the ΔBIC = 10 detection threshold, as shown in Appendix A.
This is a good approximation over the HARPSN timeseries, as can be seen in Fig. 6 of Pinamonti et al. (2018).
NASA exoplanet archive  08/04/2021: GJ 180 d P = 106.300 ± 0.129d, m_{p} sin i = 7.56 ± 1.07 M_{⊕}, GJ 229Ac P = 121.995 ± 0.161 d, m_{p} sin i = 7.268 ± 1.256 M_{⊕} (Feng et al. 2020); GJ 628 d , (AstudilloDefru et al. 2017); GJ 667Cg , (AngladaEscudé et al. 2013).
Other studies of Kepler M dwarfs, found lower number of planets per star, ~1.5 for 1 R_{⊕} < R_{p} < 4 R_{⊕} and 1. d < P < 10 (Sabotta et al. 2021, and references therein), but these numbers are still much higher than the occurrence rates we derive from the HADES RV sample.
Assuming the empyrical planet mass—radius relationships of Marcy et al. (2014) adopted by Muirhead et al. (2015).
Bonfils et al. (2013) used the definition of HZ from Selsis et al. (2007), while Tuomi et al. (2014) used the same updated definition we adopted in our work (Kopparapu et al. 2013b).
All Tables
Stellar parameters and summary of the observations of the 56 HADES targets analysed in this work.
Longperiod candidate planetary signals assumed from the modelled linear RV trends.
All Figures
Fig. 1 Overview of the HADES detected planetary systems. The published planets of the sample are shown as red circles: the symbol size is proportional to the minimum planetary mass. The conservative and optimistic limits of the habitable zone of each system (see Sect. 5.1) are shown as thick dark green and light green bands, respectively. 

In the text 
Fig. 2 Detection function map of the RV timeseries of GJ 15 A. The white part corresponds to the area in the period—minimum mass space where additional signals could be detected if present in the data, i.e. p_{i} = 1, while the black region corresponds to the area where the detection probability is negligible, i.e. p_{i} = 0. The red circle marks the position in the parameter space of the planet GJ 15 A b. 

In the text 
Fig. 3 Global HADES detection map. The color scale expresses the global detection function, p. The red circles mark the position in the parameter space of the confirmed HADES planets (see Fig. 1), while the yellow circles indicate the additional candidates presented in this study (see Table 2). 

In the text 
Fig. 4 Averaged HADES detection map. The color scale expresses the global detection function, p, averaged over the intervals Δ(P, M) defined for the computation of the occurrence rates. The red circles mark the positions in the parameter space of the confirmed HADES planets (see Fig. 1), while the yellow circles indicate the additional candidates presented in this study (see Table 2). 

In the text 
Fig. 5 Planetary occurrence rate f_{occ} as a function of the orbital period, with the corresponding 1σ uncertainties (top panel). The median cumulative planet frequency is shown in the lower panel. The red and yellow distributions correspond to the occurrence rates derived from the confirmedonly, k_{p}, and confirmed+candidates planets, k_{c}, respectively. 

In the text 
Fig. 6 Planetary occurrence rate f_{occ} as a function of the minimum mass, with the corresponding 1 σ uncertainties (top panel). The median cumulative planet frequency is shown in the lower panel. The red and yellow distributions correspond to the occurrence rates derived from the confirmedonly, k_{p}, and confirmed + candidates planets, k_{c}, respectively. 

In the text 
Fig. 7 Habitable zone HADES detection map, derived as described in Sect. 5.1. The colour scale expresses the detection function, as in Fig. 3. 

In the text 
Fig. 8 Distribution of stellar metallicities for the HADES samples and for the subsample of planethosting stars (including also planetary candidates from Sect. 3.4). 

In the text 
Fig. A.1 Detection function map of the RV timeseries of GJ 399 (as in Fig. 2). The yellow circle shows the position of the 32.9 d weak candidate signal. 

In the text 
Fig. B.1 Planetary candidate in the GJ 21 system. a) GLS periodogram of the RV timeseries after subtracting the Base Model and the GP activity signal. The vertical blue dashed line marks the orbital period of the candidate, while the vertical red dotdashed and dotted lines show the stellar rotation period and its first harmonic, respectively. b) Phasefolded for the RV curve of the planetary candidate. The red solid line represents the bestfit Keplerian model, while the red circles indicates the binned RVs with the corresponding rms. 

In the text 
Fig. B.2 Planetary candidate in the GJ 1074 system. a) GLS periodogram of the RV timeseries after subtracting the Base Model. The vertical blue dashed line marks the orbital period of the candidate, while the vertical red dotdashed line shows the stellar rotation period as derived by Suárez Μascareño et al. (2018). b) Phasefolded for the RV curve of the planetary candidate. The red solid line represents the bestfit Keplerian model, while the red circles indicates the binned RVs with the corresponding rms. 

In the text 
Fig. B.3 Planetary candidate in the GJ 9404 system. a) GLS periodogram of the RV timeseries after subtracting the Base Model. The vertical blue dashed line marks the orbital period of the candidate. b) Phasefolded for the RV curve of the planetary candidate. The red solid line represents the bestfit Keplerian model, while the red circles indicates the binned RVs with the corresponding rms. 

In the text 
Fig. B.4 Planetary candidate in the GJ 548A system. a) GLS periodogram of the RV timeseries after subtracting the Base Model. The vertical blue dashed line marks the orbital period of the candidate. b) Phasefolded for the RV curve of the planetary candidate. The red solid line represents the bestfit Keplerian model, while the red circles indicates the binned RVs with the corresponding rms. 

In the text 
Fig. B.5 Planetary candidate in the GJ 3822 system. a) GLS periodogram of the RV timeseries after subtracting the Base Model. The vertical blue dashed line marks the orbital period of the candidate. b) Phasefolded for the RV curve of the planetary candidate. The red solid line represents the bestfit Keplerian model, while the red circles indicates the binned RVs with the corresponding rms. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.