# No Sub-Saturn-mass Planet Desert in the CORALIE/HARPS Radial-velocity Sample

David P. Bennett, Clément Ranc, Rachel B. Fernandes

► To cite this version:

David P. Bennett, Clément Ranc, Rachel B. Fernandes. No Sub-Saturn-mass Planet Desert in the CORALIE/HARPS Radial-velocity Sample. The Astronomical Journal, 2021, 162, 10.3847/1538-3881/ac2a2b. insu-03748182

HAL Id: insu-03748182

https://insu.hal.science/insu-03748182v1

Submitted on 2 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# No Sub-Saturn-mass Planet Desert in the CORALIE/HARPS Radial-velocity Sample

David P. Bennett¹,², Clément Ranc³, and Rachel B. Fernandes⁴,⁵

¹ Laboratory for Exoplanets and Stellar Astrophysics, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA; bennettd@umd.edu

² Department of Astronomy, University of Maryland, College Park, MD 20742, USA

³ Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, F-75014 Paris, France

⁴ Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721, USA

⁵ Earths in Other Solar Systems Team, NASA Nexus for Exoplanet System Science, USA

Received 2021 April 7; revised 2021 September 19; accepted 2021 September 24; published 2021 November 16

## Abstract

We analyze the CORALIE/HARPS sample of exoplanets found by the Doppler radial-velocity method for signs of the predicted gap or “desert” at 10–100 $M_⊕$ caused by runaway gas accretion at semimajor axes of <3 au. We find that these data are not consistent with this prediction. This result is similar to the finding by the MOA gravitational microlensing survey that found no desert in the exoplanet distribution for exoplanets in slightly longer period orbits and somewhat lower host masses (Suzuki et al. 2018). Together, these results imply that the runaway gas accretion scenario of the core accretion theory does not have a large influence on the final mass and semimajor axis distribution of exoplanets.

*Unified Astronomy Thesaurus concepts: Exoplanets (498); Exoplanet astronomy (486); Exoplanet formation (492)*

## 1. Introduction

The runaway gas accretion process has long been considered to be a major step in the formation of gas-giant planets in the core accretion theory (Pollack et al. 1996). Several authors (Ida & Lin 2004a; Mordasini et al. 2009a, 2009b) have predicted that this runaway gas accretion process should produce a gap or “desert” in the mass distribution of planets. The abstract of Ida & Lin (2004a) states “Since planets” masses grow rapidly from 10 to 100 $M_⊕$, the gas-giant planets rarely form with asymptotic masses in this intermediate range”. They specifically indicate that this desert should exist at orbital major axes of 0.2–3 au. The Bern population synthesis group (Mordasini et al. 2009a, 2009b) predict a lower amplitude desert, a factor of 2–3 at masses between 30 and 100 $M_⊕$. The paper presenting the data that we analyze in this paper (Mayor et al. 2011) claims that these data tend to confirm this prediction. They state: “After correction of detection biases, we see even more clearly the importance of the population of low-mass planets on tight orbits, with a sharp decrease of the distribution between a few Earth masses and ~40 $M_⊕$. We note that the planet population synthesis models by Mordasini et al. (2009a, 2009b) predicted such a minimum in the mass distribution at precisely this mass range.” As we shall see, our interpretation differs somewhat from this, and it strongly contradicts a more recent version of the Bern group’s population synthesis code (Emsenhuber et al. 2020a, 2020a), which predicts a deeper desert, with a planet-occurrence rate drop of a factor of ~10.

In this paper, we present a reanalysis of the Mayor et al. (2011) Doppler radial-velocity-detected planet sample, that was recently analyzed by Fernandes et al. (2019). While Fernandes et al. (2019) focused primarily on the orbital period dependence of the planet distribution, our analysis is focused primarily on the dependence on planetary mass, so that we can investigate the proposed sub-Saturn-mass planetary desert. Before considering the Mayor et al. (2011) sample, however, we review recent results from gravitational microlensing that seem to contradict this cold sub-Saturn- prediction in Section 2. Then, in Section 3, we present the binned, completeness corrected planet distribution with an estimate of the error bars for each bin, and we show

that this data shows no significant sign of a sub-Saturn-mass planet desert. In Section 4, we present a more powerful forward modeling method to search for evidence of a planet desert in the Mayor et al. (2011) data using a planet-occurrence model that employs a Gaussian desert feature, and again, we find no evidence of a sub-Saturn-mass planet desert. We also compare our results to the Fernandes et al. (2019) analysis and recover their main result. We then consider the latest generation of planet population synthesis models from the Bern group (Emsenhuber et al. 2020a, 2020b) in Section 5, and we find that our forward Gaussian desert model does a reasonable job at characterizing the planet desert feature in their simulated planet population that is contradicted by the Mayor et al. (2011) CORALIE/HARPS data. Emsenhuber et al. (2020c) have also pointed out that deep planet desert at $20 M_⊕ < M < 200 M_⊕$ seen in their models is contradicted by the Mayor et al. (2011) data. Next, we show that a possible sub-Jupiter “radius desert” that may be seen in the Kepler data for wider orbit planets is likely to be due to the planet mass-radius relation instead of a desert in the mass distribution. In Section 6, we present our conclusions and discuss our expectations for future developments in the study of wide-orbit planet demographics.

Note that the planet desert that we investigate in this paper differs from other “deserts” that have been observed in or suggested by Kepler data. These are considered in the Appendix.

## 2. Search for a Planet Desert Beyond the Snow Line

A recent comparison of theoretical population synthesis calculations, based on the core accretion theory, with the distribution of planets detected by gravitational microlensing (Suzuki et al. 2018), orbiting beyond the snow line, found that the theoretical models predicted a strong desert at mass ratios of $10^{-4} \le q \le 4 \times 10^{-4}$ that was not reproduced by the microlensing results. Instead, the exoplanet mass-ratio distribution found by microlensing indicated a smooth, power-law distribution down to a peak at a mass ratio of $q \sim 10^{-4}$, with fewer planets at lower mass ratios (Udalski et al. 2018; Jung et al. 2019). This discrepancy between theory and observation was thought to be likely due to deficiencies in the assumptions
and approximations in the population synthesis calculations instead of problems with the fundamental concepts of the core accretion theory. One possibility that has been suggested for our own solar system is that captured planetesimals could heat the accreting gaseous envelope, which would act to slow accretion, and other possibilities are that the scale height or viscosity of the protoplanetary disk may be smaller than the values adopted for these population synthesis calculations (Dobbs-Dixon et al. 2007; Fung et al. 2014; Garaud & Lin 2007). A reduction of the assumed disk viscosity in the Bern group’s population synthesis model did indicate that this sub-Saturn-mass desert could disappear with lower viscosity (Suzuki et al. 2018), as earlier work had indicated (Szulágyi et al. 2014). Annular structures in protoplanetary disks observed by ALMA are thought by many to be due to planets in the process of formation. A comparison of these structures in a number of protoplanetary disks also indicates that the growth of gas-giant planets is slower than predicted by standard runaway gas accretion models (Nayakshin et al. 2019), so runaway gas accretion is disfavored unless these annular structures are caused by another mechanism.

