The California Legacy Survey I. A Catalog of 178 Planets from Precision Radial Velocity Monitoring of 719 Nearby Stars over Three Decades

We present a high-precision radial velocity (RV) survey of 719 FGKM stars, which host 164 known exoplanets and 14 newly discovered or revised exoplanets and substellar companions. This catalog updated the orbital parameters of known exoplanets and long-period candidates, some of which have decades-longer observational baselines than they did upon initial detection. The newly discovered exoplanets range from warm sub-Neptunes and super-Earths to cold gas giants. We present the catalog sample selection criteria, as well as over 100,000 radial velocity measurements, which come from the Keck-HIRES, APF-Levy, and Lick-Hamilton spectrographs. We introduce the new RV search pipeline RVSearch that we used to generate our planet catalog, and we make it available to the public as an open-source Python package. This paper is the first study in a planned series that will measure exoplanet occurrence rates and compare exoplanet populations, including studies of giant planet occurrence beyond the water ice line, and eccentricity distributions to explore giant planet formation pathways. We have made public all radial velocities and associated data that we use in this catalog.


INTRODUCTION
Expanding and characterizing the population of known exoplanets with measured masses, orbital periods, and eccentricities is crucial to painting a more complete picture of planet formation and evolution. A census of diverse exoplanets sheds light on worlds radically different than Earth, and can provide insight into how these planets, as well our own Solar System, formed. For instance, the mass, semi-major axis, and eccentricity distributions of giant planets can be used to constrain formation scenarios for these objects. Nielsen et al. (2019) and Bowler et al. (2020) used mass and eccentricity constraints from direct imaging surveys to show that planetary-mass gas giants likely form via core accretion (Pollack et al. 1996), while more massive brown dwarfs and other substellar companions likely form via gravitational instability in protoplanetary disks (Boss 1997). The present-day architectures and orbital properties of planetary systems can also be used to constrain their migration histories. Dawson & Murray-Clay (2013) used a sample of giant planets with minimum masses and orbits constrained by radial velocity (RV) observations to provide evidence that giant planets orbiting metal-rich stars are more likely to be excited to high eccentricities or migrate inward due to planet-planet interactions. Many related questions remain unanswered. What is the mass-period 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? How do small, close-in planets arrive at their final masses and system architectures? What is the relationship between these warm small planets and cold gas giants; are their formation processes related? These questions can only be answered with an expansive and rigorously constructed census of exoplanets with measured masses and wellconstrained orbits.
The community has made substantial progress on these fronts over the past few decades via targeted RV surveys. For instance, Bryan et al. (2016) surveyed 123 known giant hosts to study outer giant companions; they found that half of all giants have an outer companion, with tentatively declining frequency beyond 3 AU. Similarly, Knutson et al. (2014) found a 50% companion rate for transiting hot Jupiters using a sample of 51 stars. These two results suggest a planet formation process that favors giant multiplicity. On the small-planet front, Bryan et al. (2019) constructed an RV survey of 65 super-Earth hosts and found a giant companion rate of 39±7%. This suggests that these two populations are related in some way. Some questions have seen conflicting answers, requiring further work with a more expansive RV survey. For instance, Fernandes et al. (2019) studied planet occurrence as a function of orbital period by extracting the planetary minimum masses and periods, as well as completeness contours, from a catalog plot shown in Mayor et al. (2011), which presented a HARPS and CORALIE blind radial velocity survey of 822 stars and 155 planets over 10 yr (corresponding to a 4.6 AU circular orbit around a solar-mass star). The HARPS and CORALIE radial velocities were not published in Mayor et al. (2011), which 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. In sharp contrast, Wittenmyer et al. (2020) used their radial velocities from the Anglo-Australian Planet Search to construct a blind survey of 203 stars and 38 giant planets over 18 yr. They found that giant planet occurrence is roughly constant beyond the water ice line, out to almost 10 AU. The discrepancy between these two results needs to be resolved.
The California Planet Search team (CPS; Howard et al. 2010a) has conducted many RV surveys over the past three decades, in order to find exoplanets, measure their minimum masses, and characterize their orbits. Many of these surveys were designed explicitly for the purpose of studying planet occurrence. Therefore, they used stellar samples that were constructed without bias toward stars with known planets, or an increased likelihood of hosting planets, such as metal-rich stars (Gonzalez 1997). For instance, the Keck Planet Search (Cumming et al. 2008) used 8 yr of Keck-HIRES data collected from 585 FGKM stars to study the occurrence of gas giants with periods as long as the survey baseline, measured the mass-period distribution of giant planets out to 5 AU, and found an increase in gas giant occurrence near the water ice line. The Eta-Earth Survey (Howard et al. 2010b) used 5 yr of Keck-HIRES data collected from 166 Sun-like stars to measure the occurrence of planets with orbital periods less than 50 days, ranging from super-Earths to gas giants, and found both an abundance of planets within 10 day orbits and a mass function that increases with decreasing mass for close-in planets. The APF-50 Survey combined 5 yr of highcadence Automated Planet Finder data on a sample of 50 bright, nearby stars with 20 yr of Keck-HIRES data to constrain the mass function of super-Earths and sub-Neptunes, and discovered several planets of both varieties (Fulton et al. 2016).
We constructed an aggregate survey from these distinct RV surveys, known hereafter as the California Legacy Survey (CLS), in order to measure exoplanet occurrence, particularly for planets with long orbital periods. We selected every star in the CPS catalog that was observed as part of an occurrence survey, added 31 CPS stars that satisfied our stellar selection criteria (described below), and regularly observed these stars using the Keck and UCO-Lick observatories. The California Legacy Survey contains 103,991 RVs, and reaches observational baselines beyond three decades. We wrote an automated planet search pipeline to systematically recover all planets that are detectable in the CLS and to measure the search completeness of each star's RV time series. We can use these completeness contours to calculate exoplanet occurrence rates with respect to planetary and host-star properties (e.g. Cumming et al. 2008;Howard et al. 2010b).
In this paper, we present the CLS stellar sample and the 164 known exoplanets orbiting these stars, as well as 14 newly discovered and vetted exoplanets and substellar companions. In Section 2, we describe our methodology for stellar selection. In Section 3, we describe the RVs measured for this survey. In Section 4, we describe our methods for computing the stellar properties of our sample. In Section 5, we describe the methods by which we search for exoplanets in the RVs, confirm their planetary status, and characterize their orbits. In Section 6, we present the catalog of known exoplanets, and describe in detail each of the new exoplanet candidates. In Section 7, we discuss the significance of our catalog, and conclude with plans for future work.

