Mid-to-late M Dwarfs Lack Jupiter Analogs

We extract the spectra using the standard TRES (Buchhave et al. 2010) and CHIRON (Tokovinin et al. 2013) pipelines and measure radial velocities using an updated version of the method presented in Winters et al. (2020), which itself is based on Kurtz & Mink (1998). In addition to utilizing the TiO bandhead features from 7065 to 7165 Å (TRES echelle order 41 and CHIRON echelle order 44), we also consider five additional red echelle orders—36, 38, 39, 43, and 45 for TRES and 36, 37, 39, 40, and 51 for CHIRON, which represent wavelengths with low telluric contamination and high information content for low-mass M dwarfs. These orders fall within 6400–7850 Å.

We obtain an initial radial velocity estimate for each spectrum by cross-correlating with a template of Barnard's Star in the manner described in Winters et al. (2020). This method also produces an estimate of v sin i. We then shift and stack our observations to create a high-S/N template spectrum of the average slowly rotating, low-mass M dwarf observed by TRES. To stack the observations, we first create an error-weighted mean spectrum for each star. We then median combine the spectra of all 122 stars that were observed with TRES. Prior to coaddition, we mask regions with telluric depths greater than 2% in a representative TAPAS spectrum (Bertaux et al. 2014) or skyline emission greater than 5 × 10⁻¹⁷ erg s⁻¹ cm⁻² Å⁻¹ in the UVES sky emission atlas (Hanuschik 2003). Our final template neglects wavelengths at the edges of orders that do not include information from at least ten M dwarfs.

In creating this average low-mass M dwarf template, we note that the information content Q (Bouchy et al. 2001) in TRES echelle order 45 differs by a factor of two across our sample, with the highest values obtained for cold, metal-rich stars and the lowest values obtained for hot, metal-poor stars. To minimize radial velocity uncertainties caused by template mismatch, we therefore create six different average templates for TRES. Each star is initially classified by the Q we measure in order 45, dividing the sample equally into six preliminary bins. We then reclassify each star by determining the preliminary template that maximizes the cross-correlation coefficient, recreating the templates based on these refined classifications, and repeating until the classifications converge. This iterative process prevents stars from being misclassified in cases where their Q is not representative of their spectral type (for example, when noise leads to inflated estimates of Q). We perform a similar analysis for CHIRON, although we create only three templates so that we attain an acceptable S/N given the smaller number of observations available for coaddition.

We compute the cross-correlation function as the weighted sum:

$$\begin{eqnarray}&&{r}_{j}=\displaystyle \frac{\displaystyle \sum _{i}{a}_{i+j}{b}_{i}{w}_{i+j}}{\sqrt{\left(\displaystyle \sum _{i}{a}_{i+j}^{2}{w}_{i+j})(\displaystyle \sum _{i}{b}_{i}^{2}{w}_{i+j}\right)}},\end{eqnarray} \tag{ 1 }$$

where a is the observed flux, b is the template flux, and w are the variance weights (or w_i = 0 for masked telluric or skyline features). This sum is evaluated jointly across the six echelle orders. Both a and b are blaze corrected, normalized, continuum subtracted, and oversampled by a factor of 32. Our code also includes the ability to broaden b rotationally to analyze rapidly rotating stars, but there are no such stars in this study due to our Hα activity cut. The errors used to calculate w consider photon noise and read noise in a, and are scaled based on the blaze correction.

Lastly, we determine the radial velocity shift by fitting a Lorentzian to the cross-correlation peak.

We calculate theoretical radial velocity uncertainties using the equation derived in Bouchy et al. (2001):

$$\begin{eqnarray}&&\displaystyle \frac{{\delta }_{\mathrm{RV}}}{c}=\displaystyle \frac{1}{\sqrt{\displaystyle \sum _{i}{\lambda }_{i}^{2}{(\partial {b}_{i}/\partial {\lambda }_{i})}^{2}{w}_{i}}},\end{eqnarray} \tag{ 2 }$$

where b is once again the template flux and w are the variance weights. In the photon-noise limit, this equation simplifies to δ_RV = c/(Q × S/N).

Alternatively, we can estimate the radial velocity uncertainty directly from the cross-correlation function using the equation from Zucker (2003):

$$\begin{eqnarray}&&\displaystyle \frac{{\delta }_{\mathrm{RV}}}{c}=\sqrt{-{\left({N}_{\mathrm{eff}}\displaystyle \frac{r{{\prime\prime} }_{\,\max }}{{r}_{\max }}\displaystyle \frac{{r}_{\max }^{2}}{1-{r}_{\max }^{2}}\right)}^{-1}},\end{eqnarray} \tag{ 3 }$$

where we have replaced N with ${N}_{\mathrm{eff}}\equiv {({\sum }_{i}{w}_{i})}^{2}/{\sum }_{i}{w}_{i}^{2}$ to account for our error weighting. In this equation, ${r}_{\max }$ denotes the cross-correlation peak and $r{{\prime\prime} }_{\max }$ represents its second derivative.