### 3. Planet Occurrence for the CORALIE/HARPS Sample

The CORALIE/HARPS data set (Mayor et al. 2011) that we analyze in this paper was recently analyzed by Fernandes et al. (2019), who considered the distribution of planets in the 30–6000 $M_⊕$ mass range. They fit this with a broken power-law model that has a break and peak at a period of 1720 ± 430 days. When they add a power law for the mass distribution of planets in this mass range, the range of possible locations of the peak expands to 2100 ± 1200 days.

The Fernandes et al. (2019) analysis of the Mayor et al. (2011) data set included a tabulation of the survey completeness from Figure 6 of Mayor et al. (2011), and we use this tabulation of the survey completeness for our analysis. As a first step in our analysis, we present the binned $M \sin i$ distribution of the Mayor et al. (2011) sample in Figure 1, where M is the planet mass and i is the inclination of the planet’s orbital plane to the line of sight. We break up the Mayor et al. (2011) sample into four different period ranges, each spanning a factor of eight in period: 3.75–30 days, 30–240 days, 240–1920 days, and 1920–15,260 days. The boundary between the third and fourth period ranges is set to be 1920 days because this corresponds to the approximate outer semimajor axis boundary of 3 au of the Ida & Lin (2004a) prediction for a $1 M_⊙$ host star and is within 0.5σ of the break in the giant planet-occurrence rate found by Fernandes et al. (2019).

In order to relate the planets found in the CORALIE/HARPS data set to the intrinsic planetary occurrence rate, we correct for the survey completeness, which gives the colored histograms with error bars shown in Figure 1. The occurrence rate is given by

$$ \eta = \frac{1}{N_*} \sum_{j}^{n_p} \frac{1}{C(M_j \sin i_j, P_j)}, \quad (1) $$

where $N_* = 822$ is the number of stars in the CORALIE/HARPS sample, $n_p$ is the number of detected planets, and $M_j$, $i_j$, and $P_j$ are the mass, orbital inclination, and the period of the jth planet in the sample. The survey completeness for mass M, inclination i and period P is given by $C(M \sin i, P)$.

The gray histograms without error bars in Figure 1 indicate the raw, uncorrected histograms. The gray and colored histograms merge for large $M \sin i$ and small periods, as expected, because these planets can be detected with nearly 100% completeness.

In order to present the planet-occurrence rate as a function of mass (or $M \sin i$) and period, we bin the data in both $M \sin i$ and period. Because the CORALIE/HARPS data set is sensitive to a very wide range of planetary orbital periods and masses, it is necessary to select bins that span a relatively large range in both $M \sin i$ and period in order to have a reasonable number of planets in each bin. This complicates the issue of determining the error bars for the completeness corrected occurrence rate for each bin, because the completeness can vary significantly across some of the bins. Using the average completeness over the range of parameters for each bin could induce an error if the planet-occurrence rate also varies significantly over the bin size. Therefore, we use the completeness for the $M \sin i$ and P for each planet in Equation (1). This means that planets in the same bin will have different contributions to the planet-occurrence rate, $\eta$, so that we cannot use Poisson statistics to calculate the error bars for $\eta$. Therefore, we employ a bootstrap method to determine the error bars. We take the number of planets in each bin as the expected number for a Poisson distribution, and then we randomly select the planets from the list of planets in the bin. We repeat this procedure 100,000 times for each bin to determine the root-mean-square dispersion of the completeness corrected $\eta$ value, which we use as the error estimate for the $\eta$ value for each bin. For bins without any planets, we calculate the mean completeness for the bin, and then use the one-sided $1\sigma$ upper limit from Poisson statistics for zero detections to determine the $1\sigma$ upper limit.

The gray shaded regions in the upper three panels of Figure 1 indicate the mass range of the planet desert predicted by Ida & Lin (2004a). As noted by Mayor et al. (2011), there is a significant decrease in the exoplanet-occurrence rate, $\eta$, between 10 and ~40 $M_⊕$ at the shorter orbital periods, particularly in the 3.75 < P < 30 days range. (With our binning, this minimum appears at ~30 $M_⊕$.) So, planets at the bottom edge of the predicted “desert” are clearly more common, but to claim a desert, the exoplanet occurrence must increase at higher masses, and such a feature is not seen in Figure 1. Also, as mentioned in the introduction, Mayor et al. (2011) point to a minimum in the planet distribution at ~40 $M_⊕$ as being consistent with the Mordasini et al. (2009a, 2009b) prediction of a “moderate” planet desert in the range 30–100 $M_⊕$, with a reduction of a factor of 2–3 from the number of Jupiter-mass planets. Our analysis does indicate a minimum at ~30 $M_⊕$ in the 3.75 < P < 30 days bin, but our error bar estimates indicate that this is not significant even at the $1\sigma$ level. It is possible that this minimum might appear more significant with a different binning, but we will address this issue in Section 4 with a forward modeling approach that involves no binning. Also, note that there is no evidence for the predicted factor of 2–3 deficit of planets at 30–100 $M_⊕$ compared to Jupiter-mass (~318 $M_⊕$) planets (Mordasini et al. 2009a, 2009b) in any of the period ranges shown in Figure 1. There does seem to be a marginal increase of Jupiter-mass planets compared to 100–200 $M_⊕$ planets in the 240–1920 days period range, but this is only marginally significant and is not duplicated at shorter or longer periods.

In this discussion, we have avoided the distinction between planet masses, M, and the $M \sin i$ values that are measured by radial-velocity monitoring. Because small inclinations, $i \ll 1$, are disfavored geometrically, this is generally expected to have
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Planet-occurrence rate histograms for the planets in the Mayor et al. (2011) sample in four different orbital period ranges, each spanning a factor of eight in period. The histograms without any completeness corrections are shown in gray and the light gray shaded regions show the location of the “planet desert” predicted by Ida & Lin (2004a).</figcaption></figure>

a small effect, but it could conceivably be important if the occurrence rate, $\eta$, is strongly dependent on planet mass. Fortunately, this complication disappears with the forward modeling approach presented in the next section.

# 4. Forward Modeling Analysis

We now turn to a more rigorous forward modeling analysis to search for evidence of this predicted desert, attributed to the runaway gas accretion process of the core accretion theory. Forward modeling allows us to investigate the possibility of a desert with somewhat different parameters than proposed by Ida & Lin (2004a) and Mordasini et al. (2009a, 2009b), and it enables us to avoid the uncertainties relating to the binning of data and error bar estimates mentioned in the previous section.

## 4.1. Exoplanet Desert Occurrence Rate Models

