The California Legacy Survey II. Occurrence of Giant Planets Beyond the Ice line

We used high-precision radial velocity measurements of FGKM stars to determine the occurrence of giant planets as a function of orbital separation spanning 0.03-30 au. Giant planets are more prevalent at orbital distances of 1-10 au compared to orbits interior or exterior of this range. The increase in planet occurrence at $\sim$1 au by a factor of $\sim$4 is highly statistically significant. A fall-off in giant planet occurrence at larger orbital distances is favored over models with flat or increasing occurrence. We measure $14.1^{+2.0}_{-1.8}$ giant planets per 100 stars with semi-major axes of 2-8 au and $8.9^{+3.0}_{-2.4}$ giant planets per 100 stars in the range 8-32 au, a decrease in giant planet occurrence with increasing orbital separation that is significant at the $\sim$2$\sigma$ level. We find that the occurrence rate of sub-Jovian planets (0.1-1 Jupiter masses) is also enhanced for 1-10 au orbits. This suggests that lower mass planets may share the formation or migration mechanisms that drive the increased prevalence near the water-ice line for their Jovian counterparts. Our measurements of cold gas giant occurrence are consistent with the latest results from direct imaging surveys and gravitational lensing surveys despite different stellar samples. We corroborate previous findings that giant planet occurrence increases with stellar mass and metallicity.


INTRODUCTION
Expanding and characterizing the population of known exoplanets with measured masses and orbital periods is crucial to painting a more complete picture of planet formation and evolution. A census of diverse exoplanets sheds light on worlds radically different from Earth and can provide insight into how these planetsand those orbiting the Sun-formed. Ground-based ra-dial velocity (RV) surveys measure the Doppler shifts of stellar spectra to discover exoplanets and characterize their orbits and masses. These surveys have provided landmark discoveries that shaped our understanding of the formation and architectures of other worlds (e.g., Mayor & Queloz 1995;Marcy et al. 2002;Tamuz et al. 2008).
Doppler planet searches take time to accumulate the time series measurements that trace out planetary orbits. The Keck Planet Survey (Cumming et al. 2008) used eight years of RVs from Keck-HIRES (Vogt et al. 1994) to make the first broad measurement of giant planet occurrence (M sin i ≥ 0.1M J ). This survey discovered an increase in the abundance of giant planets for orbits near the water-ice line and found that about 10% of Sun-like stars have giant planets with a semimajor axes of <3 au. The survey only reported planet detections for orbital periods shorter than 2000 days, the observational baseline of the survey. Extrapolating based on the detection of partial orbits, Cumming et al. (2008) estimated that ∼20% of such stars have a giant planet orbiting within 20 au.
Other teams of astronomers have surveyed the Northern and Southern skies in parallel with the Keck search. Mayor et al. (2011) used 8 years of precise HARPS RVs supplemented by additional RVs from CORALIE to measure occurrence patterns in the population of giant planets that are similar to those described above. They found that the planet mass function is "bottom heavy". That is, low-mass planets (0.3-30 M ⊕ ) are significantly more common than giant planets, a finding consistent with measurements from Keck Observatory by . Since then, the HARPS team has continued to discover increasingly longer-period and lowermass planets (Udry et al. 2017;Rickman et al. 2019). Two other 'legacy' planet searches have contributed significantly to our knowledge of giant planets. Wittenmyer et al. (2020a) used data from a subset of the stars surveyed by the 18-year Anglo-Australian Planet Search, which has also uncovered a number of cold giant planets (Wittenmyer et al. 2017;Kane et al. 2019), to measure a significant increase in giant planet occurrence at ∼1 au and a constant occurrence for orbits in the range ∼1-6 au. Similarly, the McDonald Observatory planet search has been operating for more than 20 years using the 2.7-m Harlan J. Smith Telescope, and has contributed valuable discoveries of long-period giant planets (e.g., Robertson et al. 2012;Endl et al. 2016;Blunt et al. 2019).
We are now in the fourth decade of Doppler planet searches. As we begin to discover planets with orbital periods comparable to Saturn's, we can answer ques-tions that require a rigorous accounting of giant planets spanning a large range of orbital distances. What is the mass versus semi-major axis distribution of planets out to 10 au? How abundant are cold gas giants beyond the water-ice line, and what can this abundance tell us about planet formation across protoplanetary disks?
The California Legacy Survey (CLS, Rosenthal et al. 2021) is uniquely suited for this work. As an unbiased radial velocity survey of 719 stars over three decades, the CLS is an excellent sample for a variety of occurrence measurements, particularly for cold gas giants. In this paper, we explore giant planet occurrence as a function of orbital separation. In Section 2, we review the star and planet catalog of the California Legacy Survey. In Section 3, we describe our methods for computing planet occurrence. Section 4 describes the patterns of planet occurrence that we observe in the sample. In Section 5, we discuss our findings and their context. We summarize our work in Section 6.