Generally, we find good agreement between these two uncertainty estimators, although they differ when template mismatch is a significant source of uncertainty (e.g., in the high-S/N limit or when there is substantial rotational broadening, although the latter is not relevant for the sample discussed here). We make the conservative choice to always adopt the larger of the two.

To validate our uncertainties, we also compute the cross-correlation function for each of the six echelle orders individually and compare the resultant radial velocities using the chi-squared estimator:

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{k=1}^{6}{w}_{k}{({v}_{\mathrm{rad},k}-{v}_{\mathrm{rad}})}^{2},\end{eqnarray} \tag{ 4 }$$

where the variance weights (${w}_{k}=1/{\delta }_{k}^{2}$) are determined using the above uncertainty calculations for each order. For both TRES and CHIRON, the resulting distribution of χ² across all observations shows good agreement with a chi-squared distribution with five degrees of freedom, as we would expect for properly estimated uncertainties.

An independent source of uncertainty comes from the barycentric correction, which ideally would be calculated at the photon-weighted midpoint of the observation. As we lack information on the temporal flux distribution over an exposure, we adopt the geometric midpoint as the observation time. This approximation can produce a systematic error of up to 2 m s⁻¹ per minute of difference between the geometric and photon-weighted midpoints, as shown in Tronsgaard et al. (2019; see their Equation (8)). That work found that the offset between these midpoints is typically around 5% of the exposure length. For completeness, we therefore add a radial velocity uncertainty term in quadrature corresponding to a 5% offset to the midpoint time; however, we find that this contribution is ultimately negligible in the total error budget.

We note that our uncertainty estimates thus far do not account for the long-term instability of the spectrograph. Over the observation campaign, one zero-point offset occurred for CHIRON and two for TRES due to hardware and software changes to the instrument. To determine the amplitudes of these offsets, we create an ensemble of all our observations from each spectrograph, subtracting the error-weighted mean radial velocity of each star to put the observations on a common zero-point. We then use these ensembles to measure the amplitude of the radial velocity discontinuities and apply a correction to subsequent measurements. For CHIRON, this offset is +50 m s⁻¹ at BJD = 2,458,840 days. For TRES, these offsets are +45 m s⁻¹ at BJD = 2,458,700 days and −20 m s⁻¹ at BJD = 2,459,218 days.

To evaluate the stability of each spectrograph outside of these events, we consider the standard deviation of the most precise observations—namely, observations with uncertainty estimates between 5 and 18 m s⁻¹. If our observations were unaffected by a noise floor, we would expect the standard deviation to fall somewhere within this 5–18 m s⁻¹ range. Note that the delta degrees of freedom (ddof) in this calculation is not 1; this would underestimate the standard deviation, as we have subtracted the error-weighted mean of each star. The ddof is also not equal to the number of stars; this would overestimate the standard deviation, as all observations are used to inform the error-weighted means, not just observations with small estimated uncertainties. We make a more reasonable estimate of the ddof by assuming that each star contributes ∑w_j /∑w_i , where w_j are the variance weights of the observations with estimated uncertainties below 18 m s⁻¹ and w_i are the variance weights of all observations. The maximum value of this expression is 1, which occurs when the error-weighted mean radial velocity is fully informed by observations with small estimated uncertainties. For CHIRON, our calculation includes 137 observations, has ddof = 47.4, and results in a standard deviation of 20 m s⁻¹. For TRES, we have 48 observations with small estimated uncertainties, ddof = 16.2, and a standard deviation of 21 m s⁻¹. As these results are greater than 18 m s⁻¹, it appears that instrumental instability is indeed setting a noise floor. These results are consistent with Winters et al. (2020), who found typical rms velocity residuals of 20 m s⁻¹ using orbital solutions for bright, slowly rotating low-mass M dwarf binaries observed by TRES. As uncertainties below 20 m s⁻¹ may therefore be underestimated due to instrumental instability, we do not allow our estimates to be smaller than this floor.

Our radial velocity measurements and uncertainties for all 200 stars are available in a machine-readable format, with the format of this file shown in Table 1.

We use the metric P(χ²) to evaluate whether each radial velocity time series is consistent with an unvarying model given the uncertainties determined in the previous section. In this context, the chi-squared estimator is given by:

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i=1}^{4}\displaystyle \frac{{({v}_{\mathrm{rad},i}-\langle {v}_{\mathrm{rad}}\rangle )}^{2}}{{\delta }_{\mathrm{RV},i}^{2}},\end{eqnarray} \tag{ 5 }$$

with N_obs − 1 = 3 degrees of freedom. 〈v_rad〉 denotes the error-weighted mean of the radial velocity time series. We flag a signal as significant if the variability is inconsistent with the null hypothesis with 99% confidence (P(χ²) < 1%). In other words, we are searching for statistically significant excess radial velocity jitter; if such jitter is detected, further follow up is necessary to verify that this jitter is due to a planet and not an astrophysical or statistical false positive. Given our sample size of 200 stars, we expect an average of two statistical false positives due to random chance.