STELLAR SAMPLE SELECTION
Our goal for this study was to construct a sample of RV-observed FGKM stars and their associated planets, in order to provide a stellar and planetary catalog for occurrence studies. We want a survey that is quantifiably complete in some way, such as being volume-or magnitude-limited, so that we can perform unbiased occurrence measurements. One way to do this would be to observe every HD star within our desired range of stellar parameters, with the same cadence and a thirty year baseline. Given the constraints of finite observing time and instrumental magnitude limits, this is not possible. More importantly, there is no achievable, Platonic ideal of a quantifiably complete survey. However, we can approximate one by selecting CPS-observed stars that were originally chosen without bias toward a higher-or lower-than-average likelihood of hosting planets. Multiple CPS surveys, including the Keck Planet Search and Eta-Earth Survey, performed their stellar selection with these criteria.
We began with the Keck Planet Search sample, so that we can make direct comparisons to their results. We then supplemented those 585 stars with 135 stars that were not originally included as part of that sample, but they have since been observed by the CPS team and satisfy a set of criteria intended to ensure survey quality and statistical rigor for planet occurrence measurements. We selected these criteria to ensure data quality, both of individual measurements and stellar datasets, and proper stellar selection, without bias toward known or likely planet hosts, which would skew our occurrence measurements. We included CPS-observed stars that have at least 20 total RVs and at least 10 High Resolution Echelle Spectrometer (HIRES) RVs collected after the HIRES CCD upgrade in 2004, to guarantee enough RVs for well-constrained Keplerian fits, and have an observational baseline of at least 8 yr, which is the maximum baseline of the Cumming et al. (2008) sample at the time of publication. All stars in the Keck Planet Search sample pass these criteria, since we have collected more than 10 new HIRES RVs for each of them since 2004.
In order to ensure proper stellar selection, we did not include CPS-observed stars that were chosen for surveys that deliberately selected known planet hosts, metal-rich stars, or non-main-sequence stars, since these surveys would bias planet occurrence measurements. We excluded stars that were observed as part of the "N2K" and "M2K" surveys, which targeted metal rich stars to search for gas giants (Fischer et al. 2005;Apps et al. 2010). We excluded all massive stars that were observed as part of a search for planets orbiting subgiants (Johnson et al. 2010b), since that survey used a particular observing strategy geared solely toward detecting giant planets. We excluded all young stars that were selected for CPS observing based on photometric IR excess, since such stars were selected for an increased probability of planet occurrence (Hillenbrand et al. 2015). We excluded all stars from the "Friends of Hot Jupiters" surveys, which targeted known planet hosts (Knutson et al. 2014). For the same reason, we excluded all stars that were observed as part of Kepler, K2, TrES, HAT, WASP, or KELT transiting planet surveys (Bakos et al. 2002;Alonso et al. 2004;Pollacco et al. 2006;Pepper 2007;Borucki 2016).
This selection process left us with 719 stars. Figure 1 shows the entire CLS samples as a Venn diagram, illustrating the overlap of the Cumming et al. (2008) sam-ple with the Eta-Earth (Howard et al. 2010b) and 25 pc northern hemisphere volume-limited (Hirsch et al. 2020) samples. The 25 pc sample includes 255 G and early K dwarfs with apparent V magnitudes ranging from V ≈ 3 to V ≈ 9. These stars have a median temperature of 5360 K and a median mass of 0.86 M . The median number and duration of RV observations for this sample was 71 RVs spanning 21 yr, while the minimum number and duration of observations in the sample was 20 RVs spanning 3 yr. The architects of all three of these surveys designed them for planet occurrence studies. Therefore, they did not construct these catalogs by selecting on properties known to correlate or anticorrelate with planet occurrence. There are only 31 stars in the California Legacy Survey that do not belong to any of these three surveys but do still pass of our selection criteria. This survey has no hard constraints on distance, apparent magnitude, or color, as seen in Figure  4.  (Vogt et al. 1994) has been in operation on the Keck I Telescope since 1994 and has been used to measure stellar RVs via the Doppler technique since 1996 (Cumming et al. 2008). This technique relies on measuring the Doppler shift of starlight relative to a reference spectrum of molecular iodine, which is at rest in the observatory frame (Butler et al. 1996). We consistently set up HIRES with the same wavelength format on the CCDs for each observation and followed other standard procedures of CPS Howard et al. (2010a). With the iodine technique, starlight passes through a glass cell of iodine gas heated to 50 • C, imprinting thousands of molecular absorption lines onto the stellar spectrum, which act as a wavelength reference. We also collected an iodine-free "template" spectrum for each star. This spectrum is naturally convolved with the instrumental point spread function (PSF) and is sampled at the resolving power of HIRES (R = 55,000-86,000, depending on the width of the decker used). These spectra are deconvolved using PSF measurements from spectra of featureless, rapidly rotating B stars with the iodine cell in the light path. The final, deconvolved intrinsic stellar spectra serve as ingredients in a forward-modeling procedure from which we measure relative Doppler shifts of each iodine-in spectrum of a given star (Valenti et al. 1995). We also used this process to compute uncertainties on the Doppler shifts. The uncertainty for each measurement is the standard error on the mean of the RVs for 700 segments of each spectrum (each 2Åwide) run through the Doppler pipeline. We distinguish between "pre-upgrade" RVs (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004); ∼3 m s −1 uncertainties) and "post-upgrade" RVs (2004-present; ∼1 m s −1 uncertainties). In 2004, HIRES was upgraded with a new CCD and other optical improvements. We account in the time series modeling for different RVs zero points (γ) for data from the two different eras.
The RVs reported here stem from HIRES observations with a long history. The RVs from 1996 to 2004 are based on HIRES spectra acquired by the California & Carnegie Planet Search (CCPS) collaboration and were reported in Cumming et al. (2008). CCPS continued to observe these stars, but split into two separate collaborations: CPS and the Lick-Carnegie Exoplanet Survey (LCES). This paper principally reports results from 41,804 CPS and CCPS HIRES spectra that were obtained and analyzed by our team during 1996-2020. In addition, we have included RVs computed by our pipeline for 7530 spectra of CLS stars taken by LCES during 2008-2014. These HIRES spectra were acquired with the same instrumental setup as the CPS spectra and are publicly available in the Keck Observatory Archive. Butler et al. (2017) separately published RVs based on the same HIRES observations from CCPS, CPS, and LCES for the 1996-2014 time span. The LCES and CPS Doppler pipelines diverged in ∼2007. Tal-Or et al. (2019) uncovered the 2004 zero-point offset, which we model with two independent offsets. They also claimed two second-order systematics in the LCES 2017 dataset: a long-term drift of order < 1 m s −1 , and a correlation between stellar RVs and time of night with respect to midnight. They estimated the long-term drift by averaging the zero points of three RV-quiet stars on each night, where possible. However, by our estimates, even the quietest stars exhibit 1-2 m s −1 jitter in HIRES time series. Averaging the zero points of three such stars will likely yield a scatter of 1 m s −1 across many nights. Additionally, they did not remove planet RV signals from their data before estimating the linear correlation between RV and time of night, and it is unclear how they derived the uncertainty in that correlation.

Automated Planet Finder
The APF-Levy spectrograph is a robotic telescope near the summit of Mt. Hamilton, designed to find and characterize exoplanets with high-cadence Doppler spectroscopy (Vogt et al. 2014;Radovan et al. 2014). The facility consists of a 2.4-m telescope and the Levy Spectrometer, which has been optimized for optical Doppler shift measurements. The Doppler pipeline that was developed for Keck-HIRES also extracts RV measurements from APF spectra. Most of the APF data in the California Legacy Survey was collected as part of the APF-50 Survey (Fulton 2017), the stellar sample of which was drawn entirely from the Eta Earth sample. These two surveys have slightly different selection criteria. While both surveys have a distance cut d < 25 pc and luminosity cut M V < 3, Eta-Earth cuts on apparent magnitude V < 11, whereas APF-50 has V < 7; Eta-Earth cuts on chromospheric activity logR HK < −4.7, whereas APF-50 has logR HK < −4.95; and Eta-Earth cuts on declination > −30 • , whereas APF-50 has declination > −10 • . These stricter cuts were made to ensure higher data quality for the high-cadence APF survey.