Following earlier Doppler radial-velocity analyses (Cumming et al. 2008; Fernandes et al. 2019; Wittenmyer et al. 2020), we base our model on power laws in both mass, *M*, and period, *P*. The previous analysis of the CORALIE/HARPS sample (Fernandes et al. 2019) found a break in the planet-occurrence rate at a period of *P* = 2100 ± 1200 days, which is consistent with the *P* = 1920 days boundary between the third and fourth period ranges presented in Figure 1. (Fernandes et al. (2019) indicate that this broken power law applies to masses in the range 30–6000 *M*⊕, but Figure 1 makes it clear that the Mayor et al. (2011) sample has little sensitivity for *M* < 100 *M*⊕ at periods longer than the break.) Since the *M* sin i dependence does not appear to follow the same power law in these four period ranges, we elect to fit these four period ranges separately without any power-law break.
In order to describe a planetary desert, we add the following desert reduction factor, $D(P, M_{des})$ to the power law model:

$$D(M, M_{des}) = 1 - B \exp\left(-\frac{[\log M - \log M_{des}]^2}{2\sigma_{des}^2}\right), \quad (2)$$

where the planet desert is centered at period $M_{des}$ with a Gaussian dispersion of $\sigma_{des}$ in $\log P$. (Throughout this paper, log refers to the base-10 logarithm.) This desert reduction factor multiplies the power-law factors to yield our planet desert mass-function model:

$$\begin{aligned} f_{des}(M, P; A, B, \log M_{des}, \sigma_{des}, p, m) \\ &\equiv \frac{d^2 N_{pl}}{d \log M d \log P} \\ &= A \left(\frac{P}{1995 \text{ days}}\right)^p D(M, M_{des}) \left(\frac{M}{100 M_{\oplus}}\right)^m. \end{aligned} \quad (3)$$

We select the Gaussian function to describe the desert because we think that a smooth function is more plausible than a step function, and we do find that it provides a reasonably good description of the planet distribution predicted by the most recent version of the Bern group's planetary population synthesis model (Emsenhuber et al. 2020a, 2020b) in some period ranges, as discussed in Section 5. Note that we have selected the pivot point for the period power-law to be $10^{3.3} = 1995$ days in order to match the break in the power law when we include all periods longer than 3.75 days, as discussed in Section 4.2.

One complication for our forward modeling effort is that the completeness, $C(M \sin i, P)$, depends on the inclination angle, $i$, while intrinsic planet properties will depend only on the planet mass, $M$. Fortunately, it is well established that there is no physical process that is likely to produce a nonrandom distribution of planetary orbital planes as seen from Earth, so we can treat the inclination angles, $i$, as a random variable to integrate over. We introduce the variable $\mathfrak{M} \equiv M \sin i$ to refer to the observable parameter measured by the Doppler radial-velocity method.

We use a standard forward modeling Bayesian likelihood analysis (Alcock et al. 1996; Suzuki et al. 2016) to determine the posterior distribution of models that are consistent with the data. The likelihood function is given by

$$\begin{aligned} \mathcal{L}_{des}(A, B, \log M_{des}, \sigma_{des}, p, m) \\ = e^{-N_{exp}} \times \prod_{j}^{N_{s}} \int di \sin i f_{des}(\mathfrak{M}_j / \sin i, P_j; A, B, \\ \log M_{des}, \sigma_{des}, p, m) \times C(\mathfrak{M}_j, P_j), \end{aligned} \quad (4)$$

where

$$\begin{aligned} N_{exp} \\ = \iiint dM dP di \sin i f_{des}(M, P; A, B, \log M_{des}, \sigma_{des}, p, m) \\ \times C(M \sin i, P), \end{aligned} \quad (5)$$

is the expected number planets detected for the model with the given parameters: A, B, $\log M_{des}$, $\sigma_{des}$, p, and m. We evaluate this likelihood function with a Markov chain Monte Carlo (MCMC) for each of the period ranges: 3.75–30 days,

30–240 days, 240–1920 days, and 1920–15,260 days, with chains of 140,000–330,000 steps after 2000 burn-in steps. The results of these calculations are summarized in Figure 2. The best-fit models are shown as the black curves, while the dark and light gray shading indicates the central 68.3% and 95.4% likelihood distributions for the planet-occurrence rate as a function of mass from the MCMC calculations. These plots are shown for the median $\log P$ values for each of the period ranges. These central periods are 10.6 days for $3.75 < P < 30$ days, 84.9 days for $30 < P < 240$ days, and 679 days for $240 < P < 1920$ days.

The histograms of the exoplanet occurrence from the CORALIE/HARPS data set from the top three panels of Figure 1 are also shown in Figure 2. However, the x axis for Figure 2 is the planet mass, $M$, instead of the Doppler radial-velocity observable $\mathfrak{M} = M \sin i$. So, we shift the x-axis from $\mathfrak{M}$ to $M \approx \mathfrak{M} / \sin i_{med} = 1.1547$, where $i_{med} = \cos^{-1}(0.5)$ is the median inclination value, assuming a random orientation of the orbital planes. This has no effect on the Bayesian likelihood analysis, since we integrate over all inclination angles, $i$, which uses the measured period and $\mathfrak{M}$ values without any binning. The parameters of the best-fit and average desert power-law models are given in Table 1. The scatter in the normalization parameter A is large because the pivot point of the period power law is set to a fairly large period of 1995 days, and this results in a relatively large scatter in the normalization parameter, A. Therefore, we also include a parameter, $A_{Pm} = A(P_m/1995 \text{ days})^p$, where $P_m$ logarithmic median of the period range being considered. For the different period ranges, we have $P_m = 10.607, 84.853, \text{ and } 678.820$ days for the $3.75–30, 30–240, \text{ and } 240–1920$ days ranges, respectively. These $A_{Pm}$ values reflect the variation in the normalization in the middle of the period ranges under consideration.

Figure 2 indicates that the best-fit models do have narrow exoplanet deserts, and the deserts in the 3.75–30 days and 30–240 days period ranges are centered in the range of 10–100 days predicted by Ida & Lin (2004a). The depth of the desert features in the best-fit models, represented by the B, is fairly deep at their centers, reducing the occurrence rate by over 95% for the 3.75–30 day and 240–1920 day ranges, and by 76.7% in the 30–240 day range. However, the distribution of models consistent with the data, as represented by the gray shaded regions in Figure 2 tell a different story. The central 68.3% and 95.4% model occurrence rate distributions generally do not follow the best-fit models, although the lowest extent of the central 95.4% occurrence rate distribution for the 3.75–30 days period range does have a dip at the location of the best-fit desert model, but there is virtually no such feature in the central 68.3% occurrence rate distribution. These features are also narrower than the Ida & Lin (2004a) prediction of a desert spanning a factor of 10 in mass. The Gaussian $\sigma_{des}$ values (in $\log M$) range from $\sigma_{des} = 0.110$ for the 240–1920 days period range to $\sigma_{des} = 0.155$ for the 3.75–30 days period range. These imply full-width half-max values in the range 1.8–2.3 in $M$. The most likely explanation for this is that the very narrow desert features of the best-fit models are due to statistical noise in the observed exoplanet distribution. This would certainly explain why the central 68.3% range of the occurrence rate distribution does not show these features.