We conducted a preliminary analysis of our data in early 2021, from which we identified two stars—LHS 2899 and G 125-34—as candidate variables based on this P(χ²) < 1% cut. Based on the significance of the variability in the initial four observations, we proceeded to collect an additional 41 observations of LHS 2899 and 46 observations of G 125-34 with TRES between 2021 April and 2022 May. At the time of the preliminary analysis, we had yet to obtain a fourth observation for some stars in the sample; our follow up of LHS 2899 and G 125-34 therefore occurred concurrently with our main observation campaign, which was completed on 2022 May 11. When analyzing observations of all stars after data collection was completed, we found that the stability of the TRES zero-point worsened in 2021–2022, motivating us to augment our 20 m s⁻¹ noise floor with a dynamic floor using the radial velocity jitter in the TRES standard stars over each observing run (S. Quinn, private communication). These standards are quiet, Sun-like stars that are known to have very low intrinsic radial velocity variation based on observations from more precise instruments. We use these standards to determine the rms scatter in each observing run, which we add to our uncertainties in quadrature. In theory, we should also be able to correct for zero-point offsets using these standards; however, they are Sun-like stars and analyzed in bluer orders than our M dwarf spectra, and we find that the offsets in the M dwarfs are not entirely commensurate with those seen in the standards.

With this refinement to our uncertainty estimation, the P(χ²) in the initial four observations of LHS 2899 increased from <1% to 3%, no longer meeting our follow-up criterion. For the purposes of our statistical study, it is therefore not necessary to establish whether the variability of LHS 2899 is a planet or a false positive; however, for completeness we will nonetheless discuss it here. G 125-34 retained its P(χ²) < 1% designation and no other signals became significant with this change. P(χ²) values for all stars are given in Table 2.

For both G 125-34 and LHS 2899, the significance of the candidate variability increased with the TRES follow-up observations, with a final P(χ²) = 4.3 × 10⁻⁵ for G 125-34 and P(χ²) = 1.0 × 10⁻⁵ for LHS 2899. A periodogram analysis yielded candidate periods for orbital solutions, although the precision and quantity of the data were insufficient to prove or disprove the presence of a planet definitively. The preferred planetary solutions were a 0.22 M_J planet on a 77 day orbit for G 125-34 and a 0.70 M_J planet on a 337 day orbit for LHS 2899.

To reach definitive conclusions about the dispositions of these candidates, we obtained five observations of LHS 2899 and 13 observations of G 125-34 using the MAROON-X spectrograph (Seifahrt et al. 2018), an extreme precision radial velocity instrument on the 8.1 m Gemini-North telescope, over a two month interval. These observations were reduced by the MAROON-X team using the SERVAL pipeline (Zechmeister et al. 2018) and are shown in Figure 2. The final observation of LHS 2899 and the final three observations of G 125-34 were taken in a different MAROON-X observing run from the prior observations, with a zero-point offset of −1.5 ± 1.0 m s⁻¹ for the blue arm and +1.5 ± 1.0 m s⁻¹ for the red arm; the uncertainty in this offset is responsible for the larger error bars for these later observations. The radial velocity time series for both stars are flat at the 1 m s⁻¹ level, definitively ruling out all orbital solutions consistent with the TRES data. LHS 2899 and G 125-34 are statistical false positives, with their significance in the TRES follow-up observations likely the result of non-Gaussian outliers. There are no giant planet detections in our 200 star sample.

Zoom In Zoom Out Reset image size

For a given star, j, the probability that we do not detect a planet is P_null(μ)_j = 1 − μ × C_j , where μ is the planetary occurrence rate in a given bin and C is the survey completeness in that bin. The probability that we do not detect any planets around any of our stars is then the product:

$$\begin{eqnarray}&&{P}_{\mathrm{null}}(\mu )=\displaystyle \prod _{j}(1-\mu {C}_{j}).\end{eqnarray} \tag{ 6 }$$

To determine a 95% confidence limit on the occurrence rate, we numerically solve for μ given that P_null(μ) = 0.05. While we use this form of the equation to generate our Figure 3, note that in the limit where the binomial approximation ${(1-\mu )}^{{C}_{j}}\approx 1-\mu {C}_{j}$ applies (i.e., where μ < 1 and μ C_j ≪ 1, which is appropriate for this study), this equation can be rewritten as ${P}_{\mathrm{null}}{(\mu )=(1-\mu )}^{{N}_{\mathrm{eff}}}$, where N_eff, the effective sample size, is given by N_eff = ∑_j C_j . This form is a convenient simplification, as it only requires the average completeness instead of the completeness for each star in the sample, and it motivates our use of N_eff as the color bar axis in this figure.