Lick-Hamilton
The Hamilton Spectrograph is a high-resolution echelle spectrometer, attached to the 3 m Shane telescope on Mt. Hamilton. Beginning in 1987, and ending in 2011 with a catastrophic iodine cell failure, the Lick Planet Search program (Fischer et al. 2014) monitored 387 bright FGKM dwarfs to search for and characterize giant exoplanets. This was one of the first surveys to produce precise RVs via Doppler spectroscopy with iodine cell calibration, and yielded RVs with precision in the range 3-10 m s −1 . The Lick Planet Search overlaps heavily with the Keck Planet Search and other CPS surveys, since these surveys drew from the same bright-star catalogs.

Activity Indices
For each HIRES and APF spectrum from which we measure radial velocities, we also measure the strength of emission in the cores of the Ca II H & K lines (Svalues) following the techniques of Isaacson & Fischer (2010) and Robertson et al. (2014). There is a small, arbitrary offset between the HIRES and APF activity indices. We adopted uniform S-value uncertainties with values of 0.002 and 0.004 for HIRES and APF respectively. We provide activity indices along with our RV measurements. Missing values are the result of sky contamination and/or low SNR.

APT Photometry
We collected long-term photometric observations of the subset of our sample that were included in the APF-50 survey (Fulton 2017), in order to search for evidence of rotation-induced stellar activity. We collected these measurements with Tennessee State University's Automated Photometric Telescopes (APTs) at Fairborn Observatory as part of a long-term program to study stellar magnetic activity cycles (Lockwood et al. 2013). Most stars have photometric datasets spanning 15 -23 yr. The APTs are equipped with photomultiplier tubes that measure the flux in the Stromgren b and y bands relative to three comparison stars. We combined the differential b and y measurements into a single (b+y)/2 "passband" then converted the differential magnitudes into a relative flux normalized to 1.0. The precision in relative flux is typically between 0.001 and 0.0015. Further details of the observing strategy and data reduction pipeline are available in Henry (1999); Eaton et al. (2003);Henry et al. (2013). We make the photometric data available as a machine-readable table.

Observational Statistics
We examined the range of observing cadences and observational baselines within the CLS sample, to determine whether stars without known planets were observed with strategies that differed significantly from those for stars with known planets. Figure 2 shows the distribution of number of observations and observational baselines for three groups of stars: the entire sample, the stars around which we detected planets, and the star around which we did not detect planets. Each of these three samples has a median baseline of 21 yr. Stars with detected planet have a median of 74 observations, compared to 35 observations for stars without detected planets and 41 observations for the entire CLS sample. A factor of two in number of observations will have a small but measurable impact on planet detectability of a given data set -and therefore on its search completeness contours.  (Huber 2017) software packages to the template Keck-HIRES spectra of our stars. Specmatch takes an optical stellar spectrum as input, and by interpolating over a grid of template spectra with known associated stellar properties, returns three spectral properties and uncertainties. For stars hotter than 4700 K, we interpolated over synthetic spectra to derive spectral parameters (Petigura 2015). For stars below this threshold, we interpolated over real spectra of cool stars with wellcharacterized stellar properties, since synthetic spectral models are unreliable below this temperature (Yee et al. 2017).
Specmatch produces metallicity, effective temperature, and surface gravity when interpolating over synthetic spectra; it produces metallicity, effective temperature, and radius when interpolating over empirical spectra. Isoclassify takes effective temperature, metallicity, and surface gravity as spectral parameter inputs, and uses isochrone models and multinest Bayesian sampling (Buchner 2016) to produce estimates and uncertainties of physical parameters, in particular stellar mass. For stars cooler than 4700 K, we passed Isoclassify a wide Gaussian input prior on surface gravity, since temperature and metallicity strongly constrain the masses of cool, main-sequence stars (Johnson et al. 2017).
Almost all stars in the California Legacy Survey have both Gaia-measured parallaxes (Gaia Collaboration et al. 2016Collaboration et al. , 2018Lindegren et al. 2018) and apparent K-band magnitudes. For stars with both of these measurements available, we pass them and their uncertainties into Isoclassify as additional inputs, since taken together, they constrain stellar luminosity and therefore place tighter constraints on stellar mass. Isoclassify also returns more precise estimates of stellar radius when provided with parallax and apparent magnitude. With the inclusion of this luminosity constraint, the median precision of our stellar mass measurements is 3.6%.
In Table 2 in Appendix B, we report stellar mass, radius, surface gravity, effective temperature, and metallicity for a subsection of the CLS sample. We make this table available for the entire sample in machine-readable format, with additional columns including V -band magnitude and Gaia parallax. Figure 3 is a visualization of these stellar properties, while Figure 4 shows individual histograms for mass, metallicity, and effective temperature, as well as for the following observational properties: parallax-inverse distance, V , and B − V .