We do note that the location of the best-fit desert in the 3.75–30 day and 30–240 day period ranges do roughly correspond to the minimum at $M \sin i = 40 M_{\oplus}$ mentioned
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The desert power-law planet-occurrence rate models compared to the histograms from the top three period range panels of Figure 1. The black curves are the best-fit models, and the dark and light gray shaded regions indicate the central 68.3% and 95.4% distribution of planet-occurrence rates as a function of mass from the MCMC distributions. These indicate no significant evidence of a sub-Saturn-mass planet desert in the Mayor et al. (2011) CORALIE/HARPS data set.</figcaption></figure>

<table><caption>**Table 1**
Exoplanet Desert Power-law Model Parameters</caption><thead><tr><th rowspan="3">Parameter</th><th colspan="6">Period ranges</th></tr><tr><th colspan="2">3.75–30 days</th><th colspan="2">30–240 days</th><th colspan="2">240–1920 days</th></tr><tr><th>Best</th><th>Average</th><th>Best</th><th>Average</th><th>Best</th><th>Average</th></tr></thead><tbody><tr><td>A</td><td>0.745</td><td>4.03 ± 2.44</td><td>0.104</td><td>0.38 ± 0.32</td><td>0.105</td><td>0.178 ± 0.087</td></tr><tr><td>A<sub>Pm</sub></td><td>0.0177</td><td>0.020 ± 0.009</td><td>0.0561</td><td>0.068 ± 0.022</td><td>0.0920</td><td>0.125 ± 0.052</td></tr><tr><td>B</td><td>0.977</td><td>0.52 ± 0.26</td><td>0.767</td><td>0.48 ± 0.24</td><td>0.954</td><td>0.51 ± 0.27</td></tr><tr><td>log M<sub>des</sub></td><td>1.529</td><td>2.08 ± 0.48</td><td>1.698</td><td>2.08 ± 0.42</td><td>2.183</td><td>1.61 ± 0.44</td></tr><tr><td>σ<sub>des</sub></td><td>0.155</td><td>0.45 ± 0.28</td><td>0.137</td><td>0.55 ± 0.28</td><td>0.110</td><td>0.53 ± 0.29</td></tr><tr><td>p</td><td>0.714</td><td>0.97 ± 0.18</td><td>0.194</td><td>0.473 ± 0.23</td><td>0.123</td><td>0.30 ± 0.24</td></tr><tr><td>m</td><td>-1.229</td><td>-1.18 ± 0.15</td><td>-1.008</td><td>-0.94 ± 0.12</td><td>-0.407</td><td>-0.49 ± 0.14</td></tr></tbody></table>

by Mayor et al. (2011), particularly in their sample with periods <100 days. However, our analysis indicates that this minimum is not likely to be significant. In fact, this conclusion seems

consistent with the completeness corrected histogram shown in Figure 12 of Mayor et al. (2011), where the bin centered near $M \sin i = 40 M_⊕$ is consistent at $1\sigma$ with the next two larger
<table><caption>**Table 2**
Pure Power-law Model Parameters</caption><thead><tr><th rowspan="3">Parameter</th><th colspan="4">Period Ranges</th></tr><tr><th colspan="2">3.75–30 days</th><th colspan="2">30–240 days</th><th colspan="2">240–1920 days</th></tr><tr><th>Best</th><th>Average</th><th>Best</th><th>Average</th><th>Best</th><th>Average</th></tr></thead><tbody><tr><td>A</td><td>0.542</td><td>1.77 ± 0.62</td><td>0.081</td><td>0.22 ± 0.16</td><td>0.099</td><td>0.113 ± 0.035</td></tr><tr><td>A<sub>Pm</sub></td><td>0.0148</td><td>0.0148 ± 0.0044</td><td>0.0481</td><td>0.048 ± 0.009</td><td>0.0812</td><td>0.083 ± 0.016</td></tr><tr><td>p</td><td>0.701</td><td>0.90 ± 0.12</td><td>0.165</td><td>0.40 ± 0.23</td><td>0.179</td><td>0.26 ± 0.24</td></tr><tr><td>m</td><td>-1.251</td><td>-1.26 ± 0.13</td><td>-0.972</td><td>-0.98 ± 0.11</td><td>-0.367</td><td>-0.38 ± 0.09</td></tr></tbody></table>

bins. The large rise at lower masses is significant, as is the case with our histogram in the top panel of Figure 2, although our error bars are more conservative. However, this is not support the runaway accretion prediction that planets should grow rapidly across the 10–100 $M_⊕$ mass interval.

## 4.2. Power-law Occurrence Rate Models

Since we have concluded that there is no significant indication of a planet desert in the CORALIE/HARPS data set, it is sensible to consider simpler models that do not include the desert. The model parameters for such models are presented in Table 2, which includes both the best-fit models and the MCMC averages. These power-law models are given by Equation (3) with $B \equiv 0$ and therefore $D(M, M_{des}) \equiv 1$. As was the case for the planet desert models, the parameter $A_{Pm}$ is not a model parameter. Instead, it reflects the normalization of the occurrence rate function at the logarithmic center of the period range under consideration. So, the scatter in $A_{Pm}$ is a better indication of the uncertainty in the occurrence rate at the center of the period range under consideration than A.

We can compare these power-law models to the corresponding power-law plus desert models to see if the addition of the desert feature significantly improves the fit. We can do this by comparing the log $\mathcal{L}$ values for the models without and with the desert features. We find that the desert features improve the log $\mathcal{L}$ values by 1.526, 0.722, and 0.944 for the 3.75–30, 30–240, and 240–1920 days period ranges, respectively. Using frequentist statistics, these would imply that the pure power-law models were disfavored by $p = 0.030, 0.189, \text{ and } 0.115$, for these three period ranges. Following Equation (27) of Trotta (2008), we can compute upper bounds on the Bayes factor of 3.52, 1.17, and 1.48 for the 3.75–30, 30–240, and 240–1920 days period ranges. These correspond to weak evidence favoring the power-law plus desert model for the 3.75–30 days period range, and no significant evidence in favor of the desert additions to the power-law models at the longer period ranges. While there may be marginal evidence that the desert model is an improvement over the pure power-law model for the 3.75–30 days period range, we do not consider this to be evidence in favor of the desert feature itself based upon the MCMC results.

Finally, to compare with Fernandes et al. (2019), we also consider a model with a broken power law in the period with masses constrained to the range $0.1 M_{\text{Jup}} \le M \sin i \le 20 M_{\text{Jup}}$. This mass range is considered (Fernandes et al. 2019) to be the range of giant planets. Our broken power-law mass-function