Zoom In Zoom Out Reset image size

We conduct an injection and recovery analysis to determine C as a function of M_P sin i and orbital period. We derive our occurrence rate constraints under the assumption of circular orbits, following the precedent of past works (e.g., Howard et al. 2010b; Mayor et al. 2011). This assumption has been shown to be reasonable for e < 0.5 (Endl et al. 2002; Cumming & Dragomir 2010). While many eccentric cold Jupiters meet this criterion (Buchhave et al. 2018), more highly eccentric giant planets are known (e.g., Robertson et al. 2012; Moutou et al. 2015). Notably, giant planet formation via gravitational instability may produce large orbital eccentricities (Jennings & Chiang 2021). We consider the limitations of the circular orbit assumption in the following section.

We consider 100 periods spaced logarithmically between 10⁻¹ and 10⁵ days and 100 values of M_P sin i spaced logarithmically between 10⁻¹ and 10¹ M_J. For each mass and period, we generate 100 artificial radial velocity time series for each star using the observation times and radial velocity uncertainties of the real data, drawing phases from a random uniform distribution. For each hypothetical planet, every star is assigned a completeness C between 0 and 1, indicating the fraction of trials in which the simulated observations were variable at the P(χ²) < 1% level. The sample size N_eff is the sum of these fractions. This completeness assumes that we would gather sufficient MAROON-X follow-up observations to verify that any signal flagged at P(χ²) < 1% significance is a true planet.

Our injection and recovery results are shown in Figure 3, which also includes the period axis recast in terms of the zero-albedo equilibrium temperature, ${T}_{\mathrm{eq}}={[L/(16\pi \sigma {a}^{2})]}^{1/4}$, and instellation, $S=4\sigma {T}_{\mathrm{eq}}^{4}$. This transformation requires the masses and luminosities of each star; we adopt the masses from Winters et al. (2021), which are based on the K-band relation from Benedict et al. (2016), and estimate the bolometric luminosities from the photometry collated in Winters et al. (2021) and distances from Gaia parallaxes (Gaia Collaboration et al. 2018, 2021; Lindegren et al. 2021). In particular, we consider the bolometric corrections from both Pecaut & Mamajek (2013) and Mann et al. (2015), adopting the average of the luminosities calculated from these two methods. These masses and luminosities are included in Table 2.

There are two regimes in Figure 3, with a transition around T_eq ≈ 50 K. At short periods, the sensitivity is determined by the radial velocity semiamplitude of the orbit. The slope of the transition region between detectable and undetectable planets represents the interplay between the radial velocity semiamplitude and typical radial velocity uncertainties. At long periods, the sensitivity drops off rapidly due to the finite length of the observing campaign.

To quantify the impact of the circular orbit assumption, we repeat the analysis of the previous section but under the assumption that all planets have e = 0.5 in one trial or e = 0.8 in another. For each simulated planet, we draw the argument of periastron from a random uniform distribution. Figure 4 shows the change in sensitivities compared to the circular orbits case in terms of ΔN_eff, the change in the number of stars around which the planet is detectable. We find that ΔN_eff is generally small; the sample size is never decreased by more than 15 stars in the e = 0.5 trial and 39 stars in the e = 0.8 trial. There are other regions of parameter space in which the sample size is increased by similar amounts. Moreover, the percent decrease in effective sample size does not exceed 13% or 30%, respectively. Even in the limit where all giant planets are highly eccentric, the circular orbit assumption will produce reasonable occurrence rate constraints.

Zoom In Zoom Out Reset image size

Compared to a planet on a circular orbit, an eccentric planet achieves a higher radial velocity semiamplitude but spends more time at relative radial velocities near zero. As we approach the sensitivity limit from the detectable side, the decreased likelihood of observing the eccentric orbit in a peak/trough renders a previously detectable planet undetectable. As we approach the sensitivity limit from the undetectable side, the increased amplitude of the peaks/troughs renders a previously undetectable planet detectable.

To summarize the above reduction and analysis, we do not detect any giant planets in our sample of 200 low-mass, inactive M dwarfs. More specifically, our four measurements for each star do not vary in excess of our P(χ²) = 1% detection limit for 198 stars, and we initially flagged the remaining two stars as candidate radial velocity variables. We collected additional observations of these two candidates with TRES, but were unable to establish the provenance of the signals conclusively. We then obtained observations of the two stars using MAROON-X. The MAROON-X radial velocities show no variation at the 1 m s⁻¹ level over a two month interval, indicating that the candidate variability from the TRES observations are statistical false positives. After refuting these candidates, none of our stars exhibit radial velocity variability that exceeds our detection threshold. Nevertheless, our injection and recovery analysis indicates that we are very sensitive to a variety of giant planets (Figure 3). We can therefore place strong upper limits on the occurrence rate of these planets, as we discuss in the following section. The median mass of stars in our sample is 0.18 M_⊙, with a median radial velocity precision of 29 m s⁻¹ and a median observation baseline of 3.1 yr.