Planet Search
We developed an iterative approach to a search for periodic signals in RV data in order to generate the CLS planet catalog. We outline this algorithm, which we de-  veloped as the open-source Python package RVSearch and have made public alongside the publication of this paper. Figure 5 is a flowchart that lays out each step of the algorithm, and Figure 6 is a visualization of an example RVSearch output, where the top two panels show the final model, and each successive row shows an iterative search for each signal in the model. First, we provide an initial model, from which the iterative search begins. This initial model contains an RV data set, and a likelihood function. The natural logarithm of the latter is defined as where i is the measurement index, v i is the ith RV measurement, γ D is the offset of the instrumental dataset from which the ith measurement is drawn, and σ 2 i is the quadrature sum of the instrumental error and the stellar jitter term of the ith measurement's instrumental data set. Here, m(t i ) is the model RV at time t i , defined as where n is a given Keplerian orbit in the model, K(t|K, P, e, ω, t c ) is the Keplerian orbit RV signature at time t given RV amplitude K, period P , eccentricity e, argument of periastron ω, and time of inferior conjunction t c ,γ is a linear trend term,γ is a quadratic trend term, and t 0 is a reference time, which we defined as the median time of observation.
We used RadVel (Fulton et al. 2018) to fit Keplerian orbits. The initial likelihood model contains either a one-planet Keplerian model with undefined orbital parameters, or a predefined model including trend/curvature terms and/or Keplerian terms associated with known orbital companions. We defaulted to performing a blind search starting with the undefined single-planet model, and we only supply a predefined model if there is evidence for a highly eccentric companion whose period is misidentified by our search algorithm. Several highly eccentric stellar binaries satisfy this criterion, as do two planets: HD 120066 b (Blunt et al. 2019), and HD 80606 b (Wittenmyer et al. 2007b).
Before beginning a blind search, RVSearch determines whether the data merits a trend with curvature, a linear trend, or no trend. It does this by fitting each of these three models to the data, then performing a goodness-of-fit test to decide which model is favored. . Stellar parameter distributions. The left column shows mass, metallicity, and effective temperature, while the right column shows parallax-inferred distance, V -band magnitude, and B − V color. Black lines are histograms of the stellar parameter median values. For the left column, colored lines are 500 histograms per panel, with parameters redrawn from normal distributions with width equal to their individual measurement uncertainties. We left these redrawn parameter histograms out for the plots in the right column because distance, magnitude, and color have uncertainties that are smaller than the chosen bin size.
We measured the Bayesian Information Criterion (BIC) for each of the three models, and computed the ∆BIC between each model. RVSearch selects the linear model if it has ∆BIC = 5 with respect to the flat model, and the quadratic model if it has ∆BIC = 5 with respect to the linear model. We did not perform this test on datasets that contain eccentric companions with orbital periods greater than the data's observational baseline, since such datasets would be better fit with a long-period Keplerian orbit than with linear and parabolic trends. The Bayesian information criterion is defined as where n obs is the number of observations, k is the number of free model parameters, and ln(L) is the loglikelihood of the model in question.
Once we provide an initial model, RVSearch defines an orbital period grid over which to search, with sampling such that the difference in frequency between adjacent grid points is 1 2πτ , where τ is the observational baseline. We chose this grid spacing in accordance with Horne & Baliunas (1986), who state that, in frequency space, a Lomb-Scargle periodogram has a minimum peak width of 1 2πτ . For each dataset, we searched for periodicity between two days and five times the observational baseline. Searching out to five times the baseline only adds a few more points to the period grid, and it allows for the possibility of recovering highly eccentric, ultra-longperiod planet candidates with best-fit orbital period.
The search algorithm then computes a goodness-of-fit periodogram by iterating through the period grid and fitting a sinusoid with a fixed period to the data. We measure goodness-of-fit as the ∆BIC at each grid point between the best-fit, n+1-planet model with the given fixed period, and the n-planet fit to the data (this is the zero-planet model for the first planet search).
After constructing a ∆BIC periodogram, the algorithm performs a linear fit to a log-scale histogram of the periodogram power values. The algorithm then extrapolates a ∆BIC detection threshold corresponding to an empirical false-alarm probability of 0.1%, meaning that, according to the power-law fit, only 0.1% of periodogram values are expected to fall beyond this threshold. This process follows the detection methodology outlined in Howard & Fulton (2016).
If a periodic signal exceeds this detection threshold, RVSearch refines the fit of the corresponding Keplerian orbit by performing a maximum a posteriori (MAP) fit with all model parameters free, including eccentricity, and records the BIC of that best-fit model. RVSearch includes two hard-bound priors, which constrain K > 0 and 0 <= e < 1. The algorithm then adds another planet to the RV model and conducts another grid search, leaving all parameters of the known Keplerian orbits free so that they might converge to a more optimal solution. In the case of the search for the second planet in a system, the goodness of fit is defined as the difference between the BIC of the best-fit one-planet model and the BIC of the two-planet model at each fixed period in the grid. RVSearch once again sets a detection threshold in the manner described above, and this iterative search continues until it returns a nondetection.
This iterative periodogram search is superior to a Lomb-Scargle residual subtraction search in two key ways. First, this process fits for the instrument-specific parameters of each dataset, stellar jitter and RV-offset, as free parameters throughout the search. Second, by leaving the known model parameters free while searching for each successive planet, we allow the solutions for the already discovered planets to reach better maxlikelihood solutions that only become evident with the inclusion of another planet in the model.
Note that our search and model comparison process is not Bayesian; we do not use priors to inform our model selection, and we do not sample posteriors, beyond a grid search in period space, until we settle upon a final model. We use the BIC as our model comparison metric because it incorporates the number of free parameters as a penalty on more complex models, which, in our case, corresponds to models with additional planets.
We make RVSearch publicly available alongside this paper via a GitHub repository. See the RVSearch website for installation and use instructions.

Search Completeness
We characterized the search completeness of each individual dataset, and of the entire survey, by running injection-recovery tests. Once RVSearch completed an iterative search of a dataset, it injected synthetic planets into the data and ran one more search iteration to determine whether it recovers these synthetic planets in that particular dataset. We ran 3000 injection tests for each star. We drew the injected planet period and M sin i from log-uniform distributions, and drew eccentricity from the beta distribution described in Kipping (2013), which was fit to a population of RV-observed planets.
We used the results of these injection tests to compute search completeness for each individual dataset, and report constant M sin i / a contours of detection probability. Figure 7 shows examples of these contours and the corresponding RVs for three different stars, all early G-type: one with 25 observations, one with 94, and one with 372. We make the 10th, 50th, and 90th percentile completeness contours for each individual star available in machine-readable format.

Model posteriors
Once RVSearch returned max-likelihood estimates of the orbital model parameters for a given dataset, we sampled the model posterior using affineinvariant sampling, implemented via emcee and RadVel (Foreman-Mackey et al. 2013;Fulton et al. 2018). We sampled using the orbital parameter basis { logP K t c √ esinω √ ecosω }. We placed uniform priors on all fitting parameters, with hard bounds such that K > 0 and 0 ≤ e < 1. We fit in logP space to efficiently sample orbits with periods longer than our observational baseline, and in √ esinω and √ ecosω to minimize bias toward higher eccentricities (Lucy & Sweeney 1971). We reported parameter estimates and uncertainties as the median and ±1σ intervals.
If a dataset is so poorly constrained by a Keplerian model that emcee's affine-invariant sampler cannot efficiently sample the posterior distribution, we instead used a rejection sampling algorithm to estimate the posterior. In these cases, we used TheJoker (Price-Whelan et al. 2017), a modified MCMC algorithm designed to sample Keplerian orbital fits to sparse radial velocity measurements. We chose a flat prior on logP , with a minimum at the observing baseline and a maximum at twenty times the observing baseline. We drew orbital eccentricity from a beta prior weighted toward zero, as modeled in Kipping (2013), in order to downweight orbits with arbitrarily high eccentricity, which can be vi-able fits to sparse or otherwise underconstraining RV data sets.

False-positive vetting
We performed a series of tests to vet each planet candidate discovered by our search pipeline. The following subsections each detail one test we perform to rule out one way in which a signal might be a false-positive. We also represent this process with a flowchart in Figure 9, and include a table of all false-positive signals recovered by RVSearch in Table 6.

Stellar activity, magnetic/long-period
Many main-sequence stars, particularly F-and Gtype, have magnetic activity cycles on timescales of several to tens of years. These fluctuations in activity can cause changes in the core depths of stellar Calcium H & K lines, which manifest as apparent RV shifts (Isaacson & Fischer 2010). To evaluate whether stellar activity may be the cause of a signal recovered by our search pipeline, we measure the linear correlation between the RV signature of that signal and a measured stellar activity metric-in our case, S-values. We computed S-values for both post-upgrade HIRES and APF data by measuring the core flux of Calcium H & K lines.
If we found a periodic signal in the S-value data that has a similar period and phase similar to one of the Keplerian terms in our RV model, we searched for correlations between our RV model and S-values. If we found one periodic signal in an RV dataset, we measured its correlation with stellar activity simply as the linear correlations between the RVs of each instrument and their associated S-values. If we found multiple periodic signals, then for each signal, we subtracted the associated RV models of all other signals from the data, and measured the correlations between these residuals and the S-values. A significant linear correlation between a signal's RV residuals and the associated S-values does not necessarily mean that this signal is caused by stellar activity, even when these signals also have the same period and phase, but we took it as sufficient evidence to remove such signals from our catalog of confirmed planets.
It is important to note that our approach to vetting our planet candidates is systematic but not exhaustive, particularly with respect to stellar activity. One might use activity metrics beyond S-values and photometry, such as Hα line modulation. Furthermore, there are more sophisticated ways to deal with activity than searching for linear correlations with RVs. For instance, one might actively model stellar activity during the search process, using a Gaussian Process . RVs and completeness contours for three datasets with similar baselines, median measurement errors, and stellar jitter. The left column plots RVs with respect to time, while the right column plots injected signals in the M sin i and a plane, where blue dots are recovered injections and red dots are not. The right column also shows detection probability contours, with 50% plotted as a solid black line. From top to bottom, we show RVs and contours for HD 44420, for which we have 24 RVs; HD 97343, for which we have 94 RVs; and HD 12051, for which we have 372 RVs.
candidate parameters and catalog selection, but require case-by-case analysis for each stellar system, as activity modeling is sometimes unwarranted or even counterproductive, e.g., for low-activity stars or confirmed planets that have periods similar to their host star's activity cycle. We chose to perform uniform, after-the-fact vetting for our catalog, and invite others to perform more sophisticated modeling for individual systems of interest.