model is given by

$$ 
\begin{aligned}
f_{\text{bpl}}(M, P; A, P_{\text{br}}, p_1, p_2, m) &\equiv \frac{d^2 N_{\text{pl}}}{d \log M \, d \log P} \\
&= A \left[ \left( \frac{P}{P_{\text{br}}} \right)^{p_1} \Theta(P - P_{\text{br}}) + \left( \frac{P}{P_{\text{br}}} \right)^{p_2} \Theta(P_{\text{br}} - P) \right] M^m,
\end{aligned} 
\quad (6) $$

where $\Theta$ is the Heaviside theta function (or step function), $P_{\text{br}}$ is the period break and $p_1, p_2$, and $m$ are the power-law exponents for small periods, large periods, and the planet mass, $M$, respectively. The likelihood function then becomes

$$ 
\begin{aligned}
\mathcal{L}_{\text{bpl}}(A, P_{\text{br}}, p_1, p_2, m) = e^{-N_{\text{exp}}} \\
& \times \prod_{j}^{N_r} \int di \, \sin i \, f_{\text{bpl}}(\mathfrak{M}_j / \sin i, P_j; A, P_{\text{br}}, p_1, p_2, m) \\
& \times C(\mathfrak{M}_j, P_j),
\end{aligned} 
\quad (7) $$

and the number of expected planet detections becomes

$$ N_{\text{exp}} = \iiint dM \, dP \, di \, \sin i \, f_{\text{bpl}}(M, P; A, P_{\text{br}}, p_1, p_2, m) \times C(M \sin i, P). \quad (8) $$

Evaluating this with a MCMC, as we did with the planet desert models, we find the best-fit model, as all the range of models from the MCMC distribution. The model parameter results are shown in the first three columns of Table 3 indicating a power-law break at $P_{\text{br}} = 10^{3.20 \pm 0.17} = 1600_{-500}^{+800}$ days. Note that we have used $\log(P_{\text{br}})$ instead of $P_{\text{br}}$ as the parameter to indicate the location of the break in the period distribution in our MCMC calculations. This is equivalent to assuming that the prior distribution of the $P_{\text{br}}$ variable in uniform in $\log(P_{\text{br}})$ rather than uniform in $P_{\text{br}}$. These results are consistent with the Fernandes et al. (2019) results for the same model as determined by the EPOS program, which are reproduced in the fourth column of Table 3. (Note that the normalization parameter, $A$, has been modified from the value in the Fernandes et al. (2019) paper to account for the fact that they define the broken power law mass function over $d \ln M \, d \ln P$, whereas we use $d \log_{10} M \, d \log_{10} P$.) Thus, we confirm the main conclusion of Fernandes et al. (2019) that the giant planet-occurrence rate has a maximum in the vicinity of the snow line.

## 5. Comparison of CORALIE/HARPS Data and Gaussian Desert Model to Generation III Bern Models

In this section, we compare the CORALIE/HARPS exoplanet sample (Mayor et al. 2011) to the results of Generation
<table><caption>**Table 3**
Broken Power-law Model Parameters for $0.1 M_{\text{Jup}} \le M \sin i \le 20 M_{\text{Jup}}$</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="2">P range: 1–11000days</th><th rowspan="2">EPOS</th></tr><tr><th>Best</th><th>Average</th></tr></thead><tbody><tr><td>A</td><td>0.167</td><td>0.165 ± 0.035</td><td>0.157<sup>+0.036</sup><sub>-0.030</sub></td></tr><tr><td>log(P<sub>br</sub>)</td><td>3.32</td><td>3.20 ± 0.17</td><td>3.32<sup>+0.19</sup><sub>-0.24</sub></td></tr><tr><td>p<sub>1</sub></td><td>0.553</td><td>0.592 ± 0.085</td><td>0.70<sup>+0.32</sup><sub>-0.16</sub></td></tr><tr><td>p<sub>2</sub></td><td>-0.75</td><td>-0.99 ± 0.53</td><td>-1.20<sup>+0.92</sup><sub>-1.26</sub></td></tr><tr><td>m</td><td>-0.462</td><td>-0.475 ± 0.075</td><td>-0.46 ± 0.06</td></tr></tbody></table>

III Bern population synthesis model of planetary formation and evolution (Emsenhuber et al. 2020a, 2020b), and we also fit these model results with our Gaussian desert, power-law mass-function model for the three period ranges used in this paper: 3.75–30 days, 30–240 days, and 240–1920 days. Emsenhuber et al. (2020b) present results for different numbers of 0.01 $M_⊕$ planetary embryos in the initial conditions, and they consider their runs with 100 planetary embryos to be their “nominal” or most realistic results. However, they also indicate that their runs with 10, 20, and 50 planetary embryos yield results that are similar to the 100 embryo results, and Fernandes et al. (2019) compared the Mayor et al. (2011) CORALIE/HARPS sample to an early run of the Emsenhuber et al. (2020b) simulations with 20 planetary embryos. So, we have presented comparisons to both the Emsenhuber et al. (2020b) 20 and 100 embryo simulations in Figure 3. The Emsenhuber et al. (2020a, 2020b) simulations consider only solar mass host stars, but this is similar to the Mayor et al. (2011) sample that is dominated by solar-type stars in the mass range $0.7 M_⊙ \le M_{\text{host}} \le 1.3 M_⊙$. The Gaussian desert model parameters for these models are given in Table 4. Each Emsenhuber et al. (2020b) simulation includes 1000 simulated host stars, so it is comparable to the 822 stars in the Mayor et al. (2011) CORALIE/HARPS sample. Our analysis of the simulated planet sample is done in exactly the same way as our analysis of the CORALIE/HARPS sample, with orbital inclinations assigned randomly.

In contrast to the case with the CORALIE/HARPS sample, the Emsenhuber et al. (2020b) simulations show very significant desert features for all period ranges with maximum depths (B values) ranging from $0.838 \pm 0.168$ to $0.974 \pm 0.028$. The widths of these features (in log M) are also larger with $\sigma_{\text{des}}$ values ranging from $0.414 \pm 0.167$ to $0.804 \pm 0.116$. The one peculiar result is for the 100 embryo simulations in the 3.75–30 day period range, where center of the desert occurs above $100 M_⊕$ at $\log M_{\text{des}} = 2.475 \pm 0.236$, but this is largely because this data set has only two of the 98 simulated planets with $M \sin i > 100 M_⊕$.

Since the Mayor et al. (2011) CORALIE/HARPS sample shows no evidence for a planet desert, it is not surprising that the population synthesis results do not match the CORALIE/HARPS data. The simulations indicate the desert features that are not seen in the data, and for the 30–240 day and 240–1920 day period ranges the simulations predict a large excess of $\sim 1000 M_⊕$ planets that are not seen in the observed sample. These discrepancies were also noted by Emsenhuber et al. (2020c). However, they did their comparison using observed and simulated histograms of detected planets instead of the detection efficiency corrected histograms that we present here. The problem with their approach is that a steep rise in planet occurrence at $M \sin i \lesssim 30 M_⊕$ with a flat distribution