As shown in Figure 3, we have nearly complete sensitivity to M_P sin i > 1 M_J planets out to the water snow line. In regions of complete sensitivity, we place a 68% confidence upper limit of 0.57% and a 95% confidence upper limit of 1.5%. Jupiter- and super-Jupiter-mass planets interior to the snow line of low-mass M dwarfs are therefore exceedingly rare.

How does this compare to studies of all types of M dwarfs? The Bonfils et al. (2013) sample contains one planet in this regime: Gliese 876 b (Marcy et al. 1998; Delfosse et al. 1998). While that work published their survey sensitivity in specific bins that differ from the one currently under discussion (their Table 11), inspection of their Figure 15 indicates that they are nearly complete across this parameter range, with N_eff equal to roughly 96 stars. If the true occurrence rate of such planets around all M dwarfs was equal to our 68% confidence upper limit of 0.57%, binomial probability yields a 32% chance that Bonfils et al. (2013) would detect one planet in this bin, with a 58% chance of detecting zero and a 10% chance of detecting two or more planets. Our results are therefore in reasonable agreement with Bonfils et al. (2013) without invoking a decrease in warm Jupiter occurrence around low-mass M dwarfs relative to all M dwarfs, although our occurrence constraints are tighter given that our sample is twice as large. We also note that Gliese 876 has a mass of 0.37 M_⊙; i.e., while it is more massive than the M dwarfs in our sample, it is not a particularly massive M dwarf. That said, a decrease in warm Jupiter occurrence around low-mass M dwarfs relative to all M dwarfs is also consistent with our data.

Other constraints on the occurrence rate of giant planets around early M dwarfs are available from transit surveys, although they are limited to hot Jupiters and their bins are not as directly comparable to ours due to the observable of these studies being radius, not mass. Gan et al. (2023) found an occurrence rate of 0.27 ± 0.09% for hot Jupiters around early (0.45 ≤ M_* ≤ 0.65) M dwarfs from TESS, with this statistic defined over the ranges 7 R_⊕ ≤ R_P ≤ 2 R_J and 0.8 ≤ P ≤ 10 days. A similar occurrence rate of such planets around our low-mass M dwarfs is fully consistent with our null detection.

Planet formation theories predict an enhancement in giant planet occurrence at the water snow line (Pollack et al. 1996; Ida & Lin 2008; Morbidelli et al. 2015), a phenomenon that has been observed around Sun-like stars; specifically, Fulton et al. (2021) found that giant planet occurrence in the California Legacy Survey of FGKM stars follows a broken power law, with a peak that approximately coincides with the location of the snow line and which represents a factor of 4 increase in occurrence relative to giant planets interior to the snow line. While this sample does contain a small number of M dwarfs, it is dominated by FGK stars, with 83% of the stars in the sample having stellar masses above 0.6 M_⊙ and 98% having stellar masses above 0.3 M_⊙ (see Figure 4 of Rosenthal et al. 2021).

As an aside, we note that Fulton et al. (2021) do not make specific statements about the occurrence of giant planets at the snow line of their M dwarfs. In their Figure 7, they do present occurrence rates as a function of stellar mass for giant plants that are more massive than Saturn and that orbit between 1–5 au, including a bin for 0.3–0.5 M_⊙ M dwarfs. In this bin, they find an occurrence rate of roughly 5%, with the 68% confidence interval ranging from 5% to 12%. However, this 1–5 au bin represents much lower instellations for M dwarfs than for Sun-like stars. Using the scaling relation a_snow = 2.7(M_*/M_⊙)^1.14 au from Childs et al. (2022; adapted from Mulders et al. 2015), the snow line of a 0.4 M_⊙ star occurs at 0.95 au and decreases further to 0.68 au for a 0.3 M_⊙ star. Therefore, while the 1–5 au bin is centered on the snow line for Sun-like stars, it does not represent the occurrence rate of giant planets at the snow line of M dwarfs. Rather, it represents the occurrence rate of giant planets on wider orbits.

To investigate the occurrence rate of giant planets at the snow line of low-mass M dwarfs in our sample, we consider a bin corresponding to zero-albedo equilibrium temperatures of 100 K < T_eq < 150 K, which we subdivide into mass categories of sub-Jupiters (0.3 M_J < M_P sini < 0.8 M_J), Jupiters (0.8 M_J < M_P sin i < 3 M_J), and super-Jupiters (3 M_J < M_P sin i < 10 M_J). Assuming uniform occurrence in log(T_eq) and log(M_P sin i) across each bin, our survey yields 95% confidence upper limits of <1.5% for super-Jupiters, <1.7% for Jupiters, and <4.4% for sub-Jupiters in this region. If we combine the sub-Jupiters and Jupiters into a single 0.3 M_J < M_P sin i < 3 M_J mass bin, we obtain a 95% confidence constraint of <2.4%. While this latter bin averages over a strong gradient in our survey completeness, Fulton et al. (2021) showed that the assumption of uniform occurrence in log(M_P sin i) holds between 0.1 and 3 M_J in their survey of Sun-like stars. This assumption would therefore be well justified if the formation of all giant planets is less efficient around low-mass M dwarfs, but poorly justified if Jupiter analogs are inhibited around low-mass M dwarfs but sub-Jupiter formation remains efficient (which we discuss in Section 5.4).