Stellar activity, rotation/short-period
We only detected planet candidates that are lowamplitude and short-period enough to possibly be stellar rotation false positives in sustained, high-cadence datasets. Almost all CLS datasets that satisfy this criteria were collected as part of the APF-50 survey. We collected APT photometry of all APF-50 stars, which we can use to search for evidence of stellar rotation with moving-average smoothing and periodogram analysis. If we find strong evidence for rotation in APT photometry, or spectral S-value measurements, we discount planet candidates with periods close to the apparent rotation timescale or its harmonics.

Yearly alias
When we find a signal with a period of a year or an integer fraction of a year, we investigate whether it is an alias of long-period power, or a systematic that is correlated with the barycentric velocity at the time of observation or Doppler fitting parameters. We do this by recomputing the associated RVs using a different template observation. When another template observation was unavailable, we were able to take one using Keck-HIRES during collaborator observing nights. Templates taken in poor observing conditions or when barycentric velocity with respect to the observed star is high can produce systematic errors in the Doppler code. If a search of this new dataset returns a nondetection, or detection at a significantly different period, we conclude that this signal is an alias. Figure 8 shows the presence of yearly alias power in our survey, seen in a stack of the the final nondetection periodograms of all CLS stars.   Figure 9. Candidate vetting flowchart.

PLANET AND STELLAR/SUBSTELLAR COMPANION CATALOG
We present orbital solutions for the known planets, substellar companions, and stellar binaries that RVSearch has recovered in the California Legacy Survey. As mentioned in Section 5.1, where appropriate, we modeled long-period companions with linear or parabolic trends. We included in the appendix portions of the tables associated with each class of object: one for planets, one for stellar and substellar companions that are best modeled by Keplerian orbits, and one for stars with linear or parabolic RV trends. We also present 14 newly confirmed or significantly revised exoplanets and substellar companions. We list them and their orbital parameters in Table 1, and include individual notes on each system in Appendix A. Figure 10 shows all recovered planets in our survey, and distinguishes between known planets and new discoveries.

DISCUSSION
Through the use of high-cadence APF observations and long-baseline HIRES observations, we have expanded the population of known exoplanets along the current mass and semi-major axis boundary of detectability, as seen in Figure 10. We recovered 43 planets with M sin i < 30 M ⊕ , including four new discoveries within 1 AU. In a future paper in the California Legacy Survey series, we will leverage the decades-long-baseline datasets in which these planets were discovered, in order to constrain the probability that a host of a small planet also hosts an outer companion, as explored in Bryan et al. (2019) and Zhu & Wu (2018). We will also directly place a lower limit on the conditional occurrence of inner small planets given the presence of an outer gas giant.
In addition to expanding the population of small planets with measured M sin i , we discovered or revised the orbits of ten planets with orbital separations greater than 1 AU, six of them beyond 4 AU. We represent the model posteriors for the coldest of these planets in Figure 11, and show a gallery of some of their orbits in Figure 12. These discoveries include two new detections with incomplete orbits, HD 213472 b and HD 26161 b. Details are provided in Appendices A.3 and A.14. Using HIRES to extend the observational baseline of our survey by another decade will tighten our M sin i and orbital parameter constraints for these planets, and may reveal more cold companions beyond 10 AU.
In a future paper in the CLS series, we will use our sample of long-period planets and completeness contours to measure the mass-period planet occurrence distribution out to 10 AU, extending beyond the Keck Planet Search's limit of 5 AU (Cumming et al. 2008) and the 9 AU limit of Wittenmyer et al. (2020). This will provide novel constraints on planet occurrence beyond the water ice line, resolve the discrepancy between the results of Fernandes et al. (2019) and those of Wittenmyer et al. (2020), and provide new insight into planet formation across protoplanetary disks. Figure 13 is a visualization of the eccentricities of all planets in the California Legacy Survey. In future work, we will quantify the eccentricity distribution of gas giants in our sample and its dependence on planet mass and multiplicity, as well as the eccentricity distributions  . Contours (1-and 2-σ) of M sin i and semi-major axes for planets in the CLS sample whose semi-major axis posteriors extend beyond 10 AU. Contours for HD 26161 b have hard cutoffs due to sparsity below 7 MJ and 12 AU; these limits come from the data's baseline and RV increase to date. of brown dwarfs and other substellar companions, in order to clarify possible formation pathways. We will extend the wide-orbit population comparisons of Bowler et al. (2020) to our sample of planets and brown dwarfs within 20 AU of their hosts. We will also explore the eccentricity distribution of gas giants beyond 7 AU. As Figures 12 and 13 show, all planets recovered beyond 7 AU are eccentric with significance e > 2σ e . This may be a selection effect, as the median baseline of observations in our sample is 21 yr, which corresponds to a semi-major axis of 7.6 AU for a planet orbiting a solarmass star. It is possible that planets with orbital periods beyond our observational baselines are more easily detectable if they are eccentric. We can use injectionrecovery tests to determine whether there is a detection bias toward eccentric planets beyond observational baselines. If this phenomenon is not a selection effect, it might imply that most giant planets beyond 7 AU have undergone a scattering event or otherwise been excited to high eccentricity. Taken together, these studies will leverage this decades-long observational undertaking to provide new insights into planet formation and evolution. and others who contributed to the observations and analysis reported here. We acknowledge R. P. Butler and S. S. Vogt for many years of contributing to this dataset. This research has made use of the Keck Observatory Archive (KOA), which is operated by the W. M. Keck Observatory and the NASA Exoplanet Science Institute (NExScI), under contract with the National Aeronautics and Space Administration. We acknowledge RVs stemming from HIRES data in KOA with principal investigators from the LCES collaboration (S. S. Vogt, R. P. Butler, and N. Haghighipour). We gratefully acknowledge the efforts and dedication of the Keck Observatory staff for support of HIRES and remote observing. We are grateful to the time assignment committees of the Caltech, the University of California, the University of Hawaii, NASA, and NOAO for their generous allocations of observing time. Without their longterm commitment to radial velocity monitoring, these planets would likely remain unknown.
We thank Ken and Gloria Levy, who supported the construction of the Levy Spectrometer on the Automated Planet Finder, which was used heavily for this research. We thank the University of California This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www. cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www. cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