above $30 M_⊕$ can be converted into an apparent desert due to the lower detection efficiency for 30–100 $M_⊕$ planets compared to $>100 M_⊕$ planets.

## 6. Discussion and Conclusions

Our analysis has demonstrated that there is no evidence for the predicted sub-Saturn-mass planetary desert (Ida & Lin 2004a) in the Mayor et al. (2011) radial-velocity sample. Previously, the Suzuki et al. (2016) study demonstrated that there was no such planetary desert in the exoplanet mass ratio distribution of planets found by microlensing. These results are complementary because microlensing is very sensitive to low-mass planets (Bennett & Rhie 1996) orbiting beyond the snow line (Gould & Loeb 1992), while the radial-velocity method is currently not sensitive to planets near the lower mass limit of the predicted desert, beyond the snow line. However, the sensitivity of the radial-velocity method improves for shorter period orbits, and the sensitivity of the Mayor et al. (2011) sample extends below the predicted planet desert lower mass limit of $\sim 10 M_⊕$ for orbital periods in our first and second period bins with $P < 240$ days. In our third period bin, at $240 < P < 1920$ days, the sensitivity extends down to only $\sim 12 M_⊕$, but this period range overlaps with the sensitivity of the microlensing method. So, the combination of the radial-velocity analysis presented here and the earlier microlensing work (Suzuki et al. 2016) appears to rule out this predicted planet desert over a wide range of planetary orbits, out to $\sim 10$ au.

While the population synthesis models are often referred to as theoretical models, it might be more reasonable to consider them to be methods that attempt to interpret observational constraints in terms of the many complicated physical processes that are thought to be involved in the formation of planets. As mentioned in the introduction, Suzuki et al. (2018) have presented a number of possible modifications to the Ida & Lin (2004a, 2004b, 2005, 2008a, 2008b) and Bern group (Mordasini et al. 2009a, 2009b, 2015). These include a number of processes that could slow gas accretion or terminate it well before the planet reaches $100 M_⊕$. These include heating of the gaseous envelope by the accretion of planetesimals, a low disk viscosity, low disk scale height, or early formation of a circumplanetary disk. Runaway growth is often thought to terminate at lower masses at wider orbits, so gravitational interactions between the planets could transport lower mass planets to Jupiter-like orbits from the wider orbits where they formed.

Another possibility is that individual planetary systems might form in protoplanetary disks that allow for runaway gas accretion growth that disfavors an intermediate final planetary mass range that could be considered a desert. However, if the properties of these disks have a significant variation, then it could be that these disfavored planetary mass ranges do not line up to create a desert in the combined mass or mass ratio distribution of a large sample of planetary systems. So, it could be that the planet desert feature is simply smoothed out by large variations in the properties of protoplanetary disks. A wide variation in assumed planet accretion rates plays a role in theoretical framework developed by Adams et al. (2021). This approach is more phenomenological than the population synthesis approach, and does not seek to describe most of the detailed physical processes involved in the formation of planets. One of the basic assumptions of the Adams et al. (2021) is
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Comparison of Bern New Generation Planetary Population Synthesis calculations (red histograms) to Mayor et al. (2011) data (blue, cyan, and green histograms) and our Gaussian desert model (black curves and gray shaded regions) using the same three period ranges used in our previous analyses: 3.75–30 days, 30–240 days, and 240–1920 days. Results from simulations starting with 20 0.01 $M_⊕$ planetary embryos are shown in the left panel and results from 100 0.01 $M_⊕$ planetary embryo simulations are shown on the right.</figcaption></figure>

<table><caption>**Table 4**
Exoplanet Desert Power-law Model Parameters for New Bern Models</caption><thead><tr><th rowspan="3">Parameter</th><th colspan="6">Period Ranges for 20 Embryos</th></tr><tr><th colspan="2">3.75–30 days</th><th colspan="2">30–240 days</th><th colspan="2">240–1920 days</th></tr><tr><th>Best</th><th>Average</th><th>Best</th><th>Average</th><th>Best</th><th>Average</th></tr></thead><tbody><tr><td>A</td><td>0.04216</td><td>0.16 ± 0.12</td><td>6.644</td><td>11.8 ± 5.4</td><td>0.1789</td><td>0.132 ± 0.058</td></tr><tr><td>A<sub>Pm</sub></td><td>0.08073</td><td>0.097 ± 0.024</td><td>0.3563</td><td>0.302 ± 0.075</td><td>0.3754</td><td>0.255 ± 0.098</td></tr><tr><td>B</td><td>0.9913</td><td>0.84 ± 0.10</td><td>0.9934</td><td>0.970 ± 0.026</td><td>0.9591</td><td>0.87 ± 0.011</td></tr><tr><td>log M<sub>des</sub></td><td>1.726</td><td>1.83 ± 0.13</td><td>1.707</td><td>1.726 ± 0.070</td><td>1.509</td><td>1.55 ± 0.18</td></tr><tr><td>σ<sub>des</sub></td><td>0.245</td><td>0.41 ± 0.17</td><td>0.850</td><td>0.80 ± 0.12</td><td>0.960</td><td>0.74 ± 0.17</td></tr><tr><td>p</td><td>0.124</td><td>0.03 ± 0.17</td><td>0.927</td><td>1.12 ± 0.19</td><td>-0.687</td><td>-0.63 ± 0.14</td></tr><tr><td>m</td><td>-0.887</td><td>-0.859 ± 0.067</td><td>-0.621</td><td>-0.583 ± 0.073</td><td>-0.334</td><td>-0.22 ± 0.11</td></tr></tbody><tbody><tr><th rowspan="3">Parameter</th><th colspan="6">Period Ranges for 100 Embryos</th></tr><tr><th colspan="2">3.75–30 days</th><th colspan="2">30–240 days</th><th colspan="2">240–1920 days</th></tr><tr><th>Best</th><th>Average</th><th>Best</th><th>Average</th><th>Best</th><th>Average</th></tr><tr><td>A</td><td>0.00580</td><td>0.96 ± 0.65</td><td>12.86</td><td>8.9 ± 2.3</td><td>0.3044</td><td>0.28 ± 0.10</td></tr><tr><td>A<sub>Pm</sub></td><td>0.0877</td><td>0.074 ± 0.028</td><td>0.3941</td><td>0.274 ± 0.066</td><td>0.4348</td><td>0.39 ± 0.12</td></tr><tr><td>B</td><td>0.9993</td><td>0.84 ± 0.17</td><td>0.9820</td><td>0.946 ± 0.028</td><td>0.9996</td><td>0.97 ± 0.028</td></tr><tr><td>log M<sub>des</sub></td><td>2.608</td><td>2.48 ± 0.24</td><td>1.727</td><td>1.758 ± 0.082</td><td>1.708</td><td>1.696 ± 0.082</td></tr><tr><td>σ<sub>des</sub></td><td>0.848</td><td>0.75 ± 0.20</td><td>1.104</td><td>0.80 ± 0.12</td><td>0.757</td><td>0.74 ± 0.14</td></tr><tr><td>p</td><td>-0.081</td><td>0.42 ± 0.20</td><td>1.104</td><td>1.10 ± 0.11</td><td>-0.331</td><td>-0.34 ± 0.14</td></tr><tr><td>m</td><td>-0.896</td><td>-0.918 ± 0.087</td><td>-0.733</td><td>-0.661 ± 0.078</td><td>-0.388</td><td>-0.32 ± 0.10</td></tr></tbody></table>