Our survey indicates that Jupiter analogs—planets with instellations and masses comparable to Jupiter in our solar system—are rare around low-mass M dwarfs. Moreover, if the Fulton et al. (2021) Sun-like star result holds for low-mass M dwarfs (that is, if giant planets are four times more common near the snow line than on interior orbits), our null detection of giant planets in this region implies an even stronger constraint on warm Jupiters than that presented in the previous section. For our Jupiter bin, we find a 68% confidence upper limit of 0.66%. A fourth of this value is 0.17%; of course, it is also possible that the Sun-like occurrence rate distribution does not hold for low-mass M dwarfs, as we discuss in Section 5.4.

How do our findings at and beyond the snow line compare with other works? It is notable that while few giant planets have been detected around low-mass M dwarfs using radial velocities, others have been inferred through microlensing events (e.g., Skowron et al. 2015; Novati et al. 2018; Ryu et al. 2019). However, comparing microlensing and radial velocity statistics is challenging due to the different observables available to each technique. Clanton & Gaudi (2014b) attempted to circumvent these issues by developing a mapping between microlensing and radial velocity parameters that marginalizes over the unknown physical parameters in the microlensing sample. Using these methods, Clanton & Gaudi (2014a) synthesized radial velocity and microlensing surveys to find an occurrence rate of ${2.9}_{-1.5}^{+1.3}$% for M dwarf planets with P = 10⁰–10⁴ days and 1 M_J < M_P sin i < 13 M_J. For periods less than 10² days, their statistics are fully informed by the radial velocity results of Bonfils et al. (2013), with no contribution from the microlensing sample. As discussed in Section 5.1, Bonfils et al. (2013) yield an occurrence rate of 1% (that is, 1/96) for these warm Jupiters; subtracting this from the 2.9% found by Clanton & Gaudi (2014a) implies a roughly 1.9% occurrence rate for P = 10² − 10⁴ days and 1 M_J < M_P sin i < 13 M_J. Assuming uniform occurrence in log(P) and log(M_P sin i) across the bin, our survey obtains N_eff = 148 for giant planets in this regime, producing a 68% confidence upper limit of 0.77% and a 95% confidence upper limit of 2.0% (although note that our sensitivity drops precipitously toward the long-period edge of this bin). Given the large uncertainty in the Clanton & Gaudi (2014a) result, the lack of M_P sin i > 1 M_J planet detections in our study is not necessarily in tension with their occurrence rate. However, those authors note that their value is actually more like a lower limit on the occurrence rate, as they use lower limits in the P = 10²–10³ day range where microlensing sensitivity declines, and they claim their findings are consistent with the nearly twice as large measurement of this occurrence rate from Montet et al. (2014), another radial velocity survey. As these studies focus on more massive M dwarfs (recall that the microlensing lens distribution is centered at M_* = 0.5 M_⊙, while the mean stellar mass in the Bonfils et al. (2013) sample is 0.35 M_⊙), our lack of detections of any giant planets beyond the snow line may be indicative of a decrease in cold Jupiter occurrence around low-mass M dwarfs relative to all M dwarfs.

It is difficult to compare our results directly with Montet et al. (2014), as they report their occurrence rate only in the bin 1 M_J < M_P < 13 M_J out to 20 au. For a 0.2 M_⊙ star, 20 au corresponds to a period of 7 × 10⁴ days; our survey is insensitive to planets at such long periods. We have similar difficulty comparing with the Fulton et al. (2021) occurrence rate of 5% for 100–6000 M_⊕ planets that orbit 0.3–0.5 M_⊙ M dwarfs between 1 and 5 au; we have no sensitivity to a 100 M_⊕ planet at 5 au. While we could compare with these surveys by assuming a specific function for the distribution of planets, as done in Montet et al. (2014), we argue in Section 5.4 that the population of planets around low-mass M dwarfs likely does not follow the same functional form as around Sun-like stars. Furthermore, a comparison of occurrence rate at fixed a may not be particularly informative, as these distances represent substantially different instellations when comparing between early and late M dwarfs.

As stellar activity correlates with rotation rate and hence with radial velocity uncertainty, choosing to restrict our analysis to inactive stars does not greatly affect our sensitivity to Jupiter-sized planets at the snow line: we would not have been able to detect such planets around most active stars due to the larger uncertainties.