APPENDIX
In our appendices, we have included individual notes on each planet discovery reported in this paper; complete tables of recovered planets, Keplerian-resolved stellar binaries, and substellar companions in the California Legacy Survey; signals that RVSearch recovered and we determined to be false-positives; linear and parabolic RV trends; and excerpts from the stellar sample and RV dataset, which are available in their entirety in machinereadable format. A

. INDIVIDUAL DISCOVERIES AND REVISED ORBITS
A.1. HD 3765 HD 3765 is a K2 dwarf at a distance of 17.9 pc (Gaia Collaboration et al. 2018). Figure 14 shows the RVSearch results for this star. We recovered a signal with a period of 3.36 yr. Table 1 reports all planet parameters. There is significant periodicity in the S-value time series, but concentrated around a period of 12 yr. Furthermore, we find no correlation between the RVs and S-values. Figure 15 shows a Lomb-Scargle periodogram of the S-value time series. Thus, we label this signal as a confirmed planet, with M sin i = 0.173±0.014 M J and a = 2.108±0.033 AU. The magnetic activity cycle is too weak for RVSearch to recover, but is evident in the best-fit RV residuals. We used RadVel to model this activity cycle with a squared-exponential Gaussian process, and report MCMC-generated posteriors for both orbital and Gaussian process parameters in Figure 16 and Figure 17.

A.2. HD 24040
HD 24040 is a G1 dwarf at a distance of 46.7 pc. Figure 18 shows the RVSearch results for this star. It hosts a known gas giant (Wright et al. 2007;Feng et al. 2015) with a semi-major axis that we measured as a = 4.72 ± 0.18 AU, an orbital period of 9.53 ± 10 −4 yr, and a minimum mass M sin i = 4.09±0.22 M J . We have extended the observational baseline of our HIRES measurements to 21.7 yr, constrained the long-term trend and curvature of the RVs, and discovered a new exoplanet, a sub-Saturn (M sin i = 0.201 ± 0.027 M J ) on a 1.4 yr orbit (a = 1.30 ± 0.021 AU) that is consistent with circular. The S-values are uncorrelated with the the RVs of both planet signals, after removing the long-term trend. Figure 19 shows a Lomb-Scargle periodogram of the S-value time series. Table 1 reports all planet parameters.
In addition to the newly detected sub-Saturn, we further constrained the known linear trend in the RVs and found evidence for a curvature term as well.
RVSearch detected a curvature term with model preference ∆BIC> 10 over a purely linear trend. We measured the linear trend to be 0.00581 ± 0.00044 m s −1 d −1 , and the curvature to be −6.6×10 −7 ±1.2×10 −7 m s −1 d −2 , a 5.5σ detection. The trend and curvature parameters are slightly correlated in the posterior, but neither is correlated with any of the Keplerian orbital parameters in the model. Therefore, we kept the curvature term that RVSearch selected in our model. This long-term trend is low-amplitude enough that it may be caused by another planet in the system, orbiting beyond 30 AU. Gaia astrometry or another two decades of RVs may provide further constraints on this object.

A.3. HD 26161
HD 26161 is a G0 dwarf located at a distance of 50.0 pc. Figure 20 shows the RVSearch results for this star. Our RVs are consistent with a long-period, eccentric companion, and RVSearch detected this long-period signal. Due to the sparseness of the data and the fractional orbital coverage, traditional MCMC methods fail to return a well-sampled model posterior. Since the data underconstrains our model, we used TheJoker to sample the posterior, which is consistent with an extremely long-period gas giant with minimum mass M sin i = 13.5 +8.5 −3.7 M J , semi-major axis a = 20.4 +7.9 −4.9 AU, and eccentricity e = 0.82 +0.06 −0.05 . Table 1 reports current estimates of all orbital parameters, and Figure 21 shows their posterior distributions. A Keplerian model is significantly preferred over a quadratic trend, with ∆BIC> 15.
The Simbad stellar catalog designates HD 26161 as a stellar multiple. We used Gaia to identify a binary companion with similar parallax and within 60 arcseconds. This companion has an effective temperature identified from Gaia colors of 4053 K, and a projected separation of 562 AU. A stellar companion that is currently separated from its primary by more than 560 AU could not cause a change in RV of 100 m s −1 over 4 yr. This curve is far more likely caused by an inner planetary or substellar companion approaching periastron. Figure 22 shows a sample of possible orbits for HD 26161 b, drawn from our rejection sampling posteriors and projected over the next decade. We will continue to monitor HD 26161 with HIRES at moderate cadence, and have begun observing this star with APF. As we gather more data during the approach to periastron, we can tighten our constraints on the minimum mass, eccentricity, and orbital separation of HD 26161 b.

A.4. HD 66428
HD 66428 is a G8 dwarf found at a distance of 53.4 pc. Figure 23 shows the RVSearch results for this star. This system has one well-constrained cold Jupiter (Butler et al. 2006b) and an outer companion candidate first characterized in Bryan et al. (2016) as a linear trend. With four more years of HIRES data, we now see curvature in the RVs and a clear detection in RVSearch, and can place constraints on this outer candidate's orbit with a Keplerian model. The Keplerian orbit for the outer candidate is preferred to a parabolic trend with ∆BIC> 30. A maximum likelihood fit gives an orbital period of P = 36.4 yr. However, since we have only observed a partially resolved orbit so far, the orbit posterior in period space is wide and asymmetric. MCMC sampling produces P = 88 +153 −49 yr. Table 1 reports current estimates of all orbital parameters.
The model parameters are M sin i = 27 +22 −17 M J , a = 23.0 +19.0 −7.6 AU, and e = 0.31 +0.13 −0.13 . This orbital companion could be a massive gas giant or a low-mass star, if we only consider constraints from RV modeling. However, Bryan et al. (2016) used NIRC2 Adaptive-Optics images to place upper bounds on the mass and semi-major axis of an outer companion, at a time when it only presented as a linear trend in HIRES RVs. They found an upper bound of ≈100 M J on mass, not just M sin i, and an upper bound of ≈150 AU on a. We will continue to monitor this star with HIRES to further constrain the mass and orbit of HD 66428 c.

A.5. HD 68988
HD 68988 is a G0 dwarf found at a distance of 61 pc. Figure 24 shows the RVSearch results for this star. This system has one well-constrained hot Jupiter (Vogt et al. 2002) and an outer companion candidate that was first characterized in Bryan et al. (2016) as a partially resolved Keplerian orbit. With four more years of HIRES data, we can place tighter constraints on this outer candidate's orbit. A maximum likelihood fit gives an orbital period of 49.2 yr. However, since we have only observed a partially resolved orbit so far, the orbit posterior is wide and asymmetric in period space. MCMC sampling produces P = 61 +28 −20 yr. The model parameters are M sin i = 17.6 +2.4 −2.5 M J , a = 16.5 +4.8 −3.8 AU, and e = 0.53 +0.13 −0.09 . Table 1 reports all companion parameters.
RVSearch detects a third periodic signal, with P = 1900 days, that has the same period and phase as the peak period in the S-value time series. This signal also has a low RV amplitude, 6 m s −1 . Therefore, we designated this signal as a false-positive corresponding to stellar activity.
A.6. HD 95735 HD 95735 (GJ 411) is an M2 dwarf found at a distance of 2.55 pc. Figure 25 shows the RVSearch results for this star. This system has one known short-period super-Earth, with M sin i = 3.53 M ⊕ and an orbital period of 12.9 days. Our detection of this planet was driven by high-cadence APF data. This planet was first reported by Díaz et al. (2019), who also noted long-period power in their SOPHIE RV data, but they did not have a sufficiently long baseline or the activity metrics necessary to determine the origin of this power. With our HIRES post-upgrade and APF observations, we have an observational baseline of 14 yr, allowing us to confirm this long-period signal as a planet with M sin i = 24.7 ± 3.6 M ⊕ and an orbital period P = 8.46 yr. Table 1 reports all planet parameters. Since GJ 411 is a cool M dwarf, the Lick-Hamilton and HIRES pre-upgrade data are not reliable, because those detectors are not sufficiently high-resolution to capture a cool M dwarf's dense spectral lines (Fischer et al. 2014).
There is a long-period trend in the HIRES S-value time series, with significant power at and beyond 25 yr, but no significant power near the orbital period of the outer candidate. Therefore, we included this candidate in our catalog as a new planet candidate, to be verified and constrained with several more years of HIRES observations.
RVSearch also recovered a highly eccentric, 216 day signal, but this signal correlates with APF systematics. Therefore, we labeled it as a false positive. This systematic remained when we applied RVSearch only to the HIRES post-upgrade and APF data and left out the problematic pre-upgrade and Lick data. There is evidence for an activity cycle longer than 10,000 days, but no significant power near the period of our 3,000-day planet candidate.