that the exoplanet mass function has an approximate power-law form that appears to describe the planet distribution above ~10 $M_⊕$ according to studies using a variety of different methods (Cumming et al. 2008; Suzuki et al. 2016; Schlaufman 2018; Nielsen et al. 2019; Wagner et al. 2019). Adams et al. (2021) find that approximately power-law mass functions can be achieved through a combination of exponential decay rate for protoplanetary disk gas and a gas accretion rate roughly

proportional to planet mass if they also include a random gas accretion efficiency factor that can vary significantly between different systems.

One issue that has yet to be explored in much detail is the host mass dependence of the wide-orbit planet distribution. It is possible that the predicted planet desert does exist for some particular range of host masses, although the results presented here tend to disfavor this idea since the Mayor et al. (2011)
<figure><img src="image_6.png" /><figcaption>**Figure 4.** The predicted radius distribution for planets in the three period bins shown in Figure 2 assuming the best-fit power-law models with no desert feature. The small planet desert feature at a radius of $R \approx 9R_⊕$ is due to the mass-radius relation that reflects the pileup of planets at radii of $\sim 14 R_⊕$ due to the dominance of electron degeneracy pressure for objects in the mass range of $0.4 M_{\oplus}$ to $0.08 M_⊕$, a factor of $\gtrsim 200$ in mass.</figcaption></figure>

sample consists of FGK stars, while the Suzuki et al. (2016) sample extends to lower mass host stars. Nevertheless, Kepler observations have established that planets in short period orbits have a higher occurrence rate around low-mass stars than around solar-type stars (Mulders et al. 2015), so it would not be surprising to find a host mass dependence of the wide-orbit planet-occurrence rate. However, most of the existing demographic studies do not allow this. The exiting microlensing studies use only the planet-star mass ratio (Gould et al. 2010; Cassan et al. 2012; Suzuki et al. 2016), and the largest of the radial-velocity studies (Mayor et al. 2011) does not provide enough information for such a study. The recent radial-velocity analysrs of Wittenmyer et al. (2020) and Fulton et al. (2021) do provide the individual detection efficiencies for both the stars with and without detected planets, which is what is needed. However, the Wittenmyer et al. (2020) study and both survey arenot as sensitive to low planet masses both surveys are not as sensitive to low planet masses as the Mayor et al. (2011) study. So, this data set is less suitable for a study of the predicted sub-Saturn mass desert for wide-orbit planets.

While most of the parameters of planetary microlensing events can only be determined with data taken during the event, it is possible to determine the masses of the host stars and planets and their separation with high angular resolution follow-up observations from the Hubble Space Telescope (Bennett et al. 2015; Bhattacharya et al. 2017) or adaptive optics (AO) observations on 8–10 m class telescopes. The Keck telescopes have been the most successful at this, due an AO system that is more effective with faint guide stars (Batista et al. 2015; Beaulieu et al. 2016; Bhattacharya et al. 2018, 2020; Bennett et al. 2020; Vandorou et al. 2020; Terry et al. 2021).

A good AO correction is important for the interpretation of these follow-up observations because it is generally necessary to confirm the identification of the host star by confirming the relative proper motion between the candidate host star and the background source star (Koshimoto et al. 2020). The MOA Collaboration is in the process of extending the Suzuki et al. (2016) statistical sample, while the Keck and Hubble follow-up program is expanding with observations and host and planet-mass determinations for many more microlens planetary systems in this extended statistical sample. New radial velocity have recently been reported from the California Legacy Survey (Fulton et al. 2021). So, we expect significant developments in the understanding of wide-orbit exoplanet demographics in the next few years. Ultimately, the Nancy Grace Roman Telescope’s Galactic Exoplanet survey (Bennett & Rhie 2002; Penny et al. 2019) will perform a much more comprehensive statistical study of the occurrence rate of wide-orbit planets at semimajor axes of $a \gtrsim 0.5$ au with masses down to the mass of Mars ($0.1 M_⊕$).

We thank Fred Adams, B. Scott Gaudi, Gijs Mulders, and Ilaria Pascucci for helpful suggestions and comments on a draft version of this paper. D.P.B. and C.R. were supported by NASA through grant NASA-80NSSC18K0274 and award number 80GSFC17M0002. C.R. is supported by the ANR project COLD-WORLDS of the French *Agence Nationale de la Recherche* with the reference ANR-18-CE31-0002. The results reported herein have benefited from collaborations and/or information exchange within NASA’s Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA’s Science Mission Directorate.
# Appendix
## Planet Deserts in Kepler Data?

There have been extensive discussions in the literature of planet deserts seen in Kepler data. As mentioned in the introduction, there is the “hot-Neptune desert”, also known as the “sub-Jovian desert” (Szabó & Kiss 2011; Lundkvist et al. 2016; Mazeh et al. 2016; Owen 2019). This was an observational discovery, based primarily on Kepler data, and it is thought to be due to heating of the atmospheres of planets that come very close to their host stars. There is also the Fulton gap (Fulton et al. 2017), which is sometimes referred to as the “sub-Neptune” desert and is thought to be caused by photoevaporation. These deserts occur for only very short period orbits and are thought to be caused by heating or photoevaporation due to the close proximity of the host star, so they are not relevant to the desert that the topic of this paper, which is predicted by the runaway gas accretion scenario at much wider orbits.

There has also been some, mostly informal, suggestions that the Kepler data shows a sub-Jupiter radius desert for cooler planets at wider separations. For example, this shows up, somewhat weakly, in Figure 7 of Thompson et al. (2018) at orbital periods ≳20 days, and it also can be seen in Figure 3 of Berger et al. (2020) at an incident flux of ≲10× that received by the Earth. In this section, we show that this sub-Jupiter radius gap is likely to be real, but that it does not imply a sub-Jupiter-mass gap.