SURVEY REVIEW
The California Legacy Survey is a Doppler search for planets orbiting a well-defined sample of nearby FGKM stars conducted by the California Planet Search team (CPS; . Paper I in this series (Rosenthal et al. 2021) describes the CLS in detail, including the stellar sample, the search methodology, and the resulting planet sample upon which this paper and forthcoming works in the CLS paper series build. The CLS stellar sample was selected specifically to make the measurements reported here-planet occurrence measurements, especially of giant planets with orbits out to 10 au and beyond-and it approximates a random sample of nearby stars. In particular, stars were selected for CLS observations independent of whether planets were known to orbit them. Stars were also selected independent of their metallicity or other factors that might make them more or less likely to harbor planets.
CLS builds on Doppler measurements from the Keck Planet Search (Cumming et al. 2008), a touchstone Doppler survey of 585 stars observed with HIRES at the W. M. Keck Observatory during 1996-2004. We continued to observe those stars and an additional 134 stars at Keck Observatory through 2020. CLS also includes observations of a subset of these stars made with the Hamilton spectrometer at Lick Observatory during 1988-2011, high-cadence Keck-HIRES observations of 235 magnetically inactive stars as part of the Eta-Earth Survey , and high-cadence Lick-APF observations of 135 of those stars (Fulton et al. 2016;Hirsch et al. 2021). The average star has been observed for 22 years and has 71 RVs with a precision of ∼2 m s −1 . While these stars do not have homogeneous observing histories, our search methodology accounts for this by incorporating the search completeness of each star's individual dataset. (A Doppler survey that is completely homogeneous in the number, precision, and temporal spacing of measurements is infeasible given the three decade history of this planet searchindeed, this survey spans an era longer than the time during which extrasolar planets orbiting Sun-like stars have been known!) By the metric of "Doppler survey étendue" (number of stars surveyed × typical time series duration), CLS is the largest planet search to date at the ∼m s −1 level.
Our search methodology (described in Rosenthal et al. 2021) involves an automated, iterative periodogrambased search for Keplerian signals with uniform vetting to identify false positives. This methodology detected 177 planets orbiting the 719 stars in the CLS stellar sample. The algorithm is sensitive to orbital periods much longer than the baseline of our dataset, with the longest period signals detected as partial orbits.
The search was also sensitive to orbital segments only seen as linear and parabolic trends in an RV time series. There were only six such detections in our sample of trends that are not associated with known stellar binaries and are potentially consistent with planetary mass companions. Thus, nearly all orbital signals were resolved or partially resolved as Keplerian signals.
To characterize survey completeness for each star in the survey, we conducted injection-recovery tests of synthetic Doppler planet signals over a range of injected masses, orbital periods, and orbital geometries. Detected planets and CLS survey completeness are shown in Figure 1. We refer the reader to Rosenthal et al. (2021) for the full stellar sample and planet catalog.
The CLS stellar sample has a median metallicity of [Fe/H] = 0.0 dex, a median stellar mass of 1.0 M , and a small number of evolved stars (subgiants). These are good heuristics for verifying that we successfully constructed a blind occurrence survey, since a bias toward known giant planet hosts could manifest as a metal-rich sample (Fischer & Valenti 2005;Santos et al. 2004), a particularly massive sample, or an excess of evolved stars (Johnson et al. 2011).

METHODS
The primary goal of this work is to measure planet occurrence. Many studies of RV or transit surveys use the intuitive occurrence measurement method known as "inverse detection efficiency" (Howard et al. 2012;Petigura et al. 2013). According to this procedure, one estimates occurrence in a region of parameter space by counting the planets found in that region, with each planet weighted by the local search completeness. One can measure the search completeness map of a survey by injecting many synthetic signals into each dataset and computing the fraction of signals in a given region that are recovered by the search algorithm in use. Inverse detection efficiency is actually a specific case of a Poisson likelihood method, in which one models an observed planet catalog as the product of an underlying Poisson process and empirical completeness map (Foreman-Mackey et al. 2014). This can be done with a parametric occurrence rate density model, like a broken power law, or a non-parametric density model, with a piecewise-constant step function. In this paper, we used the Poisson likelihood method to model the occurrence of giant planets, taking measurement uncertainty into account.
We used the hierarchical Bayesian methodology outlined in Hogg et al. (2010) and Foreman-Mackey et al. (2014) to evaluate our occurrence likelihood. Given an observed population of planets with orbital and M sin i posteriors {ω} and associated survey completeness map Q(ω), and assuming that our observed planet catalog is generated by a set of independent Poisson process draws, we evaluated a Poisson likelihood for a given occurrence model Γ(ω|θ), where Γ is an occurrence density (1) The Poisson likelihood for an observed population of objects is where K is the number of observed objects, and ω k is a vector of parameters that completely describe the kth planet's orbit. In our case, the two relevant parameters are M sin i and semi-major axis a, taken from the broader set that includes eccentricity, time of inferior conjunction, and argument of periastron. The Poisson likelihood can be understood as the product of the probability of detecting an observed set of objects (the product term in Equation 2) and the probability of observing no additional objects in the considered parameter space (the exponentiated integral). Equations 1 and 2 serve as the foundation for our occurrence model but do not take into account uncertainty in measurements of planetary orbits and minimum masses. In order to do this, we used RadVel and emcee to empirically sample the orbital posteriors of each system (Fulton et al. 2018;Foreman-Mackey et al. 2013). We hierarchically modeled the orbital posteriors of each planet in our catalog by summing our occurrence model over many posterior samples for each planet. The hierarchical Poisson likelihood is therefore approximated as where N k is the number of posterior samples for the kth planet in our survey and ω n k is the nth sample of the kth planet's posterior. p(ω|α) is our prior on the individual planet posteriors. We placed linear-uniform priors on M sini and log-uniform priors on a. We used emcee to sample our hierarchical Poisson likelihood.
We used two different occurrence frameworks to model our planet population. The first is a non-parametric model across bins uniformly spaced in ln(M sini) and ln(a), with a set of steps ∆ of height θ. We define this framework with the occurrence function The second framework is a broken power law as a function of semi-major axis, defined with the function where C is a normalization constant, β is the occurrence power law index beyond the breaking point, a 0 determines the semi-major axis location of the breaking point, and β + γ is the power law index within the breaking point. This model assumes a giant planet mass function that does not change with respect to semi-major axis. We fit this model to our population in order to explore whether giant planet occurrence falls off beyond the water-ice line. Figure 2 shows occurrence rates as a function of semimajor axis for planets with masses between 30 M ⊕ and 6000 M ⊕ , derived using the non-parametric model described in §3 and assuming uniform occurrence across ln(M sin i). We confirmed the previous result from . The transparent curves represent draws from the broken power law posterior. We find that the power law index beyond the break is ∼2.5σ-separated from zero, implying an occurrence falloff beyond the water-ice line. Cumming et al. (2008) performed a power-law fit to the occurrence rates of planets orbiting only within 3 au; the light dotted blue line represents an extrapolation to wider separations.

Enhancement for giant planets
more massive than 30 M ⊕ are 2-4 times more common at orbital distances between 1-3 au relative to 0.1-0.3 au. Using our broken power law model, we find a median power law slope inside the break of 0.72 +0.16 −0.20 , which is 2 σ higher than the power law slope measured by Cumming et al. (2008) (0.26±0.1). This difference is likely caused by the single power law model being pulled to lower values due to neglecting a flattening or turnover in occurrence at long orbital periods since Cumming et al. (2008) was limited to planets orbiting inside 3 au.

Distribution of giant planets beyond 3 au
Due to low completeness beyond our observational baselines, our occurrence results beyond 10 au are highly uncertain. However, we can estimate occurrence trends with the broken power law model described in §3. Figure  3 shows the broken power law results juxtaposed with the non-parametric results, and Figure 4 presents the posteriors for the parametric model parameters. The medians and 68th percentile credible intervals for the broken power law model are listed in Table 1. Both assume uniform occurrence across ln(M sin i). We find that 99.4% of the posterior samples are consistent with a plateauing or declining occurrence rate beyond a peak around 3.6 +2.0 −1.8 au. We find that the power law index beyond the peak is β = −0.86 +0.41

Comparing sub-and super-Jovians
Figure 5 compares non-parametric occurrence rates for giant planets more and less massive than 300 M ⊕ . We find a quantitatively similar occurrence enhance-ment around 1-10 au for both the sub-Jovian-mass and Jovian-mass planets. However, we lack the sensitivity to measure the occurrence rate of sub-Jovian mass planets beyond 10 au, to assess whether they exhibit the fall-off in occurrence at large orbital separations seen when examining occurrence across both mass ranges. The sub-Jovian planets are more common than the super-Jovian planets across a wide range of separations, particularly beyond the water ice line. We find a similar enhancement for sub-Saturns below 150 M ⊕ , implying that this occurrence enhancement is independent of planet mass. Combining these two populations produces the same trends seen when we assume uniform occurrence across all masses.
We more concretely measured occurrence as a function of mass by performing a non-parametric fit to our sample within 1-5 au. Figure 6 shows occurrence as a function of M sin i within 30-3000 M ⊕ , in four steps. This figure shows that our assumption of a uniform ln(M sin i) distribution beyond the ice line is valid up to 900 M ⊕ , but the distribution falls off with ∼2σ significance above 900 M ⊕ . If this is also true beyond 5 au, where low completeness prevents us from making a similar measurement, then we may be underestimating broad giant planet occurrence in our lowest-completeness region of parameter space, beyond 10 au. This is because our only detections in that regime are more massive than 300 M ⊕ , and all but one of them are more massive than 900 M ⊕ .

Occurrence with respect to stellar mass and metallicity
In addition to measuring occurrence with respect to semi-major axis and M sin i, we measured the broad occurrence rate of giant planets more massive than 100 M ⊕ and within 1-5 au with respect to host-star mass and metallicity. We chose a lower limit of 100 M ⊕ instead of 30 M ⊕ in order to restrict our analysis to searchcomplete regions within 1-5 au, since 30 M ⊕ planets are effectively undetectable beyond 3 au. For each of these two stellar properties, we computed occurrence across six divisions, in steps of 0.2 M across 0.3-1.5 M and 0.15 dex across -0.5-0.4 dex respectively. Figure 7 shows occurrence with respect to host-star mass, while Figure 8 shows occurrence with respect to host-star [Fe/H]. Both of our measurements agree with prior results. Johnson et al. (2010), whose stellar sample was excluded from CLS due to its bias toward giant planet hosts, measured giant planet occurrence across stellar mass and found an increase in occurrence with increasing stellar mass beginning near 1 M . Wittenmyer et al. (2020b) independently found an increase in giant planet occurrence beyond 1 M . We see the same phenomenon in our sample, as presented in Figure 7. Similarly, Fischer & Valenti (2005) found that giant planet occurrence increases with increasing [Fe/H] beyond 0.1 dex, as did Reffert et al. (2015) and Jones et al. (2016). We see the same transition near 0.1 dex in Figure 8. presented a HARPS (Mayor et al. 2003) and CORALIE (Baranne et al. 1996) blind radial velocity survey of 822 stars and 155 planets over 10 years (corresponding to a 4.6 au circular orbit around a Solar-mass star). Mayor et al. (2011), who did not publish their HARPS and CORALIE RVs, measured giant planet occurrence as a function of orbital period out to 4000 days, in the range of the water-ice line. Fernandes et al. (2019) pushed out to low-completeness regimes and estimated a sharp falloff in occurrence beyond the water-ice line. They measured an integrated occurrence rate of 1.44±0.54 giant planets (0.1-20 M J ) per 100 stars for separations between 3.8 and 7.1 au. Our results indicate a much higher occurrence rate for the same planets at those separations; 15.5 +3.2 −3.0 giant planets per 100 stars. The treatment of partial orbits in Mayor et al. (2011) is unclear, and they only measured occurrence with respect to orbital period out to 3000 days (∼4 au). If Mayor et al. (2011) under-reported partial orbits beyond this period in their sample or overestimated sensitivity to partial orbits, then that could explain the large discrepancy between this work and Fernandes et al. (2019) at separations beyond 10 au.
In contrast, Wittenmyer et al. (2020a), which drew from the Anglo-Australian Planet Search (Tinney et al. 2001) to construct a blind survey of 203 stars and 38 giant planets over 18 years, found that giant planet occurrence is roughly constant beyond the water-ice line, out to almost 10 au. Wittenmyer et al. (2020a) reports an occurrence rate of 6.9 +4.2 −2.1 giant planets > 0.3 M J per 100 stars with periods between 3000 and 10,000 days (≈4-9 au). Our integrated occurrence rate in the same region of parameter space is slightly higher at 12.6 +2.6  -Mackey et al. (2016) performed an automated search for long-period transiting exoplanets in a set of archival Kepler light curves of G and K stars. For planets between 1.5-9 au and 0.01-20 M J and using a probabilistic mass-radius relationship drawn from Chen & Kipping (2016), they found an occurrence rate density of d 2 N dln(a)dln(M ) = 0.068 ± 0.019. We applied our occurrence model to the same parameter space and found

Comparison to direct imaging surveys
RV surveys have recently begun to approach baselines long enough to detect and place limits on the frequency of planets like those detected by direct imaging. One caveat is that direct imaging surveys usually target stars younger than 100 Myr, while RV surveys generally target stars older than 1 Gyr. Young planets retain significant heat from their formation and are bright in the infrared wavelengths covered by direct imaging surveys. However, young stars also tend to be active and rapidly rotating, which makes precise RV work difficult. Because of this, there is minimal overlap between planets that have been detected by direct imaging and planets that have been detected by radial velocity measurements.
We can still compare rates across these detection methods by making the assumption that giant planet occurrence does not change as host stars age beyond ∼10 Myr, once protoplanetary disks have dissipated. We compared our occurrence model to the results of two direct imaging surveys of nearby stars. Biller et al. (2013) imaged 80 stars in nearby moving groups and detected a small number of brown dwarf companions but no planetary-mass companions. They used stellar evolution and planet formation models to estimate constraints on cold giant occurrence from their nondetections and sensitivity. More recently, Nielsen et al. (2019) imaged 300 stars and detected six planets and three brown dwarfs. Figure 9 compares these results to our occurrence measurements in their respective regions of parameter space. Our measurements are compatible with the limits placed on planets with masses 1-20 M J and separations between 10-50 au by Biller et al. (2013), depending on their assumed stellar evolutionary model that determines the expected brightness of young giant planets. Our measurement for planets with masses 5-14 M J orbiting between 10-100 au is in excellent agreement with the results of Nielsen et al. (2019). The only shared quality of our modeling methods is a Poisson counting likelihood. With the caveat of small number statistics, this is a remarkable benchmark for comparing exoplanet occurrence across independent search methods.

Comparison to gravitational microlensing surveys
We compare our model to the microlensing surveys of Cassan et al. (2012) and Clanton & Gaudi (2016). Like all gravitational lensing surveys, these studies assume a broad prior for stellar type based on Galactic observations, a prior that peaks in the M dwarf range. Our planet-hosting stars have a much higher median mass than this range, but since the gravitational lensing estimates comes purely from a galactic model prior, we chose to perform this broad comparison across stellar masses with the knowledge that the mass range for the lensing numbers is poorly constrained. Figure 10 shows that our estimates agree with broad constraints from the pure lensing survey (Cassan et al. 2012). On the other hand, the constraints of Clanton & Gaudi (2016) strongly disagree with our planet occurrence measurement in the same parameter box. This may be due 5.5. Implications for planet formation Cumming et al. (2008) first identified an enhancement in the occurrence rate of giant planets beyond orbital periods of ∼300 days. We expect such enhancements based on planetary migration models (Ida & Lin 2004). The orbital period distribution in Cumming et al. (2008) predicted a smooth rise in occurrence toward longer or-bital periods, but we observed a sharp transition around 1 au, as seen in Figure 2. Ida & Lin (2008) later suggested that additional solid materials due to ices in the protoplanetary disk could augment the formation of gas giant planets and cause a rapid rise in the occurrence rate of these planets beyond the water-ice line.
If increased solids near and beyond the ice line cause a sharp rise in the occurrence rate, then we might expect this rise to be more well-defined when looking in a unit more closely related to the temperature in the protoplanetary disk. In Figure 11, we plot the occurrence rate as a function of stellar light intensity relative to Earth. The occurrence rate with respect to flux is qualitatively similar to the rate with respect to orbital separation. We do not see strong evidence that the occurrence rate enhancement is any more localized in terms of stellar light intensity relative to Earth. The decline in occurrence for fluxes less than that received by the Earth from the Sun looks more convincing, but all except for the one bin with the highest occurrence near 1 au are consistent with a constant occurrence rate.
We can separate the puzzle of gas giant formation into two components: the growth of solid cores that are large enough to undergo runaway gas accretion, and the process of gas accretion onto solid cores. It is currently unclear whether giant planet occurrence increases beyond the ice line because cores form more easily in this region, or because conditions are more favorable for rapid gas accretion onto solid cores. A number of studies (e.g. Morbidelli et al. 2015;Schoonenberg & Ormel 2017;Drążkowska & Alibert 2017) have argued that that core growth should be enhanced beyond the ice line. If the solid grain sizes and densities beyond the ice line are enhanced during the earliest stages of planet formation, it would facilitate pebble clumping that leads to planetesimal formation and also result in higher pebble accretion rates onto the growing core (e.g. Bitsch et al. 2019).
It is also possible that gas giants are more common beyond the ice line because it is easier for cores to rapidly grow their gas envelopes in this region. The rate at which the growing planet's envelope can cool and contract (hence accreting more gas) depends sensitively on its envelope opacity (e.g. Bitsch & Savvidou 2021). In a recent study, Chachan et al. (2021) used dust evolution models to study the effect of dust opacity and dust-to-gas ratio on giant planet formation in the epoch immediately following the end of core formation. They found that as the disk evolves, decreasing dust opacity beyond the water-ice line allows for higher gas accretion rates in this region. Ida et al. (2018) recently updated their models with an improved treatment of Type II migration. This mecha-  Figure 2, occurrence with respect to stellar light intensity instead of orbital separation. Here we see a similar enhancement in the occurrence rate of giant planets where the insolation flux is equal to that of Earth and tentative evidence for a fall off in occurrence just beyond that.
nism would produce a broad semi-major axis distribution with many giant planets migrating inward to separations less than 1 au. However, Fernandes et al. (2019) show that this model does not agree well with the occurrence contrast between the peak and the inner regions of these systems. Our results are in close agreement with those of Fernandes et al. (2019) for separations less than 3 au. The multi-core accretion models of Mordasini (2018) are also in good agreement with the overall shape of the semi-major axis distribution, but they underestimate the absolute occurrence rate of giant planets. This could be due to the finite number of cores injected into their simulations.
One common theme among planet formation models of gas giants is that protoplanets tend to migrate inward, all the way to the inner edge of the disk, on timescales much shorter than the gas dissipation timescale. This tends to produce an enhancement of occurrence closer to the star and/or many planets being engulfed by the host star. Jennings et al. (2018) attempt to solve this issue by simultaneously modeling the effects of photoevaporation and viscous evolution on the gas disk. They find that, depending on the dominant energy of the evaporating photons, this could clear gaps in the disk that halt Type I migration and creates a pile-up of planets at orbital separations between 0.8-4 au. They showed that this can produce very strong and narrow enhancements near certain orbital separations, but it is conceivable that the shape of the final semi-major axis distribution would actually be driven by the spectral energy distributions of host stars during the early years of their formation.
Hallatt & Lee (2020) also proposed gap formation in the protoplanetary disk shortly after the formation of gas giant planets as a mechanism to slow or halt migration at preferred orbital separations. Their model requires that the giant planets that form further out in the disk are more massive in order to reproduce the observed enhancements. We expect this to be the case if the dust content of disk envelopes is very low.
The observed enhancement in the occurrence rate of sub-Jovian planets near 1-10 au, seen in Figure 5, suggests that the processes that drive the formation and pile-up of planets at those orbital distances also apply to these lower-mass planets. It appears just as likely for a gaseous planet to undergo runaway accretion and grow into a Jovian planet as it is to halt that runaway accretion process early and remain in the sub-Saturn regime.
Unfortunately, it is difficult to extract significant constraints on planet formation models from semi-major axis distributions alone. Future planet catalogs produced by Gaia and The Roman Space Telescope will help to measure the precise shape of the occurrence enhancement around 1 au with planet samples several orders of magnitude larger, but the stellar samples will be different from ours. We plan for future works in this series to analyze the host star metallicity, eccentricity, and multiplicity distributions of our sample, in the hopes of uncovering evidence that discriminates between different planet formation models.

CONCLUSION
In this work, we utilize the catalogue of stars, RV-detected planets, and completeness contours from Rosenthal et al. (2021) to measure giant planet occurrence as a function of semi-major axis. We applied a hierarchical Bayesian technique to incorporate measured search completeness and uncertainties in our observations into uncertainties in our occurrence rates. Our results are consistent with previous studies that have found a strong enhancement in the occurrence rates of these planets around 1 au.
We find that the occurrence of planets less massive than Jupiter (30 ≤ M sin i≤ 300 M ⊕ ) is enhanced near 1-10 au in concordance with their more massive counterparts. We find that a fall-off in giant planet occurrence at larger orbital distances is favored over models with flat or increasing occurrence, with 2.5 σ confidence from our broken power-low model and with 1.5 σ confidence from our non-parametric model. Additionally, our occurrence measurements beyond 10 au are consistent with with those derived from direct imaging surveys.
Finally, we lay out the methodology and groundwork for future studies of giant occurrence as a function of planet and host-star properties. With these tools, we plan to study the occurrence rates of giant planets in particular configurations. Paper III in the CLS series will examine the relationship between giant planets and smaller companions, while Paper IV will split our sample into single-giant and multiple-giant systems and investigate the differences and commonalities between these two groups. These undertakings may provide new insight into the formation and evolution of this class of planets that played a crucial role in sculpting the final architecture of our own Solar System. This research has made use of NASA's Astrophysics Data System. This work has made use of data from the European Space Agency (ESA) mission Gaia (cosmos. esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, cosmos.esa.int/web/ gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.
Finally, the authors wish to recognize and acknowledge the significant cultural role and reverence that the summit of Maunakea has long had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.
Software: All code used in this paper is available at github.com/California-Planet-Search/rvsearch and github.com/leerosenthalj/CLSII. This research makes use of GNU Parallel (Tange 2011). We made use of the following publicly available Python modules: astropy (