A.7. HD 107148
HD 107148 is a G5 dwarf at a distance of 49.5 pc. Figure 27 shows the RVSearch results for this star. Butler et al. (2006b) reported a planet with a period of 44 days. They reported periodicity at 77 days, but determined that this was an alias of the 44 day signal. The 77 day signal is significantly stronger in our likelihood periodogram, as seen in Figure 27, and better fits the data than a 44 day Keplerian by a significant ∆BICṪhis constitutes strong evidence that the true period of this planet is 77 days. We report new orbital parameters for this planet in Table 3.
We also recovered a signal with a period of 18.3 days. There is significant periodicity in the S-value time series, a periodogram of which is shown in Figure 28. However, it is concentrated around a period of 6 yr, and there is no significant power near 18.3 days. Furthermore, we find no correlation between the RVs and S-values. Thus, we report this signal as a confirmed planet, with M sin i = 19.9 ± 3.1 M ⊕ and a = 0.1406 ± 0.0018 AU.

A.8. HD 136925
HD 136925 is a G0 dwarf, found at a distance of 47.9 pc. RVSearch detected two periodic signals in this dataset, as seen in Figure 29, at 311 days and 12.4 yr. This dataset is currently sparse, with two gaps of several years in the post-upgrade HIRES data, but there is clear long-period variation in the RVs. Keplerian modeling predicts M sin i = 0.84 M J for the giant planet.
The S-value periodogram seen in Figure 30 shows no significant power beyond 1000 days, suggesting that the long-period HD 136925 b is a real planet. There is broad power around 300 days, overlapping with the period of the inner signal. It is unclear whether this periodicity is caused by real stellar variability or is a product of sparse data. Table 1 reports current estimates of all planet parameters. We need more data in order to clarify our model, and determine whether the inner signal is caused by a planet or a product of stellar activity and sparse data. Therefore, we designated HD 136925 b as a planet, and the inner signal as a probable false positive, to be clarified with continued HIRES observing. A.9. HD 141004 HD 141004 is a G0 dwarf found at a distance of 11.8 pc. Figure 31 shows the RVSearch results for this star. Roy et al. (2021, in preparation) discovered a sub-Neptune at an orbital period of 15.5 days, with M sin i = 13.9 ± 1.5 M ⊕ , and will report on the analysis of this system in greater detail. Table 1 reports current estimates of all planet parameters.
A.10. HD 145675 HD 145675 (14 Her) is a K0 dwarf found at a distance of 17.9 pc. Figure 32 shows the RVSearch results for this star. This system has one known cold gas giant, with M sin i= 5.10 M J and an orbital period of 4.84 yr, which was first reported in Butler et al. (2003). Wittenmyer et al. (2007a) conducted further analysis with a longer observational baseline of twelve years, and noted a long-period trend. Wright et al. (2007) used additional RV curvature constraints to show that this trend must correspond to a companion with P > 12 yr and M sin i > 5 M J . The observational baseline has since increased from 12 yr to 22, and regular observations with HIRES and APF allow us to place further constraints on this long-period companion. We find M sin i = 5.8 +1.  Table 1 reports all planet parameters. Figure 33 shows a Lomb-Scargle periodogram of the HIRES S-value time series. There is strong periodicity in the HIRES S-value time series, peaking around 10 yr, but no significant power near the supposed orbital period of the long-period candidate. These S-values strongly correlate with a third Keplerian signal picked up by our search, also with a period of 10 yr, as seen in the Figure 34, therefore we designate this signal as stellar activity.
There is a potential complication owed to a stellar binary candidate. Roberts et al. (2011) conducted a direct-imaging survey of known exoplanet hosts and reported a candidate stellar companion to 14 Her, with a differential magnitude of 10.9 ± 1.0, an angular separation of 4.3", and a minimum orbital separation of 78 AU. This is a single-epoch detection, and therefore could be only a visual binary. Additionally, Rodigas et al. (2011) conducted a deep direct imaging study of 14 Her, to constrain the mass and orbital parameters of 14 Her c, which, at the time, presented only as a parabolic trend in RV data. They used the Clio-2 photometer on the MMT, which has a 9" x 30" field of view; the authors only looked at imaging data within 2", to filter out background stars. Although this deep imaging study did not mention any stellar companion, the candidate reported by Roberts et al. (2011) falls outside of their considered imaging data, which corresponds to a minimum separation of 112.8 AU. Wittrock et al. (2017) also found a null binary detection, using the Differential Speckle Survey Instrument (DSSI) at the Gemini-North Observatory. A 6 Jupiter mass object would not have been detected by the above surveys, as they were designed only to rule out stellar companions and therefore used shorter imaging exposures that would miss planetary-mass companions.
Additionally, we used Gaia DR2 to search for bound stellar companions within 10", and found no such companions. We conclude that 14 Her does not have a bound stellar companion. Therefore, we designated 14 Her c as an eccentric, long-period planet. We will continue to monitor this star with Keck/HIRES and APF, to further constrain the orbit of this planet. In the HIRES and APF data, we measured > 3σ correlations for the third Keplerian signal. The APF and HIRES linear correlations are within 3σ of each other, implying that this signal is caused by stellar activity. We find correlations between the residuals and S-values for the second signal as well, but they are significantly different for HIRES and APF. Since the period of this signal is much greater than the APF baseline of this star, we discount this second correlation as caused by the limited baseline of the data with respect to the signal.

A.11. HD 156668
HD 156668 is a K3 dwarf found at a distance of 24.4 pc. Figure 35 shows the RVSearch results for this star. This system has one known short-period super-Earth, with M sin i= 4.15 M ⊕ and an orbital period of 4.64 days. This planet was first reported by Howard et al. (2011), who also noted a long-period (P ≈ 2.3 yr) signal with insufficient RV observations or additional data for confirmation as a planet. The observational baseline has since increased from five years to fourteen, allowing us to confirm this long-period signal as a planet with M sin i = 0.167 M J and an orbital period P = 2.22 yr.
There is a strong periodicity in the HIRES S-value time series, peaking around 10 yr, but no significant power near the orbital period of the long-period candidate. If we do not model this activity, a one-year alias signal appears in the periodogram search (Fig. 35). The data do not sufficiently constrain a Keplerian fit with a 10 yr period, but we find that a linear trend models the activity well enough to remove the one-year alias from the search. We opt to include this linear trend, which we treat as a nuisance parameter.