However, two giant planets have been previously detected around 0.1–0.3 M_⊙ stars: LHS 252 b, a 0.46 M_J planet with a 204 day period (Morales et al. 2019) and GJ 83.1 b, a 0.21 M_J planet with a 771 day period (Feng et al. 2020; Quirrenbach et al. 2022). While both of these stars are within 15 pc, we measure their Hα emission to be in excess of our −1 Å cutoff and hence do not include them in our sample. In our four observations, both stars show statistically significant radial velocity variations that are consistent with their published ephemerides; however, our observing cadence is too sparse to establish that activity has not contributed to the significance of these signals. It is intriguing that the only known giant planets around low-mass M dwarfs orbit active stars. Most active stars are rotationally broadened, meaning they have large radial velocity uncertainties and it would be difficult to detect planets, even in the absence of activity-induced radial velocity jitter. There are only 26 active, low-mass M dwarfs without close binary companions in our volume-complete sample that have v sin i below 3 km s⁻¹, our detection threshold for line broadening. Two out of these 26 are already known to host a giant planet, and we measure P(χ²) < 1% for three out of the 26 (LHS 252, GJ 83.1, and also GJ 1224). In contrast, we do not detect any giant planets for our 200 inactive M dwarfs.

Is there a reason to expect giant planets to occur preferentially around active M dwarfs? A possible mechanism is the age–activity relation. Younger stars are more active (Kiraga & Stepien 2007; Medina et al. 2022) and have higher metallicities on average. There is compelling evidence that giant planets in shorter period orbits are more common around metal-rich stars than metal-poor stars (Santos et al. 2001; Fischer & Valenti 2005; Fulton et al. 2021; Rosenthal et al. 2022), although the picture is less clear for M dwarfs (Gaidos & Mann 2014). An inactive sample is therefore biased toward older, metal-poor, possibly planet-deficient stars. However, neither GJ 83.1 nor LHS 252 is particularly metal rich, with [Fe/H] estimates of −0.13 and 0.02 dex, respectively (Newton et al. 2015).

The planetary occurrence rate is often parameterized as a function of planetary mass and period (e.g., Cumming et al. 2008). As discussed above, Fulton et al. (2021) advocated for a broken power-law parameterization for the giant planets of Sun-like stars. One might hypothesize that the giant planet occurrence function has the same shape around M dwarfs as it does around Sun-like stars, with the entire distribution scaled down by some constant. However, our results suggest that this is not the case. If the activity of the known planet hosts from the previous section is a coincidence, the discrepancy in occurrence rates could instead be explained by the relatively low masses of LHS 252 b and GJ 83.1 b. LHS 252 b would fall in our sub-Jupiter bin, where we place a weaker 95% confidence upper limit of 4.4%. For 0.21 M_J planets like GJ 83.1 b, our sensitivity drops to only a handful of stars in our snow-line installation bin (N_eff = 4). It may be that while Jupiter analogs are rare around low-mass M dwarfs, lower-mass giant planets still form, in a way that deviates from the uniform distribution in log(M_P sin i) found for the giant planets of Sun-like stars. Future surveys with higher sensitivity to 0.1–0.3 M_J planets beyond the snow line will be necessary to explore this hypothesis.

However, lower-mass giant planets are unlikely to have a Jupiter-like impact on inner terrestrial planets, as Jupiter's migration profoundly shaped our inner solar system (Walsh et al. 2011). For Sun-like stars, a Jupiter-like mass is required to open a gap in the protoplanetary disk, or a Saturn-like mass in regions of low turbulent activity; smaller giant planets only partially open a gap in their disk (Baruteau & Masset 2013), meaning they operate in a different regime of migration and lack Jupiter's capacity to suppress the inward flow of pebbles. While one might intuitively expect this gap-opening mass to decrease for lower-mass stars, the simulations of Sinclair et al. (2020a) indicate that the threshold mass either increases or holds constant with decreasing stellar mass; those authors note that this behavior is the result of disks around lower-mass stars being geometrically thicker due to reduced gravity, yielding increased pressure forces. Lower-mass giant planets around M dwarfs will therefore not open Jupiter-like gaps.

Our results coupled with microlensing studies also hint at differing behavior with instellation. While the occurrence of giant planets peaks just beyond the snow line for Sun-like stars (Fulton et al. 2021), we find no giant planets at these instellations. Meanwhile, Poleski et al. (2021) analyzed 20 yr of microlensing data to conclude that, on average, every microlensing star hosts a wide-orbit giant planet (5 < a < 15 au and 10⁻⁴ < M_p/M_* < 0.033); note that the microlensing lens distribution is dominated by M dwarfs (Gould et al. 2010). There have also been reported detections of >1 M_J planets around low-mass M dwarfs at very wide orbital separations from direct imaging (Gaidos et al. 2022) and disk studies (Curone et al. 2022). These results hint that giant planet occurrence may peak at lower instellations around low-mass M dwarfs as compared to Sun-like stars. Such a distribution could indicate a pathway for giant planet formation governed by processes unrelated to the water snow line, such as disk instability (Boss 2006; Mercer & Stamatellos 2020). Other possible mechanisms to move giant planets to wide orbits include scattering by planet–planet interactions (e.g., Rasio & Ford 1996) or outward migration in resonance (e.g., Crida et al.2009).