In order to relate the masses and radii of planets, we use the mass–radius relation of Chen & Kipping (2017), and to describe the mass distribution of exoplanets we employ the best-fit pure power-law models from Table 2 for our three different period ranges: 3.75–30 days, 30–240 days and 240-1920days. The results of this calculation are shown in Table 4. This shows a clear minimum in the planetary radius distribution at $R \simeq 9R_⊕$ and a maximum at $R \simeq 14R_⊕$. This feature is dominated by planets in the longest period bin of 240 days $< P \le 1920$ days because the assumed long period mass function has the shallowest mass dependence: $M^m$ with $m = -0.367$. Note that this exercise is intended to show how a smooth mass function can still imply a “planet-radius desert” and we do not intend to claim that these exoplanet mass functions can provide a good description of the Kepler data with its much larger exoplanet sample.

The reason that the smooth, monotonic mass function is converted into a radius function with a dip and a peak is clear from Figure 3 of Chen & Kipping (2017), which displays their mass–radius relation. This relation is nearly flat between masses of $0.4 M_{Jup}$ to $0.08 M_⊕$, a factor of 209 in mass. This is a well known feature caused by the fact that gas-giant planets, brown dwarfs, and the lowest mass stars are all largely supported by electron degeneracy pressure (which is also responsible for the decreasing radius with increasing mass for white dwarfs). Thus, the observed potential sub-Jupiter desert in the Kepler planetary radius distribution is likely to be explained by the mass–radius relation rather than a feature in the exoplanet mass distribution, since the feature appears for a model that has no desert in the mass distribution.

# ORCID iDs

David P. Bennett https://orcid.org/0000-0001-8043-8413

Clément Ranc https://orcid.org/0000-0003-2388-4534

# References

Adams, F. C., Meyer, M. R., & Adams, A. D. 2021, ApJ, 909, 1

Alcock, C., Allsman, R. A., Axelrod, T. S., et al. 1996, ApJ, 461, 84

Batista, V., Beaulieu, J.-P., Bennett, D. P., et al. 2015, ApJ, 808, 170

Beaulieu, J.-P., Bennett, D. P., Batista, V., et al. 2016, ApJ, 824, 83

Bennett, D. P., Bhattacharya, A., Anderson, J., et al. 2015, ApJ, 808, 169

Bennett, D. P., Bhattacharya, A., Beaulieu, J.-P., et al. 2020, AJ, 159, 68

Bennett, D. P., & Rhie, S. H. 1996, ApJ, 472, 660

Bennett, D. P., & Rhie, S. H. 2002, ApJ, 574, 985

Berger, T. A., Huber, D., Gaidos, E., et al. 2020, AJ, 160, 108

Bhattacharya, A., Beaulieu, J.-P., Bennett, D. P., et al. 2018, AJ, 156, 289

Bhattacharya, A., Bennett, D. P., Anderson, J., et al. 2017, AJ, 154, 59

Bhattacharya, A., Bennett, D. P., Beaulieu, J. P., et al. 2021, AJ, 162, 60

Cassan, A., Kubas, D., Beaulieu, J.-P., et al. 2012, Natur, 481, 167

Chen, J., & Kipping, D. 2017, ApJ, 834, 17

Cumming, A., Butler, R. P., Marcy, G. W., et al. 2008, PASP, 120, 531

Dobbs-Dixon, I., Li, S. L., & Lin, D. N. C. 2007, ApJ, 660, 791

Emsenhuber, A., Mordasini, C., Burn, R., et al. 2020a, arXiv:2007.05561

Emsenhuber, A., Mordasini, C., Burn, R., et al. 2020b, arXiv:2007.05562

Emsenhuber, A., Mordasini, C., Mayor, M., et al. 2020c, EPSC, 14, 339

Fernandes, R. B., Mulders, G. D., Pascucci, I., et al. 2019, ApJ, 874, 81

Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109

Fulton, Benjamin J., Rosenthal, Lee J., Hirsch, Lea A., et al. 2021, ApJS, 255, 14

Fung, J., Shi, J.-M., & Chiang, E. 2014, ApJS, 782, 88

Garaud, P., & Lin, D. N. C. 2007, ApJ, 654, 606

Gould, A., Dong, S., Gaudi, B. S., et al. 2010, ApJ, 720, 1073

Gould, A., & Loeb, A. 1992, ApJ, 396, 104

Ida, S., & Lin, D. N. C. 2004a, ApJ, 604, 388

Ida, S., & Lin, D. N. C. 2004b, ApJ, 616, 567

Ida, S., & Lin, D. N. C. 2005, ApJ, 626, 1045

Ida, S., & Lin, D. N. C. 2008a, ApJ, 673, 487

Ida, S., & Lin, D. N. C. 2008b, ApJ, 685, 584

Jung, Y. K., Gould, A., Zang, W., et al. 2019, AJ, 157, 72

Koshimoto, N., Bennett, D. P., & Suzuki, D. 2020, AJ, 159, 268

Lundkvist, M. S., Kjeldsen, H., Albrecht, S., et al. 2016, NatCo, 7, 11201

Mayor, M., Marmier, M., Lovis, C., et al. 2011, arXiv:1109.2497

Mazeh, T., Holczer, T., & Faigler, S. 2016, A&A, 589, A75

Mordasini, C., Alibert, Y., & Benz, W. 2009a, A&A, 501, 1139

Mordasini, C., Alibert, Y., & Benz, W. 2009b, A&A, 501, 1161

Mordasini, C., Mollière, P., Dittkrist, K.-M., Jin, S., & Alibert, Y. 2015, IJASB, 14, 201

Mulders, G. D., Pascucci, I., & Apai, D. 2015, ApJ, 798, 112

Nayakshin, S., Dipierro, G., & Szulágyi, J. 2019, MNRAS, 488, L12

Nielsen, E. L., De Rosa, R. J., Macintosh, B., et al. 2019, AJ, 158, 13

Owen, J. E. 2019, AREPS, 47, 67

Penny, M. T., Gaudi, B. S., Kerins, E., et al. 2019, ApJS, 241, 3

Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icar, 124, 62

Schlaufman, K. C. 2018, ApJ, 853, 37

Suzuki, D., Bennett, D. P., Ida, S., et al. 2018, ApJL, 869, L34

Suzuki, D., Bennett, D. P., Sumi, T., et al. 2016, ApJ, 833, 145

Szabó, G. M., & Kiss, L. L. 2011, ApJL, 727, L44

Szulágyi, J., Morbidelli, A., Crida, A., et al. 2014, ApJ, 782, 65

Terry, S. K., Bhattacharya, A., Bennett, D. P., et al. 2021, AJ, 161, 54

Thompson, S. E., Coughlin, J. L., Hoffman, K., et al. 2018, ApJS, 235, 38

Trotta, R. 2008, ConPh, 49, 71

Udalski, A., Ryu, Y.-H., Sajadian, S., et al. 2018, A&A, 68, 1

Vandorou, A., Bennett, D. P., Beaulieu, J.-P., et al. 2020, AJ, 160, 121

Wagner, K., Apai, D., & Kratter, K. M. 2019, ApJ, 877, 46

Wittenmyer, R. A., Wang, S., Horner, J., et al. 2020, MNRAS, 492, 377