A.12. HD 164922
HD 164922 is a G9 V dwarf located at a distance of 22.1 pc. Figure 36 shows the RVSearch results for this star. It hosts two known planets: a 0.3 M J planet with an orbital period of 1207 days (Butler et al. 2006b), and a super-Earth with M sin i = 14.3 M ⊕ and an orbital period of 75.8 days. This super-Earth was reported by Fulton et al. (2016), who also reported residual power around 41.7 days but did not find it significant enough to merit candidate status. With approximately two more years of HIRES and APF data, we identified the 41.7 day signal as a strong planet candidate and confirmed the 12.5 day planet reported in Benatti et al. (2020). Both planets are of sub-Neptune mass and have eccentricity posteriors that are consistent with circular orbits. The 41.7 day planet has M sin i = 10.7 ± 1.0 M ⊕ and a semi-major axis a = 0.2294 ± 0.0031 AU. The 12.5 day planet has M sin i = 4.63 ± 0.70 M ⊕ and a semi-major axis a = 0.1024 ± 0.0014 AU. Table 1 reports all planet parameters.
To validate these candidates, we searched for periodicity in both S-value activity metrics and APT photometry. We found no evidence for stellar rotation in Svalues, but estimated a stellar rotation period of 62.1 days from our APT photometry. Figure 37 shows periodograms and a phase-folded curve from this APT analysis, and Figure 38 shows equivalent analysis for HIRES S-values. The 1 yr alias of 62.1 days is 75.8 days, but the 75.8 day planet detection is high-amplitude and clean, without an additional peak near 62 days in any of the RVSearch periodograms. Therefore, within the limits of our activity metrics and vetting process, we ruled out stellar rotation as a cause of the 41.7 day signal. Benatti et al. (2020) used multiple HARPS-N spectral activity indicators to estimate a stellar rotation period of 41.6 days, and they note that this rotation period is to be expected from empirical activity-rotation relationships. Therefore, they determined that the strong 42 day signal present in their HARPS RVs is caused by rotation. However, we find no evidence of significant 42 day periodicity in our analysis of spectral activity indicators or APT photometry, as seen in Figures 37 and 38, and both datasets reflect significant periodicity near 60 days. Since our RV detection of this planet candidate is clean and does not conflict with our activity analysis, we chose to include this signal in our catalog as a planet candidate, to be confirmed or refuted by independent analysis.

A.13. HD 168009
HD 168009 is a G1 dwarf found at a distance of 23.3 pc. Figure 39 shows the RVSearch results for this star. Roy et al. (2021, in preparation) discovered a super-Earth candidate at an orbital period of 15.5 days, with M sin i = 10.3 ± 1.1 M ⊕ , and will report on the analysis of this candidate in greater detail. Table 1 reports current estimates of all planet parameters.
RVSearch also recovered a highly eccentric 1 yr signal, but this signal correlates with APF systematics. Therefore, we labeled it as a false positive. Roy et al. (2021, in preparation) will model these systematics in greater detail.
A.14. HD 213472 HD 213472 is a G5 dwarf located at a distance of 64.6 pc. Figure 40 shows the RVSearch results for this star. There is an approximately eleven-year gap in RV observations of this star. The first post-upgrade HIRES observation was measured in 2005, shortly after the last pre-upgrade observation, and the second post-upgrade observation was measured in 2016. The 40 m s −1 difference between these two observations prompted the CPS team to begin observing HD 213472 regularly. Together with observations since 2016, and the thirteen pre-upgrade HIRES measurements, the data are consistent with a long-period, eccentric, planetary companion. Our periodogram search detects such a long-period signal. Due to the sparseness of the data, traditional MCMC methods fail to return a well-sampled model posterior. We used the rejection sampling algorithm TheJoker (Price-Whelan et al. 2017) to estimate the posterior, and found it to be unimodal. This model is consistent with a very long-period gas giant, with M sin i = 3.48 +1.10 −0.59 M J orbital period P = 46 +33 −13 yr, semi-major axis a = 13.0 +5.7 −2.6 AU, and eccentricity of e = 0.53 +0.12 −0.09 . Table 1 reports all planet parameters. Figure 41 shows the orbital parameter posteriors generated by TheJoker.
To investigate the possibility of a stellar or substellar companion, we compared this Keplerian model to a simple linear trend by computing the ∆BIC between the two max-likelihood models. The Keplerian model is significantly preferred with ∆BIC = 23.7. Additionally, we used Gaia to search for bound companions within 10", and found no such companions. Therefore, we inferred that HD 213472 b is either a planet or low-mass substellar companion, and not a wide-orbit stellar companion. Figure 42 shows a sample of possible orbits for HD 213472 b, drawn from our rejection sampling posteriors and projected over the next decade. More HIRES observations will further constrain this object's mass and orbital parameters.        We record all RVSearch-detected false positives in Table 7. The 'cause' column denotes why a signal was labeled as a false positive. 'A' refers to a long-period magnetic activity cycle, 'R' refers to stellar rotation, and 'N' refers to an annual and/or instrumental systematic. Long-period instrumental systematics are occasionally caused by offsets between dewars in the Lick data. Sev-eral of these false positives correspond to reported planets in the literature, or to stars that have been discussed extensively in the literature. We elaborate on each of these cases in the subsections below.  Table 7. False Positives (Continued) Figure 43. RVSearch summary plot for HD 115617. See Figure 6 for plot description. Note the nearly equivalent-height peaks at 1/3 and 1/4 year in panel h, corresponding to the 124 day reported planet. Panel j shows that there is residual power at 1 year after subtracting the 122 day signal, suggesting the presence of yearly systematic noise in the data.
ments as seen in Figure 52, and determined that 52 days is the likely stellar rotation period of HD 34445. This places our weak detection of the 49 day claimed planet candidate under suspicion, and we have labeled it as a false positive in our catalog. There is also evidence of semiannual HIRES systematics, as seen in Figure 54, which shows the correlation between HIRES RVs minus the giant planet signature and PSF parameters, and in Figure 53, which shows periodograms of each PSF parameter time series. Multiple PSF parameters correlate (|R| > 0.15) with the RV residuals, and multiple parameters show periodicity around one-third and one-fourth of a year. The two claimed planets at 118 and 215 days are close to one-third and one-half of a year, respectively, and show weak and equal-strength signatures in their RVSearch periodograms, as seen in Figure 50. Therefore, we have labeled these signals as false positives in our catalog.  Figure 46. HIRES post-upgrade RV and S-value activity timeseries for HD 154345. Note that the two datasets share minima and appear to be in phase when post-upgrade observations began, but have drifted completely out of phase over the following 23 years. Figure 47. RVSearch summary plot for HD 154345; see Figure 6 for description. RVSearch first recovered a strong signal at 9 years, but then recovered additional power at a similar period due to stellar activity. The final orbit fit switched the two models, so that panels e) and d) show the planetary signal, while panels c) and f) show the stellar activity signal.  Note that these two datasets share a negative long-term trend, which we believe accounts for the claimed 5,700-day planet in the system.