For Sun-like stars, Rosenthal et al. (2022) found that 41% of stars with a close-in, small planet also hosted an outer giant, in contrast to a 17.6% occurrence of outer giants overall. We are therefore particularly attentive to the sensitivity curves for the ten stars in our sample with published small planets: GJ 1214 (Charbonneau et al. 2009), GJ 1132 (Berta-Thompson et al. 2015; Bonfils et al. 2018), LHS 1140 (Dittmann et al. 2017; Ment et al. 2019), GJ 1265 & LHS 350 (Luque et al. 2018), LHS 3844 (Vanderspek et al. 2019), LTT 1445A (Winters et al.2019, 2022; Lavie et al. 2023), GJ 1061 (Dreizler et al. 2020), GJ 1057 (Bauer et al. 2020), and GJ 585 (Harakawa et al. 2022). While we do not have the radial velocity sensitivity to recover these small planets, we can provide constraints on the existence of further out giant planets to inform our understanding of exoplanetary system architectures.

In Figure 5, we show our sensitivity to giant planets around each of these small-planet hosts. For all stars, we are likely to detect a planet that is Jupiter-like in mass and instellation, with our sensitivity varying from just under 50% for LHS 350 to nearly 100% for GJ 1061. As we only have four observations of each star, the sensitivity curves are bumpier than the smooth results of Figure 3; for a given star, there are specific periods in which a planet could feasibly evade detection, although these hiding spots are averaged out when we consider the sample as a whole. For the purposes of Figure 5, we have rerun our analysis with 1000 samples in mass and period to allow the nonhomogeneity of the sensitivity curves to be more readily appreciated. While we cannot definitively rule out the existence of a giant planet beyond the snow line for any individual system, Jupiter-analog companions to inner terrestrial planets—i.e., planetary system architectures like our own—cannot be commonplace for low-mass M dwarfs. If 41% of these stars had a planet that was Jupiter-like in mass and instellation (a percentage motivated by the outer giant companion occurrence rate around Sun-like stars from Rosenthal et al. 2022), there is a 98% chance we would have detected this planet around at least one of the ten stars. An abundance of lower-mass giant companions beyond the snow line (e.g., M_P sin i ≤ 0.3 M_J) is not ruled out by our observations: if all ten stars had a 0.3 M_J planet at the Jovian instellation, there is a 78% chance that we would have detected at least one, but this chance drops to 43% when we lower the incidence to 41%.

Zoom In Zoom Out Reset image size

What does a lack of outer Jovians imply for the inner terrestrial planets in these systems (notably, the only terrestrial planets amenable to atmospheric study for at least the next decade)? While simulations on this topic are beyond the scope of this work, one can look to the large body of literature discussing the impact of Jupiter on the evolution of our own solar system. Models such as Walsh et al. (2011) suggest that the migration of Jupiter and Saturn led to the truncation of the planetesimal disk at 1 au, ultimately limiting the size of the terrestrial planets. Meanwhile, Izidoro et al. (2015) suggest that Jupiter acted as a dynamical barrier to the inward migration of gas giant cores, preventing them from migrating into the inner solar system from beyond the snow line and explaining why our solar system lacks a super-Earth (although note Bryan et al.2019, who found that the presence of super-Earths correlates positively with the presence of an outer Jovian). Similarly, Bitsch et al. (2021) argue that a Jovian beyond the snow line can block the inward flow of water-rich pebbles, resulting in a drier inner stellar system. Systems without an outer Jovian may therefore have wetter (and potentially, larger) small planets; notably, Luque & Pallé (2022) recently identified a population of water worlds among the small planets that transit M dwarfs (although note Ment & Charbonneau 2023, who found that volatile-rich small planets are less common around mid-to-late M dwarfs). Giant planets are thought to also set the terrestrial water budget in our own solar system, with Raymond & Izidoro (2017) arguing that Jupiter and Saturn were responsible for delivering Earth's small amount of water through planetesimal scattering.

Childs et al. (2022) note another implication of a lack of outer Jovians: without such a planet, a stellar system is unlikely to have an asteroid belt or a mechanism of delivering asteroids to inner terrestrial planets, and these asteroid impacts may be necessary for the origin of life (e.g., Osinski et al. 2020). Using modeling from Childs et al. (2019), they further argue that lower-mass giant planets (specifically, they replaced Jupiter with a planet of mass 0.14 M_J) are unable to sustain suitable asteroid bombardment; more massive giant planets are needed. In the case of Earth, large asteroid impacts are thought to be responsible for the production of a reducing atmosphere (Sinclair et al. 2020b), which ultimately led to the emergence of prebiotic chemistry (Benner et al. 2019).