# The TESS-Keck Survey. XVI. Mass Measurements for 12 Planets in Eight
Systems111Based on observations obtained at the W. M. Keck Observatory, which
is operated jointly by the University of California and the California
Institute of Technology.

Joseph M. Akana Murphy NSF Graduate Research Fellow  Department of Astronomy
and Astrophysics, University of California, Santa Cruz, CA 95064, USA  Natalie
M. Batalha Department of Astronomy and Astrophysics, University of California,
Santa Cruz, CA 95064, USA  Nicholas Scarsdale Department of Astronomy and
Astrophysics, University of California, Santa Cruz, CA 95064, USA  Howard
Isaacson Department of Astronomy, University of California, Berkeley, CA
94720, USA  Centre for Astrophysics, University of Southern Queensland,
Toowoomba, QLD, Australia  David R. Ciardi NASA Exoplanet Science Institute-
Caltech/IPAC, Pasadena, CA 91125, USA  Erica J. Gonzales NSF Graduate Research
Fellow  Department of Astronomy and Astrophysics, University of California,
Santa Cruz, CA 95064, USA  Steven Giacalone Department of Astronomy,
University of California, Berkeley, CA 94720, USA  Joseph D. Twicken NASA Ames
Research Center, Moffett Field, CA 94035, USA  SETI Institute, Mountain View,
CA 94043, USA  Anne Dattilo Department of Astronomy and Astrophysics,
University of California, Santa Cruz, CA 95064, USA  Tara Fetherolf UC
Chancellor’s Fellow  Department of Earth and Planetary Sciences, University of
California, Riverside, CA 92521, USA  Ryan A. Rubenzahl NSF Graduate Research
Fellow  Department of Astronomy, California Institute of Technology, Pasadena,
CA 91125, USA  Ian J. M. Crossfield  Department of Physics and Astronomy,
University of Kansas, Lawrence, KS 66045, USA  Courtney D. Dressing Department
of Astronomy, University of California, Berkeley, CA 94720, USA  Benjamin
Fulton NASA Exoplanet Science Institute/Caltech-IPAC, Pasadena, CA 91125, USA
Andrew W. Howard Department of Astronomy, California Institute of Technology,
Pasadena, CA 91125, USA  Daniel Huber Institute for Astronomy, University of
Hawai ‘ i, Honolulu, HI 96822, USA  Stephen R. Kane Department of Earth and
Planetary Sciences, University of California, Riverside, CA 92521, USA  Erik
A. Petigura Department of Physics and Astronomy, University of California, Los
Angeles, CA 90095, USA  Paul Robertson Department of Physics and Astronomy,
University of California, Irvine, CA 92697, USA  Arpita Roy Space Telescope
Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  Department
of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA
Lauren M. Weiss Department of Physics and Astronomy, University of Notre Dame,
Notre Dame, IN 46556, USA  Corey Beard Department of Physics and Astronomy,
University of California, Irvine, CA 92697, USA  Ashley Chontos Henry Norris
Russell Fellow  Department of Astrophysical Sciences, Princeton University,
Princeton, NJ 08540, USA  Institute for Astronomy, University of Hawai ‘ i,
Honolulu, HI 96822, USA  Fei Dai NASA Sagan Fellow  Division of Geological and
Planetary Sciences, California Institute of Technology, Pasadena, CA, 91125,
USA  Department of Astronomy, California Institute of Technology, Pasadena, CA
91125, USA  Malena Rice Heising-Simons 51 Pegasi b Postdoctoral Fellow
Department of Physics, Massachusetts Institute of Technology, Cambridge, MA
02139, USA  Kavli Institute for Astrophysics and Space Research, Massachusetts
Institute of Technology, Cambridge, MA 02139, USA  Department of Astronomy,
Yale University, New Haven, CT 06511, USA  Judah Van Zandt Department of
Physics and Astronomy, University of California, Los Angeles, CA 90095, USA
Jack Lubin Department of Physics and Astronomy, University of California,
Irvine, CA 92697, USA  Sarah Blunt Department of Astronomy, California
Institute of Technology, Pasadena, CA 91125, USA  Alex S. Polanski Department
of Physics and Astronomy, University of Kansas, Lawrence, KS 66045, USA  Aida
Behmard NSF Graduate Research Fellow  Division of Geological and Planetary
Science, California Institute of Technology, Pasadena, CA 91125, USA  Paul A.
Dalba Heising-Simons 51 Pegasi b Postdoctoral Fellow  Department of Astronomy
and Astrophysics, University of California, Santa Cruz, CA 95064, USA  SETI
Institute, Carl Sagan Center, Mountain View, CA 94043, USA  Michelle L. Hill
Department of Earth and Planetary Sciences, University of California,
Riverside, CA 92521, USA  Lee J. Rosenthal Department of Astronomy, California
Institute of Technology, Pasadena, CA 91125, USA  Casey L. Brinkman Institute
for Astronomy, University of Hawai ‘ i, Honolulu, HI 96822, USA  Andrew W.
Mayo Department of Astronomy, University of California, Berkeley, CA 94720,
USA  Emma V. Turtelboom Department of Astronomy, University of California,
Berkeley CA 94720, USA  Isabel Angelo Department of Physics and Astronomy,
University of California, Los Angeles, CA 90095, USA  Mani L. Bhaumik
Institute for Theoretical Physics, University of California, Los Angeles, CA
90095, USA  Teo Močnik Gemini Observatory/NSF’s NOIRLab, Hilo, HI 96720, USA
Mason G. MacDougall Department of Physics and Astronomy, University of
California, Los Angeles, CA 90095, USA  Daria Pidhorodetska Department of
Earth and Planetary Sciences, University of California, Riverside, CA 92521,
USA  Dakotah Tyler Department of Physics and Astronomy, University of
California, Los Angeles, CA 90095, USA  Molly R. Kosiarek NSF Graduate
Research Fellow  Department of Astronomy and Astrophysics, University of
California, Santa Cruz, CA 95064, USA  Rae Holcomb Department of Physics and
Astronomy, University of California, Irvine, CA 92697, USA  Emma M. Louden
Department of Astronomy, Yale University, New Haven, CT 06511, USA  Lea A.
Hirsch University of Toronto, Mississauga, ON L5L 1C6, Canada  Jay Anderson
Space Telescope Science Institute, Baltimore, MD 21218, USA  Jeff A. Valenti
Space Telescope Science Institute, Baltimore, MD 21218, USA

###### Abstract

With _JWST_ ’s successful deployment and unexpectedly high fuel reserves,
measuring the masses of sub-Neptunes transiting bright, nearby stars will soon
become the bottleneck for characterizing the atmospheres of small exoplanets
via transmission spectroscopy. Using a carefully curated target list and more
than two years’ worth of APF-Levy and Keck-HIRES Doppler monitoring, the TESS-
Keck Survey is working toward alleviating this pressure. Here we present mass
measurements for 11 transiting planets in eight systems that are particularly
suited to atmospheric follow-up with _JWST_. We also report the discovery and
confirmation of a temperate super-Jovian-mass planet on a moderately eccentric
orbit. The sample of eight host stars, which includes one subgiant, spans
early-K to late-F spectral types ($T_{\mathrm{eff}}$ $=$ 5200–6200 K). We
homogeneously derive planet parameters using a joint photometry and radial
velocity modeling framework, discuss the planets’ possible bulk compositions,
and comment on their prospects for atmospheric characterization.

Exoplanets (498), radial velocity (1332)

††facilities: APF, Hale (PHARO), Keck:I (HIRES), Keck:II (NIRC2),
TESS.††software: astropy (Astropy Collaboration et al., 2013, 2018), celerite2
(Foreman-Mackey, 2018a), exoplanet (Foreman-Mackey et al., 2020), isoclassify
(Huber et al., 2017; Berger et al., 2020), matplotlib (Hunter, 2007), numpy
(Harris et al., 2020), pandas (pandas development team, 2020), pymc3
(Salvatier et al., 2016b), Python 3 (Van Rossum & Drake, 2009), RadVel (Fulton
et al., 2018), scikit-learn (Pedregosa et al., 2011), scipy (Virtanen et al.,
2020), SpecMatch-Emp (Yee et al., 2017), SpecMatch-Syn (Petigura et al.,
2017), starry (Luger et al., 2019), tessla (Akana Murphy, 2023), theano
(Theano Development Team, 2016).

ʻ ‘

##  1 Introduction

The _Kepler_ spacecraft (Borucki et al., 2010) taught us about the Milky Way
Galaxy’s intrinsic planet radius distribution for planets interior to 1 AU
(Howard et al., 2012; Batalha et al., 2013; Fressin et al., 2013; Petigura et
al., 2013; Fulton et al., 2017; Bryson et al., 2021). However, the long-stare,
single-field nature of the survey precluded ground-based Doppler mass
measurements for all but the brightest host stars. NASA’s TESS mission (Ricker
et al., 2014), on the other hand, with its detections of transiting exoplanets
orbiting bright, nearby stars across the full sky, is allowing us to better
understand the observed exoplanet mass and radius distribution. Of particular
interest to theories of planet formation and evolution are the masses of sub-
Neptunes, whose seeming diversity in bulk density has challenged our post-
_Kepler_ interpretations of small planet formation and evolution (e.g., Luque
& Pallé 2022). Thus far, TESS is responsible for discovering nearly 100
planets smaller than 4 $R_{\oplus}$ that also have robust mass
measurements.222Data accessed via the NASA Exoplanet Archive on 2023-Mar-28.
For planets with better than 50% and 15% fractional measurement precision in
mass and radius, respectively.

Measurements of atmospheric properties are key to understanding the interior
composition of sub-Neptunes (Rogers & Seager, 2010) as they lie at the
confluence of theoretical isocomposition curves in the mass-radius plane
(Valencia et al., 2007; Adams et al., 2008; Zeng et al., 2019; Otegi et al.,
2020). Constraints of atmospheric composition can also inform theories of
their formation and evolution histories (e.g., Madhusudhan 2019 and references
therein; Kite et al. 2020). For these reasons, transit spectra of sub-Neptunes
are extremely valuable. However, precise knowledge of the planets’ surface
gravities is required in order to interpret these data due to the degeneracy
between surface gravity and atmospheric mean molecular weight (Batalha et al.,
2019). Space-based transit photometry delivers precise planet radii, making
planet mass the dominant source of uncertainty in the surface gravity
calculation. Therefore, precise mass measurements and the substantial
investments of ground-based resources that they require remain the critical
first step in the effort to understand the physical drivers of sub-Neptune
diversity.

###  1.1 The TESS-Keck Survey: Planet atmospheres

The TESS-Keck Survey (TKS; Chontos et al. 2022), a multi-semester Doppler
monitoring campaign of promising TESS planet candidates with the Keck-HIRES
and APF-Levy spectrographs, is working to provide the precise planet mass
measurements required by future efforts in atmospheric characterization, among
other investigations (e.g., Scarsdale et al. 2021; Lubin et al. 2022). TKS
science falls along four main axes: (1) planet bulk composition, (2) system
architectures and dynamics, (3) planet atmospheres, and (4) evolved systems.
The systems presented in this work were all observed as members of science
case three (SC3), planet atmospheres.

The goal of the TKS SC3 program is to measure precise masses for transiting
planets in TESS systems that are particularly amenable to atmospheric follow-
up. The SC3 target list was constructed using input from a quantitative
selection function in addition to hand-tuning based on results from the TESS
Follow-up Observing Program (TFOP; e.g., a target would be dropped despite a
favorable selection function value if reconnaissance spectroscopy revealed the
system to be an eclipsing binary).333TFOP contributions to the target
selection process are acknowledged in Chontos et al. (2022). This work does
not make use of proprietary TFOP information beyond what was acknowledged by
Chontos et al. (2022). Details of the target selection procedure can be found
in Scarsdale et al. (2021) and Chontos et al. (2022). The latter contains the
complete TKS target list.

In short, the quantitative selection function we used to identify potential
SC3 targets strikes a balance between favorable prospects for atmospheric
characterization and Doppler observing cost. The function is the ratio of a
planet candidate’s expected transmission spectroscopy metric (TSM; Kempton et
al. 2018), a _JWST_ signal-to-noise ratio (S/N) proxy, and the estimated Keck-
HIRES exposure time required to achieve a 5$\sigma$ mass measurement. TSM is
defined as

$\text{TSM}=(\text{scale
factor})\times\frac{R_{\mathrm{p}}^{3}T_{\mathrm{eq}}}{M_{\mathrm{p}}R_{*}^{2}}\times
10^{-J/5},$ (1)

where “scale factor” is a normalization constant that depends on planet
radius,

$\text{scale factor}=\begin{cases}0.19&\text{for }R_{\mathrm{p}}<1.5\text{
}R_{\mathrm{\oplus}}\\\ 1.26&\text{for }1.5<R_{\mathrm{p}}<2.75\text{
}R_{\mathrm{\oplus}}\\\ 1.28&\text{for }2.75<R_{\mathrm{p}}<4\text{
}R_{\mathrm{\oplus}}\\\ 1.15&\text{for }4<R_{\mathrm{p}}<10\text{
}R_{\mathrm{\oplus}}.\\\ \end{cases}$

$R_{\mathrm{p}}$ and $M_{\mathrm{p}}$ are in Earth units, $T_{\mathrm{eq}}$ is
in Kelvin (and assumes zero Bond albedo and full day-night heat
redistribution), $R_{*}$ is in solar units, and $J$ is the host star’s
apparent magnitude in $J$-band. The scale factor values account for all unit
conversions. Since TSM depends on planet mass, we used the mass-radius
relation from Chen & Kipping (2017) to translate planet radius values from the
TESS object of interest (TOI) catalog (Guerrero et al., 2021) into preliminary
mass estimates. The Keck-HIRES exposure time required to achieve a 5$\sigma$
mass measurement was estimated using the methods in Plavchan et al. (2015). A
high value for this ratio indicates a more favorable target.

In order to encourage a sample of planets that was spread evenly over
parameter space, TOIs were divided into bins in stellar $T_{\mathrm{eff}}$,
planet radius, and planet instellation flux. Selection function values for
planets in the same $T_{\mathrm{eff}}$-$R_{\mathrm{p}}$-$S_{\mathrm{p}}$ bin
were compared against one another. The top five highest ranking planets in
each bin were then considered as candidates for the final TKS SC3 target list.
Though our binning technique attempted to select a sample that spanned a wide
range of host star $T_{\mathrm{eff}}$, since Keck-HIRES is not optimized for
observing cooler stars,444Keck-HIRES measures stellar radial velocities using
a warm cell of molecular iodine (Butler et al., 1996), which imprints
absorption lines on the stellar spectrum between $\sim$5000–6000 Å. Therefore,
Keck-HIRES is generally less efficient at measuring the radial velocities of M
dwarfs compared to G dwarfs, for example. our final SC3 target list is
comprised primarily of planets orbiting G dwarfs. Our target list also focuses
on sub-Neptunes since they offer reasonable expected Doppler observing costs
(compared to super-Earths) but are still not giant planets, for which the
literature already contains numerous atmospheric measurements.

At the start of the survey in 2019B, we identified 20 TESS systems with at
least one high-value planet candidate for atmospheric characterization
according to our sample selection procedure. After more than two years of
Doppler monitoring, the majority of these systems are now either already
published, e.g., HD 63935 (Scarsdale et al., 2021), HD 191939 (Lubin et al.,
2022), or the subject of publications in preparation by TKS collaborators. The
eight systems presented in this work constitute the remaining systems of the
TKS SC3 target list.

###  1.2 Targets in this work

In order of increasing TOI number, the systems presented in this work are:

  * •

HIP 8152 (TOI-266): a G dwarf hosting two sub-Neptunes.

  * •

HD 42813 (TOI-469): an early-K dwarf hosting one sub-Neptune.

  * •

HD 25463 (TOI-554): a late-F dwarf hosting a sub-Neptune and a super-Earth on
opposite sides of the radius valley (Fulton et al., 2017; Van Eylen et al.,
2018).

  * •

TOI-669: a G dwarf hosting one sub-Neptune.

  * •

HD 135694 (TOI-1247): an early-K dwarf hosting one sub-Neptune.

  * •

HIP 9618 (TOI-1471): a G dwarf hosting two sub-Neptunes, each with $P>20$ d.
The system is also host to a massive, distant companion as seen by a linear
trend in the radial velocities. The nature of the companion is uncertain.

  * •

HD 6061 (TOI-1473): an early-G dwarf hosting one sub-Neptune. The host star
also appears to be gravitationally bound to a mid-M dwarf companion (TIC
600433892). The two stars have a sky-projected separation of about 200 AU.

  * •

TOI-1736: a slightly evolved G star hosting one sub-Neptune and one
nontransiting, super-Jovian-mass planet on a moderately eccentric orbit. The
system is also host to a massive, distant companion as seen by a linear trend
in the radial velocities. The nature of the companion is uncertain.

Why are these systems attractive targets for atmospheric observations? While
not every planet presented here has an extraordinarily high TSM value—Kempton
et al. (2018) suggest that “good” targets for atmospheric characterization
have TSM $>50$, which is not true for three of the 11 transiting planets in
this work—Batalha et al. (2023) make clear that the most informative samples
for inferring population-level characteristics are not necessarily composed of
the best individual targets for atmospheric characterization. Furthermore,
Batalha et al. (2023) note that planets are often chosen for Doppler and
subsequent atmospheric follow-up because they are extreme in some way. This
novelty bias systematically disfavors planets which are in fact the Galaxy’s
most common products of planet formation. To this end, six of the 11
transiting planets presented here land on the mode of the sub-Neptune mass-
radius distribution. As noted above, the majority of the planets in the TKS
SC3 program orbit G dwarfs, stellar hosts which are currently underrepresented
in the set of atmospheric targets for _JWST_ , as much focus remains on small
planets orbiting cool stars.555While planets transiting cool stars are
typically more efficient targets for transmission spectroscopy (owing in part
to the larger planet-star radius ratio), the radius distribution of planets
orbiting M dwarfs is distinct from that of FGK dwarfs (Dressing & Charbonneau,
2013). This implies differences in the dominant channel(s) of planet formation
and evolution as a function of stellar properties, and, consequently, that the
atmospheric characteristics of planets around M dwarfs may not be
representative of planets with Sun-like hosts. Finally, this work presents a
large sample of planets with homogeneously derived physical properties,
mitigating the effects of potential systematic biases from the data analysis.
We discuss the planets’ prospects for atmospheric characterization further in
§13.

The paper is organized as follows: We summarize the TESS 2-min cadence
observations in §2. We present high resolution imaging of the host stars in §3
and describe our stellar characterization in §4. We discuss our Doppler
observations and data reduction in §5. We discuss our light curve inspection,
cleaning, and initial transit modeling in §6. We search for radial velocity
trends and nontransiting companions in §7. We examine stellar activity in §8.
In §9 we describe our joint photometry, radial velocity, and stellar activity
modeling framework. We present the results of this modeling in §10. In §12 we
discuss possible bulk compositions for the planets and place them in the mass-
radius diagram. In §13 we discuss the planets’ prospects for atmospheric
characterization. We conclude in §14. We note that the times of observations
labeled in Barycentric Julian Date (BJD) or Barycentric TESS Julian Date
(BTJD; $\mathrm{BTJD}=\mathrm{BJD}-2457000$; i.e., the TESS and Doppler
observations) were measured using the using the Barycentric Dynamical Time
standard (TDB; e.g., Eastman et al., 2010).

##  2 TESS photometry

Table 1: Summary of 2-min cadence TESS observations

System | Sectors | Observing start/end
---|---|---
|  | (UT)
HIP 8152 | 3, 30 | 2018-Sep-20/2020-Oct-21
HD 42813 | 6, 33 | 2018-Dec-11/2021-Jan-13
HD 25463 | 5, 32, 43, 44 | 2018-Nov-15/2021-Nov-06
TOI-669 | 9, 35 | 2019-Feb-28/2021-Mar-07
HD 135694 | (14 total) | 2019-Jul-18/2023-Jan-18
HIP 9618 | 17, 42, 43 | 2019-Oct-07/2021-Oct-12
HD 6061 | 17, 57 | 2019-Oct-07/2022-Oct-29
TOI-1736 | (6 total) | 2019-Nov-02/2022-Dec-23

Note. — According to data available on MAST as of 2023-Mar-07. Systems are
listed in order of increasing TOI number, starting with HIP 8152 (TOI-266).
The start and end dates of the TESS observing baseline are listed in the
“Observing start/end” column, but the systems were not necessarily observed
continuously during this period. HD 135694 has 2-min cadence TESS light curves
from a total of 14 sectors: 14, 15, 16, 20, 21, 22, 26, 40, 41, 47, 49, 53,
56, and 60. TOI-1736 has 2-min cadence TESS light curves from a total of six
sectors: 18, 19, 25, 52, 58, and 59.

Of the 12 planets characterized in this work, 11 are detected in transit by
TESS. Table 1 summarizes the 2-min cadence TESS observations for each system
as they were available on the Mikulski Archive for Space Telescopes (MAST) on
2023-Mar-07 (i.e., up to and including TESS Sector 60). Each system was
observed in at least two sectors, with HD 135694 being observed in 14. The
photometry was processed by the TESS Science Processing Operations Center
pipeline (SPOC; Jenkins et al., 2016). All of the TESS data used in this paper
can be found in MAST: http://dx.doi.org/10.17909/y06k-3f04 (catalog
10.17909/y06k-3f04).

For all of our targets, there are no individual sources from _Gaia_ Data
Release 3 (DR3; Gaia Collaboration et al., 2016, 2022) within 20″ that cause
$>1\%$ dilution, nor does the combined flux of all DR3 sources within that
radius cause $>1\%$ dilution for any target. Furthermore, the SPOC data
products we use are already corrected for dilution from _Gaia_ sources per
_Gaia_ DR2 (Gaia Collaboration et al., 2018). In §3 we present high resolution
imaging observations that rule out significant dilution from unresolved
companions.666The results of TFOP imaging data, in addition to the results of
other reconnaissance observations, were used to inform TKS target selection,
as acknowledged by Chontos et al. (2022). This paper formally reports the
results of imaging observations of targets without such data already in the
literature. We discuss our light curve inspection and cleaning in §6 and our
transit modeling in §9.1.

##  3 High resolution imaging

To ensure that the planet transits were not subject to dilution from sources
not resolved by _Gaia_ , we used high resolution imaging (HRI) to place
contrast limits on potential nearby companions. HD 42813, TOI-669, and HIP
9618 all have high resolution images in the literature that rule out dilution
from nearby companions. We summarize the results of these observations in
§3.1. For the remaining five systems, HIP 8152, HD 25463, HD 135694, HD 6061,
and TOI-1736, we present new observations from Palomar-PHARO (Hayward et al.,
2001) and Keck-NIRC2 (Wizinowich et al., 2000). The Palomar-PHARO and Keck-
NIRC2 observations were obtained under the programs of PIs D. R. Ciardi and E.
J. Gonzales, respectively. A summary of the imaging observations from this
work can be found in Table 2 and sensitivity curves are shown in Figure 1.

Table 2: Imaging observations from this work

System | Instrument | Observation date | Filter | $t_{\mathrm{exp}}$ | $N_{\mathrm{exp}}$ | Resolution | Contrast at $0\arcsec.5$
---|---|---|---|---|---|---|---
|  | (UT) |  | (s) |  | (FWHM) | ($\Delta$ mag)
HIP 8152 | Palomar-PHARO | 2018-Dec-22 | Br-$\gamma$ | 9.9 | 15 | $0\arcsec.11$ | 6.0
HD 25463 | Keck-NIRC2 | 2020-Sep-09 | Br-$\gamma$ | 0.2 | 9 | $0\arcsec.05$ | 7.5
HD 135694 | Keck-NIRC2 | 2020-May-28 | Br-$\gamma$ | 1.0 | 9 | $0\arcsec.05$ | 7.0
HD 6061 | Keck-NIRC2 | 2020-May-28 |  $J$-cont | 1.2 | 9 | $0\arcsec.04$ | 7.0
HD 6061 | Keck-NIRC2 | 2020-May-28 | Br-$\gamma$ | 1.5 | 18 | $0\arcsec.05$ | 7.5
HD 6061 | Palomar-PHARO | 2020-Dec-05 |  $H$-cont | 1.4 | 15 | $0\arcsec.08$ | 7.1
HD 6061 | Palomar-PHARO | 2020-Dec-05 | Br-$\gamma$ | 1.4 | 15 | $0\arcsec.09$ | 6.7
TOI-1736 | Keck-NIRC2 | 2020-Sep-09 | Br-$\gamma$ | 0.5 | 9 | $0\arcsec.05$ | 6.5

Note. — $J$-cont: $\lambda_{0}=1.213$ $\mu$m and $\Delta\lambda=0.020$ $\mu$m.
$H$-cont: $\lambda_{0}=1.668$ $\mu$m and $\Delta\lambda=0.018$ $\mu$m.
Br-$\gamma$: $\lambda_{0}=2.169$ $\mu$m and $\Delta\lambda=0.0323$ $\mu$m.
Systems are listed in order of increasing TOI number.

Figure 1: AO imaging results for HIP 8152, HD 25463, HD 135694, HD 6061, and
TOI-1736 from our observations with Palomar-PHARO and Keck-NIRC2. Contrast
curves are shown in black with $1\sigma$ error envelopes in purple. The images
themselves are shown as the postage stamp insets. HIP 8152, HD 25463, HD
135694, and TOI-1736 all appear single. Keck-NIRC2 observations of HD 6061
were taken in both Br-$\gamma$ and $J$-cont, but only the former is shown
here. In the image of HD 6061, TIC 600433892 can be seen in the lower right
corner at a separation of $\approx 3\arcsec$. TIC 600433892 is fainter than HD
6061 by $6.2$ mag in the TESS bandpass, meaning that its dilution of HD 6061
b’s transits is about a 0.1% effect (i.e., much less than the uncertainty on
the stellar radius). HD 6061 and TIC 600433892 have consistent distances and
proper motions according to _Gaia_ DR3, meaning that the two stars are almost
certainly gravitationally bound. At a distance of 67 pc, their on-sky
separation of $3\arcsec$ translates to a sky-projected separation of about 200
AU. Figure 2: A 2MASS _JHK s_ color-color diagram. The dwarf branch, giant
branch, and brown dwarf loci are shown with green, blue, and red hashes,
respectively. Black dashed lines represent the direction of reddening induced
by extinction ($A_{V}$). HD 6061 and TIC 600433892 are overplotted as the
black circles with 1$\sigma$ error bars. TIC 600433892 is consistent with
being an M4/5V dwarf. We find that HD 6061 is consistent with being a G0
dwarf, which agrees with the classification from Cannon & Pickering (1993).

###  3.1 Literature Observations

HD 42813 and TOI-669 have Keck-NIRC2 observations from Schlieder et al.
(2021). HD 42813 was observed with Keck-NIRC2 on 2019-Mar-25 using the narrow-
band Br-$\gamma$ filter ($\lambda_{0}=2.169$ $\mu$m and $\Delta\lambda=0.032$
$\mu$m) and an integration time of 30 s. The star appears single and the Keck-
NIRC2 observation rules out companions of $\Delta 7.4$ mag at $5\sigma$
confidence at a separation of 0″.5. The TOI-669 Keck-NIRC2 observation was
taken on 2019-Jun-09 using the $K$ filter ($\lambda_{0}=2.196$ $\mu$m and
$\Delta\lambda=0.336$ $\mu$m) and an integration time of 10 s. The star
appears single, and the Keck-NIRC2 observation rules out companions of $\Delta
7.6$ mag at $5\sigma$ confidence at a separation of 0″.5.

Osborn et al. (2023) reported the discovery and confirmation of HIP 9618 b and
c using space-based photometry from TESS and CHEOPS (Benz et al., 2021) along
with radial velocity observations from CAFE, HARPS-N, and SOPHIE (Aceituno et
al., 2013; Cosentino et al., 2012; Perruchot et al., 2008). To rule out
dilution from nearby sources, the authors observed HIP 9618 with a variety of
optical speckle and near-infrared (NIR) adaptive optics (AO) instruments,
including Keck-NIRC2. The HRI shows no evidence of stellar companions within
1″.

###  3.2 Palomar-PHARO and Keck-NIRC2 observations

####  3.2.1 HIP 8152 (TOI-266)

Palomar Observatory HRI observations of HIP 8152 were made with the PHARO
instrument on the 5.1 m Hale telescope. Palomar-PHARO has a pixel scale of
$0\arcsec.025$ pix-1 for a total field of view of about $25\arcsec$.
Observations of HIP 8152 were taken on 2018-Dec-22 in the narrow-band
Br-$\gamma$ filter. Observations were acquired using the natural guide star AO
system P3K (Dekany et al., 2013) in the standard 5-point quincunx dither
pattern with steps of $5\arcsec$. Each dither position was observed three
times, with $0\arcsec.5$ positional offsets between each observation, for a
total of 15 frames. For HIP 8152, each frame had an integration time of 9.9 s,
amounting to a total on-source time of 149 s. No stellar companions were
detected.

####  3.2.2 HD 25463 (TOI-554)

Keck Observatory HRI observations of HD 25463 were made with the NIRC2
instrument on the 10 m Keck II telescope. Keck-NIRC2 was used in the narrow-
angle mode with a pixel scale of approximately $0\arcsec.01$ pix-1 and a full
field of view of about $10\arcsec$. Observations of HD 25463 were taken on
2020-Sep-09 in the narrow-band Br-$\gamma$ filter. Observations were acquired
using the natural guide star AO system in the standard three-point dither
pattern to avoid the lower left quadrant of the detector, which is typically
noisier than the other three quadrants. The dither pattern has a step size of
$3\arcsec$. Each dither position was observed three times, with $0\arcsec.5$
positional offsets between each observation, for a total of nine frames. For
HD 25463, each frame had an integration time of 0.2 s, amounting to a total
on-source time of 1.8 s. No stellar companions were detected.

####  3.2.3 HD 135694 (TOI-1247)

Keck-NIRC2 observations of HD 135694 were taken on 2020-May-28 following the
methods described in §3.2.2. Images were taken in the narrow-band Br-$\gamma$
filter. Each frame had an integration time of 1.0 s, amounting to a total on-
source time of 9 s. No stellar companions were detected.

####  3.2.4 HD 6061 (TOI-1473)

Keck-NIRC2 observations of HD 6061 were taken on 2020-May-28 following the
methods described in §3.2.2. Images were taken in both the $J$-continuum
($\lambda_{0}=1.213$ $\mu$m and $\Delta\lambda=0.020$ $\mu$m) and Br-$\gamma$
narrow-band filters. Observations of HD 6061 were taken in multiple filters
due to the visual observation of a nearby diluting source (TIC 600433892,
separation of $\approx 3\arcsec$) in order to further ascertain colors and the
likelihood of the nearby stellar object being bound. For the Br-$\gamma$
observations, six images were taken at each dither position, for a total of 18
frames. Each frame had an integration time of 1.2 s in $J$-cont and 1.5 s in
Br-$\gamma$, amounting to a total on-source time of 11 s in $J$-cont and 27 s
in Br-$\gamma$.

HD 6061 was also observed with Palomar-PHARO on 2020-Dec-05 following the
methods described in §3.2.1. HD 6061 was observed in both the Br-$\gamma$ and
$H$-continuum ($\lambda_{0}=1.668$ $\mu$m and $\Delta\lambda=0.018$ $\mu$m)
narrow-band filters. Each frame (in both Br-$\gamma$ and $H$-cont) had an
integration time of 1.4 s, amounting to a total on-source time of 21 s in each
filter. We discuss the nature of the stellar companion in §3.4. Other than TIC
600433892, no other stellar companions were detected.

####  3.2.5 TOI-1736

Keck-NIRC2 observations of TOI-1736 were taken on 2020-Sep-09 following the
methods described in §3.2.2. Images were taken in the narrow-band Br-$\gamma$
filter. Each frame had an integration time of 0.5 s, amounting to a total on-
source time of 4.5 s. No stellar companions were detected.

###  3.3 Palomar-PHARO and Keck-NIRC2 reduction

Both the Palomar-PHARO and the Keck-NIRC2 data were reduced using the same
methods. The science frames were flat-fielded and sky-subtracted. The flat
fields were generated from a median average of dark subtracted flats taken on-
sky. The flats were normalized such that the median value of the flats is
unity. The sky frames were generated from the median average of the dithered
science frames; each science image was then sky-subtracted and flat-fielded.
The reduced science frames were combined into a single co-added image using an
intrapixel interpolation that conserves flux, shifts the individual dithered
frames by the appropriate fractional pixels, and median-coadds the frames. The
final resolutions of the combined dithers were determined from the full-width
half-maximum (FWHM) of the point spread functions in the corresponding filter.

The sensitivities of the final combined AO image were determined by injecting
simulated sources azimuthally around the primary target every 20° at
separations of integer multiples of the central source’s FWHM (Furlan et al.,
2017). The brightness of each injected source was scaled until standard
aperture photometry detected it with $5\sigma$ significance. The resulting
brightness of the injected sources relative to each target set the contrast
limits at that injection location. The final $5\sigma$ limit at each
separation was determined from the average of all of the determined limits at
that separation. The uncertainty on the limit was set by the root-mean-square
(RMS) dispersion of the azimuthal slices at a given radial distance. The final
sensitivity curves are shown in Figure 1. For all targets, no stellar
companions were detected within $1\arcsec$.

###  3.4 TIC 600433892: a stellar companion to HD 6061

HD 6061 was observed in multiple filters with both Keck-NIRC2 and Palomar-
PHARO due to the presence of a nearby stellar companion. In Br-$\gamma$, the
companion, TIC 600433892, has a separation of $3\arcsec.06\pm 0\arcsec.20$ and
a position angle of $213\pm 1\degree$ E of N. TIC 600433892 is fainter than HD
6061 by $6.2$ mag in the TESS bandpass, meaning that its dilution of HD 6061
b’s transits, approximately a 0.1% effect, is negligible. According to _Gaia_
DR3, HD 6061 and TIC 600433892 have consistent distances to 1$\sigma$
($67.69\pm 0.07$ pc and $66.3\pm 1.4$ pc, respectively) and consistent proper
motions to 3$\sigma$ ($\mu_{\alpha}=-9.72\pm 0.01$ mas/yr and
$\mu_{\delta}=-10.09\pm 0.01$ mas/yr for HD 6061, $\mu_{\alpha}=-10.7\pm 0.5$
mas/yr and $\mu_{\delta}=-10.6\pm 0.2$ mas/yr for TIC 600433892). This implies
that the two stars are almost certainly gravitationally bound. At a distance
of 67 pc, $3\arcsec$ translates to a sky-projected separation of about 200 AU.

Following the methods of Ciardi et al. (2018), relative photometry was
conducted on the Keck-NIRC2 $J$-cont image and the Palomar-PHARO $H$-cont and
Br-$\gamma$ images to deblend the infrared magnitudes of the two stars (where
Br-$\gamma$ is taken to have a central wavelength that is sufficiently close
to $K_{s}$). The resulting Two Micron All Sky Survey (2MASS; Skrutskie et al.,
2006) $JHK_{s}$ color-color diagram suggests that TIC 600433892 is an M4/5V
dwarf (Figure 2). Following the reasoning in Ciardi et al. (2018), it is
unlikely that TIC 600433892 is a heavily reddened ($A_{V}>6$ mag using an
$R=3.1$ extinction law) early-F or late-A background star, given that the
entire line-of-sight extinction through the Galaxy is only $A_{V}\approx 2$
mag (Schlafly & Finkbeiner, 2011).

##  4 Determination of stellar properties

###  4.1 Stellar template observations

We used the High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) on
the 10 m Keck I telescope at the W. M. Keck Observatory on Maunakea to obtain
iodine-free spectra of each system at high resolution and S/N, which were used
to produce a deconvolved stellar spectral template (DSST) for each host. The
exposure parameters for each template are summarized in Table 3. Triple-shot
exposures of rapidly rotating B stars were taken with the iodine cell in the
light path immediately before and after the high-resolution templates were
collected in order to precisely constrain the instrumental point-spread
function (PSF). The data collection and reduction followed the methods of the
California Planet Search (CPS) as described in Howard et al. (2010).

###  4.2 Stellar characterization

We performed an initial stellar characterization of each host star using
SpecMatch-Emp (Yee et al., 2017) to constrain stellar effective temperature
($T_{\mathrm{eff}}$), metallicity ([Fe/H]), and stellar radius ($R_{*}$)
directly from the iodine-free Keck-HIRES template spectra. SpecMatch-Emp fits
stellar spectra between 5000 and 5800 Å in 100 Å segments using a linear
combination of spectral templates from a library of over 400 precisely
characterized FGKM stars.

To estimate the posteriors of the fundamental stellar parameters, we used
isoclassify (Huber et al. 2017; Berger et al. 2020) in grid mode with the
allsky dust map, which is an extinction model obtained via a combination of
the models from Drimmel et al. (2003), Marshall et al. (2006), and Green et
al. (2019). isoclassify infers marginal posteriors for stellar properties by
integrating over a grid of MIST isochrones (Choi et al., 2016). To inform the
isoclassify analysis, we input priors stemming from our SpecMatch-Emp results,
parallaxes from _Gaia_ DR3, and 2MASS $JHK_{s}$ magnitudes.777In the case of
HD 6061, the 2MASS $JHK_{s}$ magnitudes were deblended to account for the flux
from TIC 600433892 (see §3.4). Following Tayar et al. (2022), to account for
model-dependent systematic uncertainties, we inflated the errors on each host
star’s mass and radius by adding an additional 5% and 4% uncertainty,
respectively, in quadrature with the measurement error reported by
isoclassify. The final stellar parameters are summarized alongside planet
parameters in Appendix A.

Table 3: Keck-HIRES template observations

System | Date | $t_{\mathrm{exp}}$ | Decker | Airmass | S/N | $N_{\mathrm{exp}}$
---|---|---|---|---|---|---
| (UT) | (s) |  |  | (pix-1) |
HIP 8152 | 2019-Aug-18 | 476 | B3 | 1.34 | 210 | 1
HD 42813 | 2019-Oct-31 | 314 | B1 | 1.21 | 217 | 2
HD 25463 | 2019-Aug-18 | 26 | B3 | 1.52 | 200 | 3
TOI-669 | 2020-Jan-04 | 1692 | B3 | 1.28 | 214 | 1
HD 135694 | 2020-Mar-09 | 180 | B1 | 1.66 | 211 | 2
HIP 9618 | 2020-Jan-30 | 263 | B3 | 1.42 | 210 | 2
HD 6061 | 2019-Dec-28 | 180 | B3 | 1.15 | 213 | 2
TOI-1736 | 2020-Aug-11 | 187 | B3 | 1.59 | 211 | 1

Note. — B1 decker: 3″.5 $\times$ 0″.574, $R=$ 60,000. B3 decker: 14″$\times$
0″.574, $R=$ 60,000. S/N measured at 5500 Å. $N_{\mathrm{exp}}>1$ means that
consecutive exposures were taken and then combined to produce the final
template spectrum. For these cases, the $t_{\mathrm{exp}}$, airmass, and S/N
reported in this table are the median values across the $N_{\mathrm{exp}}$
observations. All template observations were acquired with a moon separation
of $>30\degree$.

##  5 Doppler follow-up

###  5.1 Keck-HIRES

We obtained high-resolution spectra of each target with Keck-HIRES to measure
precise radial velocities (RVs). RVs were determined following the procedures
of Howard et al. (2010). In brief, a warm cell of molecular iodine was placed
at the entrance slit during the RV observations (Butler et al., 1996). The
superposition of the iodine absorption lines on the stellar spectrum provides
both a fiducial wavelength solution and a precise, observation-specific
characterization of the instrument’s point spread function (PSF). As part of a
forward model, the spectrum is divided into about 700 pieces between
$\sim$5000–6000 Å, with each piece being 2 Å in width. For each piece, the
product of the DSST and the Fourier Transform Spectrograph (FTS) iodine
spectrum is convolved with the PSF to match the iodine-in observation. As one
of the free parameters, an RV for each piece of spectrum is produced. The
pieces are weighted using all observations of the star to produce a single RV
for each observation. Our Keck-HIRES Doppler observations are summarized in
Table 4 and the RV measurements can be found in Table 5.

###  5.2 APF-Levy

####  5.2.1 Data reduction and cleaning

For the brighter targets in our sample ($V<9.25$ mag) we also obtained high-
resolution spectra with the Levy spectrograph mounted on the 2.4 m Automated
Planet Finder telescope (APF; Vogt et al. 2014) at Lick Observatory on Mt.
Hamilton near San José, California. Though mounted on a much smaller
telescope, APF-Levy is complementary to Keck-HIRES in both latitude and
observing cadence. In the case of TOI-1736, APF-Levy observed periastron
passage for the giant planet (TOI-1736 c) while the system was inaccessible
from Maunakea. With its queue-based observing schedule and lower
oversubscription rate compared to Keck, APF-Levy can also typically observe
targets with higher cadence than Keck-HIRES.

The standard reduction pipeline used to compute RVs from APF-Levy spectra
follows the methods of Howard et al. (2010). As with our Keck-HIRES
observations, spectra were obtained with a warm cell of molecular iodine in
the light path. We used the Keck-HIRES DSSTs to compute RVs instead of
acquiring independent iodine-free template spectra with APF-Levy. Keck-HIRES
DSSTs have been shown to serve as effective replacements for APF-Levy
templates in the CPS Doppler reduction pipeline (e.g., Dai et al. 2020;
MacDougall et al. 2021; Dalba et al. 2022; Lubin et al. 2022) and provide an
efficient alternative to the long exposures that would otherwise be required
to achieve similar S/N on an iodine-free APF-Levy template.

To avoid using low quality APF-Levy RVs in our analysis, for each system we
inspected the distribution of APF-Levy RV errors as a function of S/N at 5500
Å. We placed a conservative maximum RV error threshold of three times the
median RV error for each target. Observations that landed above the error
threshold were removed and the APF-Levy RVs were recomputed using the cleaned
data set. For HD 135694, this resulted in removing four APF-Levy spectra, all
with $\sigma_{\mathrm{RV}}>6.3$ m/s. For HIP 9618, we removed six spectra, all
with $\sigma_{\mathrm{RV}}>5.8$ m/s. For HD 6061, we removed three spectra,
all with $\sigma_{\mathrm{RV}}>21.9$ m/s. For TOI-1736, we removed 13 spectra,
all with $\sigma_{\mathrm{RV}}>6.6$ m/s. The APF-Levy Doppler observations
used in our analysis are summarized in Table 4 and the RV measurements can be
found alongside the Keck-HIRES RVs in Table 5.

####  5.2.2 The case of HD 25463

For one system, HD 25463, the reduction methods we used to measure velocities
from the APF-Levy spectra were slightly different from the methods of Howard
et al. (2010) due to the star’s rapid rotation (for HD 25463, we measure a
sky-projected stellar rotational velocity of $v\sin i_{*}$ $=$ $11.6\pm 1.0$
km/s using SpecMatch-Syn; Petigura et al. 2017). For ease of reference, we
will refer to the methods in Howard et al. (2010) as the “default” reduction
pipeline. To measure the radial velocity of a star from a spectrum, the
default Doppler pipeline breaks the spectrum into small chunks and fits
stellar absorption lines chunk-by-chunk. The size of each chunk is determined
by a fixed pixel width. For Keck-HIRES, this pixel width translates to a chunk
width of about 2 Å in wavelength space. However, because APF-Levy has higher
spectral resolution than Keck-HIRES, this fixed pixel width translates to a
smaller chunk width in wavelength space. For reference, the W decker on APF-
Levy has $R=95,000$ (Vogt et al., 2014) while the B5 decker on Keck-HIRES has
$R=45,000$ (Vogt et al., 1994), where these deckers are typical for
observations of HD 25463. Using the default Doppler reduction pipeline on APF-
Levy spectra therefore results in less spectral information being contained in
each chunk than when it is applied to Keck-HIRES spectra.

This difference in the wavelength space width of each chunk is typically not
an issue for inactive, slowly-rotating stars (as is evident in the consistency
between the APF-Levy and Keck-HIRES RVs for a representative system such as
HIP 9618). However, for more rapid rotators ($v\sin i_{*}$ $\gtrsim 10$ km/s),
line broadening can conspire with the smaller chunk width to cause
catastrophic errors in the APF-Levy RV measurement process. This failure
happens because single stellar absorption lines become too broad to fit within
a single chunk. We observe this failure mode for HD 25463 when trying to
measure RVs from the APF-Levy spectra via the default method. To circumvent
this failure, for HD 25463’s APF-Levy spectra we compute RVs by fitting entire
echelle orders simultaneously instead of fitting small chunks in series. This
method also does not depend on an iodine-free template spectrum. Instead, we
simultaneously solve for the stellar template using all of the iodine-in
spectra. Save for these changes, the rest of the reduction is similar to the
default method.

We refer to this alternative method of computing the APF-Levy RVs as the
iGrand method. For completeness, Appendix B contains figures comparing the
default APF-Levy RVs to the iGrand RVs for HD 25463. It is clear that the
default APF-Levy RVs are inconsistent with the contemporaneous Keck-HIRES
measurements (the default APF-Levy RVs show nearly 100 m/s of scatter). In
contrast, the spread and uncertainties of the iGrand velocities are more in
line with expectations for a star of this magnitude ($V=6.9$ mag) and spectral
type ($T_{\mathrm{eff}}$ $=6200$ K).

As we did for the other targets that were observed with APF-Levy, we removed
low quality APF-Levy spectra of HD 25463 by setting a maximum RV error
threshold of three times the median APF-Levy iGrand RV error. This resulted in
removing 10 APF-Levy spectra, all with $\sigma_{\mathrm{RV}}>18.8$ m/s.

Table 4: Summary of RV observations

System | Instrument | First/last observation |  $N$ RVs (unbinned) | Median $t_{\mathrm{exp}}$ | Median S/N | Typical decker
---|---|---|---|---|---|---
|  | (UT) |  | (s) | (pix-1) |
HIP 8152 | Keck-HIRES | 2019-Aug-14/2022-Jul-25 | 94 (94) | 683 | 213 | C2
HD 42813 | Keck-HIRES | 2019-Sep-17/2022-Feb-22 | 71 (71) | 429 | 219 | B5
HD 25463 | Keck-HIRES | 2019-Aug-14/2022-Feb-22 | 97 (265) | 37 | 213 | B5
APF-Levy | 2019-Aug-02/2020-Nov-16 | 124 (152) | 592 | 108 | W
TOI-669 | Keck-HIRES | 2019-Nov-07/2021-Nov-24 | 61 (62) | 897 | 150 | C2
HD 135694 | Keck-HIRES | 2020-Jan-04/2022-Jul-09 | 79 (80) | 295 | 214 | B5
APF-Levy | 2019-Oct-25/2023-Jan-31 | 117 (150) | 1200 | 68 | W
HIP 9618 | Keck-HIRES | 2020-Jan-21/2022-Jul-25 | 61 (61) | 290 | 214 | B5
APF-Levy | 2020-Feb-19/2022-Nov-28 | 127 (206) | 1800 | 80 | W
HD 6061 | Keck-HIRES | 2019-Dec-16/2021-Aug-30 | 64 (64) | 205 | 212 | B5
APF-Levy | 2020-Jan-03/2022-Nov-22 | 56 (69) | 1200 | 69 | W
TOI-1736 | Keck-HIRES | 2020-Aug-02/2022-Jan-19 | 77 (79) | 226 | 203 | B5
APF-Levy | 2021-Feb-18/2023-Jan-28 | 257 (282) | 1800 | 85 | W

Note. — RVs are binned by 8 hrs. Keck-HIRES B5 decker: 3″.5 $\times$ 0″.861,
$R=$ 45,000. Keck-HIRES C2 decker: 14″$\times$ 0″.574, $R=$ 45,000. APF-Levy W
decker: 1″$\times$ 3″, $R=$ 95,000. S/N is measured at 5500 Å. All
observations were acquired with a moon separation of $>30\degree$.

Table 5: Radial Velocities and $S_{\mathrm{HK}}$ Values

System name | Time | RV | RV Unc. | $S_{\mathrm{HK}}$ |  $S_{\mathrm{HK}}$ Unc. | Inst.
---|---|---|---|---|---|---
| (BJD) | (m/s) | (m/s) |  |  |
HIP 8152 | 2458710.099141 | -4.38 | 1.41 | 0.164 | 0.002 | HIRES
$\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$

Note. — The RV and $S_{\text{HK}}$ measurements presented in this paper. Only
the first row of the table (which is sorted by system and then by observation
date) is shown here to inform its contents and format. BJD is reported using
the TDB standard (e.g., Eastman et al., 2010). Model-specific instrumental
offsets have not been applied to the RV values. The RV errors listed here
represent measurement uncertainty and have not been added in quadrature with
the corresponding instrument jitter values resulting from our models of the
data (see Appendix A). This table is available in its entirety online in
machine-readable format.

##  6 Light curve inspection and cleaning

Before applying our joint analysis of the photometry and RVs, we first
inspected and cleaned the TESS data. Using lightkurve (Lightkurve
Collaboration et al., 2018), we downloaded all of the TESS Presearch Data
Conditioning Simple Aperture Photometry (PDCSAP; Smith et al., 2012; Stumpe et
al., 2012, 2014) 2-min cadence data for each target, excluding data with NaN
values or data quality flags. We then normalized the data on a sector-by-
sector basis. We also applied the following analysis to the Simple Aperture
Photometry (SAP; Twicken et al. 2010; Morris et al. 2020) light curves for
each target. While the best-fitting transit parameters were nearly identical
between fits to the PDCSAP and SAP data, we generally found the SAP data
contained obvious spacecraft systematics and required more outliers to be
rejected.

###  6.1 Transit search

For each system, we searched for transits in the TESS PDCSAP light curve using
the box least squares method (BLS; Kovács et al. 2002). The signals reported
by the SPOC were identified in the transiting planet search pipeline
component, which employs an adaptive, noise-compensating matched filter
(Jenkins, 2002; Jenkins et al., 2010, 2020). We recovered all SPOC-reported
signals in the TOI catalog as of 2022-Oct-04, with a median S/N of 23 across
all of our BLS detections and with each detection having S/N $\gtrsim 10$.
After recovering the SPOC-reported signals, we masked the planet transits and
re-ran our BLS search but failed to find any other candidates.

In the case of TOI-554.02, our BLS search recovers the candidate’s shallow
transits with a slightly lower significance (S/N $\approx 8$), motivating, in
part, a more thorough investigation of the purported transit signal (see
§6.1.1). In addition, there are two instances where we identify transiting
planet candidates whose properties disagree with entries in the TOI catalog.
These are the so-called “duotransit” planets orbiting HIP 8152 and HIP 9618,
which we discuss in §6.1.2.

####  6.1.1 Statistical validation of TOI-554.02 (HD 25463 c)

In the hierarchy of exoplanet detection, statistical validation is typically
an intermediate step taken between planet candidacy and confirmation888For the
purposes of this discussion, we take planet “confirmation” to mean that the
planet’s mass has been measured to some fiducial precision. For the mass-
radius diagram in Figure 10, we show planets from the NASA Exoplanet Archive
whose mass measurements have better than 50% fractional precision. where
astrophysical false positive scenarios are systematically ruled out (e.g.,
Borucki et al., 2012; Morton et al., 2016). In the case of transiting planet
candidates, validation is used to statistically exclude the possibility that
the purported transit signal is in fact, for example, a background eclipsing
binary star system. Since most of the transiting planets in this work are at
least marginally ($\gtrsim 2.5\sigma$) detected with RVs, we bypass the
statistical validation step as the measurement of their host star’s Doppler
signal confirms their planetary nature. However, in two cases, we measure only
an upper limit on the planet mass. The first, HIP 9618 c, was externally
validated and confirmed by Osborn et al. (2023) using a combination of TESS
and CHEOPS photometry and CAFE, HARPS-N, and SOPHIE RVs (see §10.6.2 for
details). The second, TOI-554.02 (HD 25463 c) has not yet been confirmed, so
we take additional measures to statistically validate this planet.

A new Threshold Crossing Event (TCE) with $P=3.04$ d was detected by the
transit search of the SPOC Sectors 1–46 2-min light curve for HD 25463 (aka
TOI-554). An initial limb-darkened transit model was fitted (Li et al., 2019)
and a suite of diagnostic tests were conducted to help determine whether or
not the signal was planetary in nature (Twicken et al., 2018). The transit
signature passed all the diagnostic tests presented in the SPOC Data
Validation reports. The TESS Science Office (TSO) reviewed the vetting
information and issued an alert for TOI-554.02 on 2022-Apr-20.

The planet candidate is small ($R_{\mathrm{p}}$ $\approx 1.3$ $R_{\oplus}$
from the SPOC report), and the pipeline only detects its transit with S/N
$=8.5$. However, due to the candidate’s short orbital period and the system’s
four sectors of photometry, TESS has observed 29 purported transits. After
removing HD 25463 b’s transits from the light curve (as identified by our
initial BLS search of the system), we re-ran BLS, but TOI-554.02’s transit
signal was not immediately apparent. We narrowed the BLS period grid to look
for signals short of 10 d (down from 100 d) and increased the number of grid
points (by a factor of 2). We identified a peak in the BLS power spectrum with
S/N $\approx 8$ that corresponded to the SPOC-reported signal for TOI-554.02.
While each individual transit is not entirely obvious by eye, the phase-folded
transit shows a clear decrease in flux. We masked transits associated with
TOI-554.02 and re-ran the BLS search but found no additional transit-like
events.

As discussed above, since our RV observations only place an upper limit on the
mass of TOI-554.02 (see §10.3), we independently analyzed the TESS photometry
and our Keck-NIRC2 HRI with the planet validation framework TRICERATOPS
(Giacalone et al., 2021) to rule out astrophysical false positive scenarios
that might be responsible for TOI-544.02’s purported transit signal.
TRICERATOPS validates planets by simulating astrophysical false positives
arising from gravitationally bound stellar companions, chance-aligned
foreground or background stars, and known nearby stars that are blended with
the target in the TESS data. The marginal likelihoods of these false positive
scenarios are calculated and compared to that of the scenario where the signal
is caused by a planet transiting the target star. This calculation yields two
quantities: the false positive probability (FPP; the overall probability that
the signal is caused by something other than a planet transiting the target
star) and the nearby false positive probability (NFPP; the probability that
the signal is caused by a known nearby star that is blended with the target in
the TESS data). In order for a planet to be considered validated, it must
achieve ${\rm FPP}<0.015$ and ${\rm NFPP}<0.001$. To account for the intrinsic
stochasticity in its calculation, we ran TRICERATOPS 50 times on the same
dataset, obtaining ${\rm FPP}=(4.6\pm 0.2)\times 10^{-4}$ and ${\rm
NFPP}=(5.7\pm 0.9)\times 10^{-6}$. We find that the dominant contributor to
FPP is the STP scenario, which involves a gravitationally bound stellar
companion that hosts a transiting planet, although we note that this scenario
is unlikely due to the absence of evidence for a stellar companion in our
iodine-free spectra or RV data. Regardless, these values are sufficiently
small to consider the planet statistically validated. We also note that these
results are independent of the fact that we confirm HD 25463 b using RVs,
which makes it even more likely that TOI-554.02 is a true planet (Lissauer et
al., 2012; Guerrero et al., 2021). Hereafter we refer to TOI-554.02 as HD
25463 c.

####  6.1.2 Duotransit systems

HIP 8152 (TOI-266) and HIP 9618 (TOI-1471) both host two transiting planets,
with the sub-Neptunes HIP 8152 c and HIP 9618 c each having two transits in
widely time-separated sectors (these planets constitute a “duotransit”
scenario; Osborn et al., 2022). CHEOPS, in tandem with our Keck-HIRES RVs of
the systems, recently confirmed the correct period of HIP 8152 c ($P=19.61$ d;
via private communication with the CHEOPS team, point of contact H. Osborn)
and HIP 9618 c ($P=52.56$ d; Osborn et al., 2023). SPOC did not correctly
identify the transit signals of these planets in the TESS data: For HIP 8152,
TOI-266.02 is spuriously reported with $P=6.19$ d and $T_{\mathrm{c}}=1392.10$
BTJD. For HIP 9618, as last updated on 2022-Apr-20, TOI-1471.02 is listed as
having $P=683.33$ d (the time difference between the transit in Sector 17 and
the transit in Sector 42) and $T_{\mathrm{c}}=1779.19$ BTJD (correct). After
masking the transits of HIP 8152 b and HIP 9618 b as identified by BLS, we re-
ran our BLS search for both systems. In each case, the BLS power spectrum
contained peaks with comparable significance (S/N $\approx 10$) at the aliases
of the period allowed by the two widely time-separated transits. We masked the
planet c transits by hand in each system and ran another BLS search, but found
no additional transit-like events.

For HIP 9618, the SPOC pipeline originally excluded all data points in the
Sector 17 light curve beyond 1787.72 BTJD due to a high level of scattered
light from Earth, resulting in the exclusion of a second transit of HIP 9618 b
in Sector 17 near 1788 BTJD. This initially caused the pipeline to match the
first Sector 17 transit of HIP 9618 b with the Sector 17 transit of HIP 9618
c, and to report that TOI-1471.01 had $P=11.8$ d. It was not until later that
TFOP follow-up revealed these two transits were actually of different depth
and duration. To include the 1788 BTJD transit of HIP 9618 b in their
analysis, Osborn et al. (2023), hereafter O23, re-extract aperture photometry
for HIP 9618 starting from the 2-min cadence target pixel files. In place of
the PDC algorithm, they then use a custom light curve detrending method
similar to Vanderburg et al. (2019) in order to remove spacecraft systematics.
For the sake of homogeneity in our analysis of each system, we forgo
replicating their custom light curve extraction and detrending, meaning that
this work does not include the second transit of HIP 9618 b in Sector 17. We
note that our measured transit parameters for HIP 9618 b are all consistent
with the values reported by O23, and the primary reason for the difference in
the size of our uncertainties on the radius of HIP 9618 b ($\pm 0.04$
$R_{\oplus}$ from O23 and $\pm 0.13$ $R_{\oplus}$ from this work) is the
difference in the reported uncertainty on our stellar radius measurements
($\pm 0.005$ $R_{\odot}$ from O23 and $\pm 0.03$ $R_{\odot}$ from this work,
where our error estimate has been inflated according to Tayar et al. 2022). We
compare our results for HIP 9618 with those from O23 in detail in §10.6.2.

###  6.2 Light curve cleaning and initial transit fitting

After inspecting the TESS data for planet transits, we cleaned the photometry
with an outlier rejection scheme. First, for each sector we smoothed the
normalized TESS PDCSAP data in bins of 0.3 days with a cubic Savitzky-Golay
filter (Savitzky & Golay, 1964) and iteratively removed out-of-transit,
$>3\sigma$ outliers until convergence. We used the SPOC-reported orbital
period, time of transit, and transit duration to mask the planet transits,
save for HIP 8152 c and HIP 9618 c, since their orbital properties are
incorrect in the TOI catalog (see §6.1.2). For HIP 8152 c, we used the transit
duration of a photometry-only fit to the TESS data (the same as the model
described below) with a narrow Gaussian prior on the externally-confirmed
period of $P=19.61$ d. For HIP 9618 c, we used the transit duration from O23.
Figure 3 illustrates the results of the Savitzky-Golay filtering for HD
42813’s Sector 6 PDCSAP data. Across all systems, the number of outliers
removed per sector by the Savitzky-Golay filtering was $74\pm 14$. For each
system, this outlier rejection excluded $\lesssim 0.5\%$ of all of the
available TESS data. In each case, our iterative Savitzky-Golay filtering
routine converged in three iterations, save for HD 6061, which converged in
four iterations.

Next, we performed an additional outlier rejection step by fitting an initial,
photometry-only transit plus Gaussian process model (GP; e.g., Rasmussen &
Williams, 2006) to the data and iteratively removing 7$\sigma$ outliers about
the fit. The transit model was implemented with a quadratic limb darkening law
(Kipping, 2013) from starry (Luger et al., 2019) and the GP, used to remove
low-frequency stellar variability and instrumental systematics, was
constructed in celerite2 (Foreman-Mackey, 2018a). Following Kipping (2013),
the limb darkening coefficients are parameterized as
$q_{1}\equiv(u_{1}+u_{2})^{2}$ and $q_{2}\equiv 0.5u_{1}(u_{1}+u_{2})^{-1}$,
where $u_{1}$ and $u_{2}$ are the usual quadratic limb darkening coefficients.
The transit model is parameterized using $\ln P$, $T_{\mathrm{c}}$, $\ln
R_{\mathrm{p}}/R_{*}$, $b$, and $\ln T_{\mathrm{dur}}$. The parameters and
priors of this initial, photometry-only model are generally the same as for
the final, joint model of the photometry and RVs (Table 6). The main
difference between the two is that the joint model does not assume a circular
orbit and explicitly uses $\sqrt{e}\cos\omega$ and $\sqrt{e}\sin\omega$
instead of $\ln T_{\mathrm{dur}}$. For this initial, photometry-only model, we
placed a broad Gaussian prior on $\ln T_{\mathrm{dur}}$, the center of which
was the logarithm of the transit duration as reported in the TOI catalog when
accessed on 2022-Oct-04, and whose width was $\ln 10$ d. This initial transit
model also assumed no information about the stellar mass, since by employing a
circular orbit and fitting in terms of $T_{\mathrm{dur}}$ we imply a stellar
density. We elaborate on the differences between this initial, photometry-only
model and our joint model in §9.1.

The kernel of the GP used to flatten the light curve is in the form of an
overdamped stochastic harmonic oscillator (SHO). The power spectral density
(PSD) of the SHO kernel can be written as

$S(\omega_{f})=\sqrt{\frac{2}{\pi}}\frac{S_{0}\omega_{0}^{4}}{(\omega_{f}^{2}-\omega_{0}^{2})^{2}+\omega_{0}^{2}\omega_{f}^{2}/Q},$
(2)

where $\omega_{f}$ is the angular frequency, $\omega_{0}$ is the undamped
fundamental angular frequency, $S_{0}$ is the power at $\omega_{0}$, and $Q$
is the quality factor of oscillation. Following the reparameterization for the
SHO PSD from the celerite2 documentation, we define

$\rho=\frac{2\pi}{\omega_{0}},$ (3) $\tau=\frac{2Q}{\omega_{0}},$ (4)

and

$\eta=\sqrt{S_{0}\omega_{0}Q},$ (5)

where $\rho$ is interpreted as the undamped fundamental period of the
oscillator, $\tau$ is the characteristic timescale of the damping, and $\eta$
scales the amplitude of the GP (i.e., $\eta^{2}$ populates the diagonal of the
GP covariance matrix). Rewriting Equation 2 in terms of $\rho$, $\tau$, and
$\eta$, we have

$S(\omega_{f})=8\sqrt{2}\pi^{3/2}\frac{\eta^{2}}{\tau\rho^{2}}\Big{[}\big{(}\omega_{f}^{2}-(\frac{2\pi}{\rho})^{2}\big{)}^{2}+(\frac{2\omega_{f}}{\tau})^{2}\Big{]}^{-1}.$
(6)

The parameters and priors for this GP kernel are the same as used for the GP
that flattens the light curve in the joint model (see §9.1 and Table 6). We
note that a lower bound of 1 d was placed on $\rho$ and $\tau$ to prevent the
GP from overfitting the transits (see Figure 5).

For HIP 8152, TOI-669, HD 135694, and HIP 9618, no 7$\sigma$ outliers were
identified about this initial photometry-only model. HD 42813, HD 25463, and
HD 6061 each had one 7$\sigma$ outlier that was removed, and TOI-1736 had
three. For the systems for which we identified these outliers, we repeated the
initial transit model fitting with the outliers removed and found no remaining
outliers about the fit. The maximum a posteriori (MAP) values from our initial
models, save for $T_{\mathrm{dur}}$, were used as the starting values in the
MAP optimization routine for the corresponding parameters in the joint models.

Figure 3: An example of our Savitzky-Golay filtering procedure for HD 42813’s
Sector 6 PDCSAP photometry. The black points are the PDCSAP data, the green
line is the data after being smoothed by the Savitzky-Golay filter, and the
blue points are the in-transit data which are not subjected to the outlier
rejection. Outliers are marked in red.

###  6.3 Search for transit timing variations

For completeness, we searched the TESS data for any signs of transit timing
variations (TTVs; e.g., Hadden & Lithwick, 2017). We used the best-fitting
transit times and orbital periods from the initial photometry-only transit
model (above) as references for the expected transit times. We performed a MAP
fit of the photometry that was analogous to the initial transit model, but
now, for each planet, $\ln P$ and $T_{\mathrm{c}}$ were replaced with free
parameters for the midpoint of each individual transit. We placed a Gaussian
prior on each of the observed transit times centered at the expected time with
a width of 1 d.

For each of the 11 transiting planets in our sample, we found that the maximum
of the absolute difference between the observed and expected transit time
($O-C$) was $<20$ min and the median of these maximum values was about 1 min
across all planets. None of the $O-C$ time series show an obvious trend or
sinusoidal variation. HD 25463 b, TOI-669 b, and HD 135694 b each had a
maximum absolute value of $O-C$ between 10 and 20 min, and the scatter in
$O-C$ for each of these three planets was about 8 min. These systems may
warrant further investigation to determine whether the differences in the
observed and expected transit times are significant. However, in the absence
of a clear periodic TTV signal, we leave this work to future investigations.
For the three multi-transiting planet systems in our sample (HIP 8152, HD
25463, and HIP 9618), we note that none of the planet pairs have a near-
integer period ratio, so TTVs may not be expected for these systems a priori.
Given the lack of obvious evidence for TTVs in each system, we exclude them in
our joint model.

##  7 Search for RV trends and nontransiting companions

With all of the transits accounted for and the photometry cleaned, next we
conducted a systematic search for long-term RV trends and the full orbits of
nontransiting999We take “nontransiting” to mean that we did not observe a
transit in the TESS photometry. planetary signals in the RV time series. Long-
term RV trends are indicative of massive, distant companions, which are more
common for FGK hosts with close-in small planets (our sample) than for other
stars (Zhu et al., 2018; Bryan et al., 2019). Our analysis identifies two
systems with linear RV trends (HIP 9618 and TOI-1736) and one nontransiting,
super-Jovian-mass planet on a moderately eccentric orbit (TOI-1736 c).

###  7.1 RV trends

First, we attempted to determine which systems required a linear RV trend. We
used the Akaike Information Criterion (AIC; Akaike 1974) to choose between
models with and without a linear trend. The AIC is defined as

$\mathrm{AIC}=2k-2\ln\mathcal{\hat{L}},$ (7)

where $k$ is the number of free parameters in the model and
$\mathcal{\hat{L}}$ is the maximum of the likelihood function with respect to
the model parameters. In general, a lower AIC value is considered more
favorable. Let
$\Delta\mathrm{AIC}_{i}\equiv\mathrm{AIC}_{i}-\mathrm{AIC}_{\mathrm{min}}$,
where AICi is the AIC of the $i$th model under consideration and AICmin is the
lowest AIC value of all models considered. Burnham & Anderson (2004) provide
the following guidelines in interpreting $\Delta$AIC values:

  * •

If $\Delta\mathrm{AIC}_{i}<2$, the two models are nearly indistinguishable.

  * •

If $2<\Delta\mathrm{AIC}_{i}<10$, the $i$th model is disfavored.

  * •

If $\Delta\mathrm{AIC}_{i}>10$, the $i$th model is essentially ruled out.

When two models had $\Delta\mathrm{AIC}<4$, we chose the simpler model (e.g.,
even if including a linear RV trend reduces the AIC, if
$\Delta\mathrm{AIC}<4$, we adopted the model without a trend). There are only
two systems that demand a linear RV trend, HIP 9618 and TOI-1736. For these
systems we also attempted to include a quadratic term in addition to the
linear trend, and while the AIC could not rule out models with curvature,
there was no evidence to justify its inclusion.

###  7.2 Nontransiting companions

With linear RV trends either excluded or identified, we next used RVSearch
(Rosenthal et al., 2021) to search for the full orbits of nontransiting planet
candidates in the RV time series. RVSearch employs an iterative Generalized
Lomb-Scargle (GLS; Lomb, 1976; Scargle, 1982; Zechmeister & Kürster, 2009)
periodogram analysis to search for significant periodicity in the RV
residuals. Significance is determined following the detection methodology of
Howard & Fulton (2016), where an empirical false alarm probability (FAP)
threshold of 0.1% is computed via the Bayesian Information Criterion (BIC;
Schwarz 1978).

The BIC is defined as

$\mathrm{BIC}=k\ln n-2\ln\mathcal{\hat{L}},$ (8)

where $n$ is the size of the data and $k$ and $\mathcal{\hat{L}}$ are the same
as in Equation 7. While we could have also used the BIC to determine whether
or not to include linear trends in our RV models, simulation studies suggest
that for finite sample sizes, the BIC may be at risk of selecting very poor
models (Burnham & Anderson, 2004; Vrieze, 2012). In our analysis, the two
comparison statistics typically agreed and could be used relatively
interchangeably.

Before applying RVSearch to our RV time series, we removed the signals of the
transiting planets and, for HIP 9618 and TOI-1736, the linear RV trends.
RVSearch identified no signals above the 0.1% FAP threshold in the RV
residuals for all systems except HIP 8152, HD 135694, and TOI-1736. However,
for all but TOI-1736, it seems that the signals identified are related to the
RV window function. We discuss our interpretation of the detections below.

In the case of HIP 8152, RVSearch identified an eccentric signal ($e\approx
0.7$, $K\approx 4.5$ m/s) at $P=122$ d. This is likely the second harmonic of
the yearly observing alias (i.e., $122\times 3=366$), and we do not interpret
it as planetary in nature. Visually, it is clear that the Keck-HIRES
observations of HIP 8152 can be roughly grouped into three observing seasons
(see Figure 6, right). This is probably contributing to the power in the RV
window function around 365 d (Figure 7). It should also be mentioned that
$P=122$ d is an alias commonly seen in archival Keck-HIRES RV time series
(Rosenthal et al., 2021). Furthermore, at periastron passage, a planet with
$P=122$ d and $e=0.7$ would have a very close ($<0.01$ AU) encounter with the
orbit of HIP 8152 c ($P_{\mathrm{c}}$ = $19.6053\pm 0.0003$ d), suggesting
that such an architecture is not stable.

For HD 135694, RVSearch identified a moderately eccentric signal ($e\approx
0.3$, $K\approx 3.5$ m/s) at 45.6 d. However, like HIP 8152, there is
significant power related to the yearly alias in the periodogram of HD
135694’s RV window function (in this case, at $2\times 365$ d; Figure 20) and
the supposed period is likely a harmonic of this signal ($365.25/45.60=8.00$).
Therefore, we also interpret the 45.6 d signal as an artifact of our RV
sampling. If the $P=45.6$ d signal truly is a planet, however, it would not
cross orbits with HD 135694 b ($P_{\mathrm{b}}$ = $15.92346\pm 0.00002$ d).

The moderately eccentric orbit of the nontransiting super-Jovian, TOI-1736 c,
near $P=570$ d is visible in TOI-1736’s RV time series (Figure 25, right). For
completeness, we conducted a blind search for the orbit of TOI-1736 c after
removing the transiting planet, TOI-1736 b, and the system’s linear RV trend.
RVSearch recovers the orbit of TOI-1736 c with $P_{\mathrm{c}}$ $=573.6$ d,
$T_{\mathrm{c}}$ $=2272.4$ BTJD, $K_{\mathrm{c}}=195$ m/s,
$e_{\mathrm{c}}=0.37$, and $\omega_{\mathrm{c}}=162\degree$. Models of the RVs
that either replaced the linear RV trend with the partial orbit of an even
longer period giant planet or included a curvature term in addition to the
linear trend were not preferred by the AIC. After removing the linear RV trend
and the orbits of planets b and c, RVSearch failed to identify any other
signals above the 0.1% FAP threshold in the RV residuals.

##  8 Stellar activity considerations

Stellar activity mitigation is a key component of RV mass measurements for
small planets, especially when the stellar rotation period or one of its
harmonics is close to the period of the planet in question (e.g., Vanderburg
et al., 2016). Most of the hosts in our sample show little Ca II H and K
emission, implying that they are relatively inactive—this is in part why they
were chosen for Doppler monitoring (Chontos et al., 2022). Using our Keck-
HIRES spectra, we measure $\log_{10}R^{\prime}_{\text{HK}}$ (Middelkoop, 1982;
Noyes et al., 1984) for each system and find a median value across all eight
hosts of $-5.00$. For reference, over its magnetic cycle the Sun oscillates
between $\log_{10}R^{\prime}_{\text{HK}}$ $=-5.05$ and $-4.84$ at the solar
minimum and maximum, respectively (Meunier et al., 2010).

With each Keck-HIRES spectrum we also measured $S_{\text{HK}}$ values, which
trace Ca II H and K emission strength (Isaacson & Fischer 2010; Isaacson et
al. in prep). While photometry can act as a proxy for stellar activity (e.g.,
Aigrain et al., 2012; Haywood et al., 2014; Grunblatt et al., 2015), if
photometric and spectroscopic monitoring are not contemporaneous, the
connection between the time-varying activity signal during the two sampling
periods can be unclear (Kosiarek & Crossfield, 2020). Because the RV and
$S_{\text{HK}}$ measurements are simultaneous, they offer a real-time view of
the star’s behavior and serve as a useful supplement to the TESS photometry.

###  8.1 Correlated $S_{\text{HK}}$ values and RV residuals?

As a first step toward understanding the connection between stellar activity
and our RV measurements, we examined the correlation between the Keck-HIRES
$S_{\text{HK}}$ values and the Keck-HIRES RVs after the planetary-attributed
RV signals were removed. By “planetary-attributed,” we mean the RV signals
from transiting planets and, for HIP 9618 and TOI-1736, linear RV trends that
we assume are caused by distant giant companions. We attributed the strong
linear trends in the RV time series for HIP 9618 and TOI-1736 to distant
giants rather than stellar activity because for these quiet stars (which have
$\log_{10}R^{\prime}_{\text{HK}}$ = $-4.99\pm 0.05$ and
$\log_{10}R^{\prime}_{\text{HK}}$ = $-5.02\pm 0.05$, respectively) we would
expect the amplitude of the stellar activity signal to be on the order of a
few m/s (e.g., Wright, 2005; Wright et al., 2008). However, the change in RV
over the observing baselines for these systems is closer to 100 m/s and the
corresponding change in $S_{\text{HK}}$ value is $<0.02$, so it does not seem
like stellar activity could be responsible for the large RV trends. For
TOI-1736, we also removed the RV signal of TOI-1736 c—a nontransiting, massive
planet on a moderately eccentric orbit near $P=570$ d—because it is clearly
planetary.

With planetary signals removed from the RV time series, we calculated both the
Spearman rank-order correlation coefficient ($r_{\mathrm{Spearman}}$) and the
Pearson correlation coefficient ($r_{\mathrm{Pearson}}$) for the Keck-HIRES
$S_{\text{HK}}$ values and the Keck-HIRES RV residuals (e.g., Press et al.
1992). There were only two systems, HD 42813 and HD 6061, where the $p$-value
for either the Spearman or Pearson test was $<0.05$ (Figure 4). These
$p$-values may not be trustworthy given the relatively small sizes of the data
sets ($N_{\mathrm{HIRES}}=$ 71 and 64 for HD 42813 and HD 6061, respectively)
and concerns regarding $p$-value testing in general (e.g., Colquhoun 2014),
but the apparent correlation between the RV residuals and $S_{\text{HK}}$
values in these systems spurred further investigation. According to Ca II H
and K emission, HD 42813 is relatively inactive
($\log_{10}R^{\prime}_{\text{HK}}$ = $-4.98\pm 0.05$), while HD 6061 shows
signs of moderate activity ($\log_{10}R^{\prime}_{\text{HK}}$ = $-4.76\pm
0.05$). For the other systems in our sample, while there appears to be no
correlation between the RVs and the Keck-HIRES $S_{\text{HK}}$ values, we
still conducted a holistic examination of stellar activity.

Figure 4: HD 42813 (left) and HD 6061 (right) are the only two systems where
the $p$-value of either the Spearman or Pearson test (and in their cases,
both) was $<0.05$ for the Keck-HIRES RV residuals and $S_{\text{HK}}$ values.
In each panel the Keck-HIRES RVs are shown in blue and a linear least-squares
fit to the data is plotted in orange. We note that for HD 42813, we removed
the one observation with $S_{\text{HK}}$ $>0.19$ and refit the data, but still
find $p_{\mathrm{Spearman}}$ and $p_{\mathrm{Pearson}}$ are both $<0.05$. The
correlation between the RV residuals and the $S_{\text{HK}}$ values suggests
that stellar activity might be manifesting itself in the RVs, but we emphasize
that $p$-value testing can be unreliable (e.g., Colquhoun 2014). This test is
just one point of reference in our broader investigation into stellar activity
contamination in the RVs.

###  8.2 Periodogram analysis

Periodograms are a powerful tool for identifying periodic signals in time
series data, though, as we saw with our search for nontransiting planets in
§7, the physical interpretation of peaks in their power spectra should be
treated with care. Caution should also be exercised when searching for signals
across complementary data sets—Kosiarek & Crossfield (2020) find that their
periodogram and autocorrelation analyses of solar photometry correctly
identify the solar rotation period less than half of the time. With these
caveats in mind, we computed GLS periodograms to search for signs of
periodicity related to stellar activity in the out-of-transit (OoT) PDCSAP
TESS photometry and the Keck-HIRES $S_{\text{HK}}$ values for each system. As
we did for the periodograms of the RV residuals from RVSearch, we also
compared these to the periodogram of the RV window function to place purported
signals in context with our imperfect time sampling. Periodograms for each
system can be found in Appendix A. FAPs were calculated for these periodograms
following Baluev (2008).

In general, we do not see an obvious stellar rotation period in the TESS
PDCSAP (or SAP) photometry for any of our hosts. HD 6061 is the only system
whose PDCSAP light curve seems to exhibit rotational modulation by eye. The
periodogram analysis from Fetherolf et al. (2022) finds that HD 6061’s Sector
17 PDCSAP light curve is well fit by a single sinusoid with a period of
$4.8\pm 0.4$ d and an amplitude of about 0.2 ppt. This $P\approx 5$ d signal
coincides with the strong peak we see in our own periodogram of the Sector 17
and 57 OoT PDCSAP photometry (Figure 24). However, the activity-rotation
relation from Noyes et al. (1984) suggests that the rotation period of this
early-G dwarf should be closer to $P_{\mathrm{rot}}\approx 16$ d and the PDC
algorithm is known to suppress stellar activity signals with $P\gtrsim 10$ d.
A rotation period is not clear in HD 6061’s SAP light curve, which seems to be
heavily impacted by spacecraft and/or detector systematics (the light curve
has sharp ramps and a low-frequency trend). In the end, perhaps the $P\approx
5$ d signal is a harmonic of the rotation period, or an artifact of the
interplay between the PDC algorithm and the true astrophysical signal (if
any).

Similar to the case for the TESS photometry, none of the systems exhibit a
clear and obvious activity signal in the GLS periodogram of their Keck-HIRES
$S_{\text{HK}}$ values. For every system, the highest peak in the
$S_{\text{HK}}$ periodogram is either the nightly alias or related to the
yearly alias. We comment on other, seemingly inconsequential features of each
system’s $S_{\text{HK}}$ periodogram in §10.

###  8.3 The case of HD 6061

While most of the stars in our sample are relatively inactive
($\log_{10}R^{\prime}_{\text{HK}}$ $\lesssim-5.0$), HD 6061 is the only host
that would sit firmly among the “active” stars ($-5.0<$
$\log_{10}R^{\prime}_{\text{HK}}$ $<-4.3$) in the activity-rotation analysis
of Mamajek & Hillenbrand (2008). HD 6061 is a G0 dwarf (Cannon & Pickering,
1993) with moderate Ca II H and K emission ($\log_{10}R^{\prime}_{\text{HK}}$
= $-4.76\pm 0.05$). The HD 6061 Keck-HIRES RV residuals and $S_{\text{HK}}$
measurements are strongly correlated, indicating that a stellar activity
signal may be contaminating the RVs. What, then, is the star’s rotation period
and how does it compare to the orbital period of HD 6061 b ($P_{\mathrm{b}}$
$=5.25$ d)?

As mentioned in §8.2, there appears to be some sort rotational modulation in
HD 6061’s TESS photometry (with $P\approx 5$ d), but its connection to the
stellar rotation period is unclear. The activity-rotation relation from Noyes
et al. (1984) suggests $P_{\mathrm{rot}}\approx 16$ d, so perhaps the signal
in the TESS photometry is a harmonic of the true rotation period. Using
SpecMatch-Syn, we find $v\sin i_{*}$ $=$ $2.4\pm 1.0$ km/s. After combining
this with our stellar radius measurement ($R_{*}$ $=$ $1.03\pm 0.03$
$R_{\odot}$) and marginalizing over the inclination of the stellar spin axis,
HD 6061’s projected rotational velocity implies
$P_{\mathrm{rot}}=14^{+13}_{-9}$ d. While HD 6061’s true rotation period
remains uncertain, all of these clues suggest it is reasonable to expect that
$P_{\mathrm{rot}}$ or its harmonics are in the neighborhood of the orbital
period for planet b.

The GLS periodograms of the HD 6061 observations do not point to a clear and
obvious stellar rotation period, but they do appear to hint at unresolved
signals. After accounting for instrumental offsets and removing the Keplerian
signal of HD 6061 b, there are several peaks in the GLS periodogram of the RV
residuals. The highest peak is located at $P=5.5$ d and rises above the 1% FAP
level. There are also peaks at the 10% FAP level near $P=7$, 10, and 16 d. In
the GLS periodogram of the Keck-HIRES $S_{\text{HK}}$ values, the highest peak
between $P=2$ d and $P=100$ d (contributions from the window function dominate
the power spectrum beyond this range) is a peak at $P=13$ d that reaches the
10% FAP level. If the stellar rotation period is somewhere between $P=$ 12–17
d, this would seem to agree with the Noyes et al. (1984) estimate, our $v\sin
i_{*}$ measurement, and the peaks in the periodograms of the RV residuals and
$S_{\text{HK}}$ values.

To summarize, for each system we explored the possibility of stellar activity
contaminating the RV time series. We checked for a correlation between the
Keck-HIRES RV residuals and the $S_{\text{HK}}$ values. We also searched for
signals in the GLS periodograms of the TESS photometry, $S_{\text{HK}}$
values, and RV residuals. HD 6061 is the only system that seems to show an
activity signal. Though the principal period of the activity signal is not
entirely obvious, various estimates seem to suggest that the stellar rotation
period is in the neighborhood of $P=$ 12–17 d. In §9.3 we describe our formal
approach for including a Gaussian process model of stellar activity in our
joint model of the photometry, RVs, and $S_{\text{HK}}$ values.

##  9 Joint photometry, radial velocity, and activity modeling

Here we describe our method for deriving planet properties. In short, we used
a custom analysis pipeline based on the Python package exoplanet (Foreman-
Mackey, 2018b) to jointly model each system’s photometry, radial velocities,
and, if necessary, Keck-HIRES $S_{\text{HK}}$ stellar activity indicators. We
applied this framework homogeneously to each system in our sample. Our code
and worked examples are publicly available online (Akana Murphy, 2023).

We summarize our joint model of each system in Tables 6 and 7. All model
parameters had relatively broad priors, save for the stellar mass and radius,
whose informed Gaussian priors stemmed from our high-resolution spectroscopy
and isochrone modeling (see §4). The likelihood function of the joint model is
the product of the likelihood of the transit model and the RV model and, if
applicable, the $S_{\text{HK}}$ model, all of which assume Gaussian residuals.

###  9.1 Transits

We parameterize the transit portion of the joint model in terms of $\ln P$,
$T_{\mathrm{c}}$, $\ln R_{\mathrm{p}}/R_{*}$, $b$, and $\sqrt{e}\cos\omega$
and $\sqrt{e}\sin\omega$. As in our initial transit modeling, we use the
quadratic limb darkening law from Kipping (2013). When modeling photometry
alone, orbital eccentricity, argument of periastron, and impact parameter can
be highly degenerate for low to moderate S/N transits (Petigura, 2020). This
$e$-$\omega$-$b$ degeneracy can lead to multimodal MAP solutions and create
regions of very high curvature on the posterior surface (i.e., the dreaded
“funnel” geometry known to plague hierarchical models; Neal, 2003). One of the
main advantages that our joint model has over separate models of the
photometry and RVs is that in most cases, the RVs are able to quickly rule out
highly eccentric orbits for the transiting planets, thereby restricting the
$e$-$\omega$-$b$ phase space and alleviating this degeneracy.

We note that when _both_ the planet’s transit and RV signals are low to
moderate S/N, our joint model can still fall victim to the $e$-$\omega$-$b$
degeneracy because the RVs are not able to rule out cases of high $e$. For
example, when fitting our joint model to HD 42813, we found that a funnel
would form at moderate impact parameter ($b\gtrsim 0.7$) and moderate
eccentricity ($e\gtrsim 0.2$) because the planet’s RV detection is not
significant enough ($K/\sigma_{K}\approx 2.5$) to exclude orbits with large
$e$ and small $K$. In this case, we fixed $e\equiv 0$, which removed the
funnel and improved sampling reliability and performance. We encountered a
similar situation for HD 25463. In any case, RV-only models of these systems
show that the orbits are consistent with being circular (see §11).

Assumptions of circular orbits, parameterizing with transit duration (which
implies a stellar density), and importance sampling can be used to derive
constraints on $e$ and $\omega$ when combined with a known stellar density
from spectroscopy, for example. This strategy circumvents the $e$-$\omega$-$b$
degeneracy entirely (e.g., MacDougall et al., 2021). This is the approach we
used for our initial photometric model (see §6.2) when removing outliers so as
to avoid the $e$-$\omega$-$b$ degeneracy when fitting for the initial MAP
solution. While the parameterization in MacDougall et al. (2021) offers a
robust method of modeling the photometry alone in the presence of this
degeneracy, we chose to fit a joint model for simplicity rather than fitting
the photometry and RVs in series. As a sanity check, the posteriors of the
transit parameters resulting from of our joint model were all 1$\sigma$
consistent with the corresponding MAP values we found from the photometry-only
fit.

As we did with our initial transit modeling (§6.2), we fit the transit model
simultaneously with a GP using a kernel in the form of an overdamped SHO
(Equation 6) in order to flatten the light curve. To prevent the GP from
absorbing part of the transit signal, we enforced that the GP’s undamped
period ($\rho$) and damping timescale ($\tau$) must both be $>1$ d. For each
system, we also visually inspected each transit to ensure that the GP’s
prediction was sufficiently smooth across the transit duration. Figure 5
illustrates the simultaneous transit and GP fitting for HD 42813 b’s second
transit in Sector 6.

Figure 5: An example of our simultaneous transit and GP fitting for the second
transit of HD 42813 b in Sector 6. The PDCSAP data are shown in black and data
in 30-min bins are shown in red. The GP prediction across the transit (plus a
small global offset fit to the data) is shown as the green dashed line, and
the best-fitting transit model is shown as the blue line. We visually
inspected each transit across all systems to ensure that the GP prediction did
not absorb any of the transit signal.

###  9.2 Radial velocities

To describe the spectroscopic orbits of transiting planets, we used $\ln P$,
$T_{\mathrm{c}}$, $\sqrt{e}\cos\omega$ and $\sqrt{e}\sin\omega$, and $\ln K$,
where all but $\ln K$ were shared with the transit model. For each RV
instrument we also included an offset ($\gamma$) and RV jitter term
($\sigma$), where the latter is added in quadrature with the pointwise RV
measurement errors. As mentioned in §7, for each system we calculated the AIC
of models that included or excluded a linear RV trend, $\dot{\gamma}$. HIP
9618 and TOI-1736 are the only two systems where the AIC ruled out models
without a trend. For these systems, we also tried adding a quadratic term to
the background trend, but the AIC did not support including the curvature.

TOI-1736 is the only system for which our adopted joint model includes a
nontransiting planet. We treated the spectroscopic orbit of the nontransiting
planet, TOI-1736 c, in the same way as was done for transiting planets, save
for the fact that we broadened the Gaussian priors on $\ln P$ and
$T_{\mathrm{c}}$. The initial guesses for $\ln P$ and $T_{\mathrm{c}}$ were
taken from our RVSearch results for the system (see §7), and their priors were
given a width of $\ln$ 50 d and 100 d, respectively. The parameters and priors
for the RV signals of nontransiting companions in our joint model are found at
the top of Table 7.

Table 6: Joint model of the photometry and RVs

Parameter | Symbol | Units | Prior | Notes
---|---|---|---|---
_Light curve parameters_
Light curve mean offset | $\mu_{\mathrm{phot}}$ | ppt |  $\mathcal{N}$(0, 10) |
Log photometric jitter |  $\ln$ $\sigma_{\mathrm{phot}}$ |  $\ln$ ppt |  $\mathcal{N}$($\ln$ $s_{\mathrm{phot}}$, 2) | A
_RV instrument parameters_
Offset for RV instrument $i$ | $\gamma_{i}$ | m/s |  $\mathcal{U}$[-250, 250] |
Log jitter for RV instrument $i$ | $\ln\sigma_{\mathrm{RV},\>i}$ | m/s |  $\mathcal{N}$($\ln s_{\mathrm{RV,\>i}}$, 2)[$\ln 0.1$, $\ln 20$] | A
_Stellar parameters_
Limb-darkening parameter 1 | $q_{1}$ |  |  $\mathcal{U}$[0, 1] | B
Limb-darkening parameter 2 | $q_{2}$ |  |  $\mathcal{U}$[0, 1] | B
Stellar mass | $M_{*}$ | $M_{\odot}$ |  $\mathcal{N}$($M_{*}$, $\sigma_{M_{*}}$)[0, 3] | C
Stellar radius | $R_{*}$ | $R_{\odot}$ |  $\mathcal{N}$($R_{*}$, $\sigma_{R_{*}}$)[0, 3] | C
_Transiting planet parameters_
Log orbital period |  $\ln$ $P$ |  $\ln$ d |  $\mathcal{N}$($\ln$ $P_{\mathrm{TOI}}$, 1) | D
Time of inferior conjunction | $T_{\mathrm{c}}$ | d |  $\mathcal{N}$($T_{\mathrm{c},\>\mathrm{TOI}}$, 1) | D
Log occultation fraction |  $\ln$ $\frac{R_{\mathrm{p}}}{R_{*}}$ |  |  $\mathcal{N}$($\ln$ $\frac{R_{\mathrm{p}}}{R_{*}}_{\mathrm{TOI}}$, $\ln$ 10) | D
Impact parameter | $b$ |  |  $\mathcal{U}$[0, 1] |
$\sqrt{e}\cos({\omega})$ | $\xi_{1}$ |  |  $\mathcal{D}$($\xi_{1}$, $\xi_{2}$)[0, 1], VE($e|\bm{\theta}$) | E
$\sqrt{e}\sin({\omega})$ | $\xi_{2}$ |  |  $\mathcal{D}$($\xi_{1}$, $\xi_{2}$)[0, 1], VE($e|\bm{\theta}$) | E
Log RV semi-amplitude |  $\ln$ $K$ | m/s |  $\mathcal{N}$($\ln$ $s_{\mathrm{RV}}$, $\ln$ 50) | A
_Light curve GP hyperparameters_
Log GP amplitude | $\ln\eta_{\mathrm{phot}}$ |  $\ln$ ppt |  $\mathcal{N}$(0, 10) | F
Log GP undamped period | $\ln\rho_{\mathrm{phot}}$ |  $\ln$ d |  $\mathcal{N}$($\ln$ 10, $\ln$ 50)[$\ln$ 1, $\ln$ 200] | F
Log GP damping timescale | $\ln\tau_{\mathrm{phot}}$ |  $\ln$ d |  $\mathcal{N}$($\ln$ 10, $\ln$ 50)[$\ln$ 1, $\ln$ 200] | F

Note. — $\mathcal{N}$(X, Y) refers to a Gaussian distribution with mean X and
standard deviation Y. $\mathcal{N}$(X, Y)[A, B] refers to a bounded Gaussian
with mean X, standard deviation Y, and hard bounds at A and B.
$\mathcal{U}$[X, Y] refers to a uniform distribution inclusive on the interval
X and Y.

A: $\sigma_{\mathrm{phot}}$ is treated as a uniform pointwise flux measurement
error. $s_{\mathrm{phot}}$ refers to the sample standard deviation of the
PDCSAP light curve flux. $s_{\mathrm{RV,\>i}}$ refers to the same for the RVs
of instrument $i$.

B: The parameterization $q_{1}\equiv(u_{1}+u_{2})^{2}$ and $q_{2}\equiv
0.5u_{1}(u_{1}+u_{2})^{-1}$, where $u_{1}$ and $u_{2}$ are the usual quadratic
limb darkening coefficients, follows the prescription by Kipping (2013).

C: The bounded Gaussian priors on stellar mass and radius have centers and
widths corresponding to our derivation of the stellar parameters in §4.

D: For some parameter, $x$, $x_{\mathrm{TOI}}$ refers to the value of that
parameter as reported in the TOI catalog when accessed on 2022-Oct-04. The TOI
catalog contains erroneous orbital properties for HIP 8152 c and HIP 9618 c,
but the correct orbital ephemerides are known from CHEOPS observations.

E: $\mathcal{D}$($\xi_{1}$, $\xi_{2}$)[0, 1] refers to a uniform distribution
over the unit disk (i.e. $\sqrt{\xi_{1}^{2}+\xi_{2}^{2}}\leq 1$).
VE($e|\bm{\theta}$) refers to the mixture distribution from Van Eylen et al.
(2019) which is used as a prior on $e$ and whose hyperparameters,
$\bm{\theta}$, are fixed to the posterior medians from that work.

F: The hyperparameters of the GP used to flatten the light curve, which has a
kernel whose power spectral density is in the form of a stochastic harmonic
oscillator (SHO; see Equation 2). $\rho_{\mathrm{phot}}$ and
$\tau_{\mathrm{phot}}$, the undamped period and damping timescale of the SHO,
respectively, are forced to be $>1$ d to prevent the GP from overfitting low
signal-to-noise transits.

=2cm

Table 7: Additional parameters for nontransiting planets and joint GP modeling
of the RVs and Keck-HIRES $S_{\text{HK}}$ values

Parameter | Symbol | Units | Prior | Notes
---|---|---|---|---
_Nontransiting planet parameters_
Log orbital period |  $\ln$ $P$ |  $\ln$ d |  $\mathcal{N}$($\ln$ $P_{\mathrm{NT}}$, $\ln$ 50) | A
Time of inferior conjunction | $T_{\mathrm{c}}$ | d |  $\mathcal{N}$($T_{\mathrm{c},\>\mathrm{NT}}$, 100) | A
$\sqrt{e}\cos({\omega})$ | $\xi_{1}$ |  |  $\mathcal{D}$($\xi_{1}$, $\xi_{2}$)[0, 1], VE($e|\bm{\theta}$) | B
$\sqrt{e}\sin({\omega})$ | $\xi_{2}$ |  |  $\mathcal{D}$($\xi_{1}$, $\xi_{2}$)[0, 1], VE($e|\bm{\theta}$) | B
Log RV semi-amplitude |  $\ln$ $K$ | m/s |  $\mathcal{N}$($\ln$ $s_{\mathrm{RV}}$, $\ln$ 50) | C
_RV trends_
Linear RV trend | $\dot{\gamma}$ | m/s/d |  $\mathcal{N}$(0, 10) |
_$S_{\text{HK}}$ instrument parameters_
Offset for Keck-HIRES $S_{\text{HK}}$ values | $\gamma_{S_{\mathrm{HK}},\>\mathrm{HIRES}}$ |  |  $\mathcal{U}$[0, 1] |
Log jitter for Keck-HIRES $S_{\text{HK}}$ values | $\ln\sigma_{S_{\mathrm{HK}},\>\mathrm{HIRES}}$ |  |  $\mathcal{N}$($\ln s_{S_{\mathrm{HK,\>HIRES}}}$, 2) | C
_Instrument-specific GP hyperparameters_
GP amplitude of rotation term for Keck-HIRES $S_{\text{HK}}$ values | $\eta_{\mathrm{rot},\>S_{\mathrm{HK}},\>\mathrm{HIRES}}$ |  | Inv-$\Gamma$($\alpha=0.85$, $\beta=0.004$) | D
Log GP amplitude of exponential decay term for Keck-HIRES $S_{\text{HK}}$ values | $\ln\eta_{\mathrm{dec},\>S_{\mathrm{HK}},\>\mathrm{HIRES}}$ |  |  $\mathcal{N}$(0, 10) |
Log GP amplitude of rotation term for RV instrument $i$ | $\ln\eta_{\mathrm{rot},\>\mathrm{RV},\>i}$ |  $\ln$ m/s |  $\mathcal{N}$(0, 10) |
Log GP amplitude of exponential decay term for RV instrument $i$ | $\ln\eta_{\mathrm{dec},\>\mathrm{RV},\>i}$ |  $\ln$ m/s |  $\mathcal{N}$(0, 10) |
_RV and $S_{\text{HK}}$ shared hyperparameters for GP rotation term_
Log GP rotation period | $\ln P_{\mathrm{rot}}$ |  $\ln$ d |  $\mathcal{N}$($\ln P_{\mathrm{rot,\>guess}}$, $\ln 1.5$)[$\ln 1$, $\ln 60$] | E
Log quality factor of secondary mode | $\ln Q_{0}$ |  |  $\mathcal{N}$(0, 2) |
Log quality factor offset between primary and secondary modes | $\ln\delta Q$ |  |  $\mathcal{N}$(0, 2) |
Fractional amplitude of secondary mode relative to primary mode | $f$ |  |  $\mathcal{U}$[0.01, 1] |
_RV and $S_{\text{HK}}$ shared hyperparameters for GP exponential decay term_
Log undamped period of exponential decay term | $\ln\rho_{\mathrm{dec}}$ |  $\ln$ d |  $\mathcal{N}$($\ln$ 10, $\ln$ 50)[$\ln$ 1, $\ln$ 100] |
Quality factor of exponential decay term | $Q_{\mathrm{dec}}$ |  | $\equiv 1/\sqrt{2}$ | F

Note. — Notation in this table mirrors that in Table 6.

A: $P_{\mathrm{NT}}$ and $T_{\mathrm{c},\>\mathrm{NT}}$ refer to initial
guesses for the values of the period and time of inferior conjunction,
respectively, for a nontransiting planet’s orbit. These were taken as the MAP
values from our RVSearch results. Each prior’s width is such that we do not
expect the initial guesses to bias the best-fitting nontransiting planet
orbital parameters.

B: $\mathcal{D}$($\xi_{1}$, $\xi_{2}$)[0, 1] refers to a uniform distribution
over the unit disk (i.e. $\sqrt{\xi_{1}^{2}+\xi_{2}^{2}}\leq 1$).
VE($e|\bm{\theta}$) refers to the mixture distribution from Van Eylen et al.
(2019) which is used as a prior on $e$ and whose hyperparameters,
$\bm{\theta}$, are fixed to the posterior medians from Van Eylen et al.
(2019).

C: $s_{\mathrm{RV}}$ refers to the sample standard deviation of the RV
timeseries across all instruments. $s_{S_{\mathrm{HK,\>HIRES}}}$ is the sample
standard deviation of the Keck-HIRES $S_{\text{HK}}$ activity indices.

D: Inv-$\Gamma$ refers to the inverse Gamma distribution, the parameters of
which have been chosen to define the tails of the distribution such that
$p(x<0.001$ $)<0.01$ and $p(x>1$ $)<0.01$. This prior helps keep the amplitude
of this GP component positive though with a lighter tail near zero as opposed
to a Gamma distribution.

E: $P_{\mathrm{rot,\>guess}}$ is chosen using a periodogram analysis on a
case-by-case basis. The default prescription is to assign the value to the
period with the peak power in a GLS periodogram of the TESS photometry, though
this can be superseded in the context of signals in the RV residuals and/or
Keck-HIRES $S_{\text{HK}}$ values.

F: Fixing $Q_{\mathrm{dec}}\equiv 1/\sqrt{2}$, gives this (overdamped) SHO the
same power spectral density (PSD) as stellar granulation (Harvey, 1985;
Kallinger et al., 2014).

###  9.3 Gaussian process modeling of stellar activity

GPs are a popular tool for modeling correlated noise in RV data due to stellar
activity (e.g., Robertson et al. 2013; Haywood et al. 2014; Grunblatt et al.
2015; Kosiarek & Crossfield 2020). To further investigate contamination from
stellar activity in our RV time series beyond our exploratory analysis in §8,
we added a multidimensional GP to our joint model. We refer to this GP as
multidimensional because it is fit to the RVs and Keck-HIRES $S_{\text{HK}}$
values simultaneously. Each instrument (APF-Levy RVs, Keck-HIRES RVs, and
Keck-HIRES $S_{\text{HK}}$) is assigned its own GP kernel which shares all
hyperparameters with the other kernels, save for the GP amplitude (which we
denote with $\eta$). In addition to the GP, the Keck-HIRES $S_{\text{HK}}$
values are also fit with an offset and jitter term.

While APF-Levy and Keck-HIRES are both iodine-based RV instruments (i.e., they
measure RVs in the same spectral region) we use different amplitude
hyperparameters for their RV GP kernels. This is done as a catch-all to
account for systematic differences in the manifestation of the stellar
activity signal in their RV time series (e.g., perhaps those pertaining to
their difference in spectral resolution, differences in long-term spectrograph
stability, etc.). The practice of using independent GP amplitude
hyperparameters for separate RV instruments that cover a similar spectral
range is commonplace in the literature (e.g., Grunblatt et al., 2015; Kosiarek
et al., 2019, 2021).

The kernel of the multidimensional GP we used to model the stellar activity
signal in the RVs and Keck-HIRES $S_{\text{HK}}$ values is a mixture of three
terms, each of which has a PSD in the form of an SHO (the kernel is the sum of
celerite2’s SHOTerm and RotationTerm; Foreman-Mackey 2018a). The first term is
an overdamped oscillator, meant to describe exponentially-decaying behavior
such as spot evolution, and is the same as the kernel of the GP used to
flatten the light curve (Equation 6). The only difference is that we fix the
quality factor to $Q\equiv\frac{1}{\sqrt{2}}$ (which effectively sets the
characteristic damping timescale, $\tau$), since this gives the SHO the same
PSD as stellar granulation (Harvey, 1985; Kallinger et al., 2014). Plugging
into and rearranging Equation 6, for instrument $i$, we have

$S_{\mathrm{dec},\>i}(\omega_{f})=16\sqrt{2}\pi^{5/2}\frac{\eta^{2}_{\mathrm{dec},\>i}}{\rho^{3}}\Big{[}\big{(}\omega_{f}^{2}-(\frac{2\pi}{\rho})^{2}\big{)}^{2}\\\
+(\frac{2\sqrt{2}\pi\omega_{f}}{\rho})^{2}\Big{]}^{-1}.$ (9)

The second and third terms of the kernel are underdamped SHOs, with
fundamental frequencies corresponding to the stellar rotation period and its
first harmonic, respectively. For instrument $i$, the PSDs of these terms can
be written as

$S_{P_{\mathrm{rot}},\>i}(\omega_{f})=\sqrt{\frac{2}{\pi}}\frac{S_{1,\>i}\>\omega_{1}^{4}}{(\omega_{f}^{2}-\omega_{1}^{2})^{2}+\omega_{1}^{2}\omega_{f}^{2}/Q_{1}^{2}}$
(10)

and

$S_{P_{\mathrm{rot}}/2,\>i}(\omega_{f})=\sqrt{\frac{2}{\pi}}\frac{S_{2,\>i}\>\omega_{2}^{4}}{(\omega_{f}^{2}-\omega_{2}^{2})^{2}+\omega_{2}^{2}\omega_{f}^{2}/Q_{2}^{2}}.$
(11)

The hyperparameters of $S_{P_{\mathrm{rot}},\>i}$ and
$S_{P_{\mathrm{rot}}/2,\>i}$ are related via

$\displaystyle Q_{1}$ $\displaystyle=\frac{1}{2}+Q_{0}+\delta Q$ (12)
$\displaystyle\omega_{1}$ $\displaystyle=\frac{4\pi
Q_{1}}{P_{\mathrm{rot}}\sqrt{4Q_{1}^{2}-1}}$ (13) $\displaystyle S_{1,\>i}$
$\displaystyle=\frac{\eta^{2}_{\mathrm{rot},\>i}}{(1+f)\omega_{1}Q_{1}}$ (14)

and

$\displaystyle Q_{2}$ $\displaystyle=\frac{1}{2}+Q_{0}$ (15)
$\displaystyle\omega_{2}$ $\displaystyle=\frac{8\pi
Q_{1}}{P_{\mathrm{rot}}\sqrt{4Q_{1}^{2}-1}}$ (16) $\displaystyle S_{2,\>i}$
$\displaystyle=\frac{f\eta^{2}_{\mathrm{rot},\>i}}{(1+f)\omega_{1}Q_{1}},$
(17)

where $\eta_{\mathrm{rot},\>i}$ is the amplitude of
$S_{P_{\mathrm{rot}},\>i}+S_{P_{\mathrm{rot}}/2,\>i}$ relative to
$S_{\mathrm{dec},\>i}$, $Q_{0}$ is the quality factor minus $1/2$ for the
oscillator at $P_{\mathrm{rot}}/2$, $\delta Q$ is the difference between the
quality factors of the oscillators at $P_{\mathrm{rot}}$ and
$P_{\mathrm{rot}}/2$, $P_{\mathrm{rot}}$ is the primary period of variability
(meant to represent the stellar rotation period), and $f$ is the fractional
amplitude of the SHO at $P_{\mathrm{rot}}/2$ relative to the SHO at
$P_{\mathrm{rot}}$.

Putting it all together, the PSD of the GP kernel for instrument $i$ is the
sum of a term describing exponentially decaying behavior
($S_{\mathrm{dec},\>i}$) and a term describing periodic behavior
($S_{\mathrm{rot},\>i}$):

$\displaystyle S_{i}(\omega_{f})$
$\displaystyle=S_{\mathrm{dec},\>i}(\omega_{f})+\big{(}S_{P_{\mathrm{rot}},\>i}(\omega_{f})+S_{P_{\mathrm{rot}}/2,\>i}(\omega_{f})\big{)}$
(18)
$\displaystyle=S_{\mathrm{dec},\>i}(\omega_{f})+S_{\mathrm{rot},\>i}(\omega_{f}).$
(19)

The GP parameters and priors are summarized in Table 7. We experimented with
variants of this kernel (e.g., removing the exponentially-decaying term,
removing the first-harmonic term, removing the underdamped oscillators and
adding a second overdamped oscillator), but found that this kernel was best at
describing both the exponentially-decaying and periodic behavior.

###  9.4 Posterior estimation

We use No-U-Turn sampling (NUTS; Hoffman & Gelman, 2014), an adaptive form of
Hamiltonian Monte Carlo (HMC; Duane et al., 1987; Neal, 2012) implemented with
exoplanet and pymc3 (Salvatier et al., 2016a), to estimate the posterior
distributions of the parameters in our joint model. HMC sampling uses the
gradient of the posterior to help inform Markov transitions, enabling more
efficient exploration of high-dimensional posterior surfaces than brute-force,
guess-and-check methods like Metropolis-Hastings (Metropolis et al., 1953;
Hastings, 1970).

For each system, a NUTS sampler ran 8 parallel chains with each chain taking
at least 8000 “tuning” steps before drawing 6000 samples. Samples drawn during
the tuning period were discarded, similar to how various Markov Chain Monte
Carlo (MCMC) methods discard burn-in samples. The chains were concatenated to
produce a total of $N=4.8\times 10^{4}$ samples from the marginal posteriors
of each model parameter.

During the tuning stage the NUTS sampler optimizes hyperparameters, such as
step size, to meet a targeted sample acceptance rate (in our case, 90%) as it
explores the posterior surface. This can help prevent the sampler from taking
too large of a step while exploring a funnel on the posterior surface (Neal,
2003), where the gradient calculation can otherwise diverge and lead to biased
inference (Betancourt & Girolami, 2013). We found that the posterior
geometries of the models of some systems were more prone to regions of high
curvature than others, which prompted us to increase the number of tuning
steps to prevent divergences (hence why each chain for each sampler took “at
least” 8000 tuning steps).

In addition to being conscious of the number of tuning steps, for many joint
model parameters that are strictly non-negative (e.g., planet orbital period,
occultation fraction, RV semi-amplitude, etc.) we fit and explored the
posterior of the natural logarithm of the parameter of interest rather than
the parameter itself. We employed this parameterization because imposing hard
bounds on the domains of model parameters can encourage the formation of a
funnel.

We assessed convergence of the HMC sampling through multiple diagnostic
statistics. Vehtari et al. (2021) pointed out serious flaws with the standard
Gelman-Rubin statistic ($\hat{R}$; Gelman & Rubin 1992), which is
conventionally used to determine convergence for iterative stochastic
algorithms like MCMC. Following their prescription, we instead assessed
convergence by verifying a sufficiently small ($<1.001$) _rank-normalized_
$\hat{R}$ for each model parameter. In brief, a rank-normalized $\hat{R}$
statistic is computed by calculating $\hat{R}$ on the normalized, rank-
transformed chains of the parameter, rather than the values of the parameter
itself. To ensure the chains could offer reliable confidence intervals, we
also calculated the rank-normalized bulk and tail effective sample sizes from
Vehtari et al. (2021) for each of the marginal posteriors (roughly, the
effective sample sizes are the number of “independent” samples obtained in the
bulk and tails of the posterior). Vehtari et al. (2021) recommend that the
effective sample size should be larger than 400 in both the bulk and the tails
of the posterior. For every parameter we find the minimum between the bulk and
tail effective sample sizes was comfortably larger than the recommended
minimum threshold (typically we find $N_{\mathrm{eff}}\gtrsim 10^{4}$).

##  10 Joint modeling results

Here we describe the results of our stellar characterization, joint modeling,
and posterior estimation for each system. In this section we include figures
of the joint modeling results and periodograms for HIP 8152 to inform their
general format, but the rest can be found in Appendix A. Table 8 contains a
brief summary of the physical properties for all 12 planets. Tables of all
measured and derived planet properties, as well as stellar properties, can be
found for each system in Appendix A.

Figure 6: Our joint modeling results for HIP 8152. Left: The photometric
model. The top panel shows the PDCSAP light curve as the black points, with
the GP used to flatten the light curve plotted as the green line. Triangles
mark transits for planet b (blue) and c (orange). The middle panel shows the
flattened light curve in black and the best-fitting transit models for planet
b and c in blue and orange, respectively. Residuals are shown below. Phase-
folded light curves and residuals are shown in the bottom panels for planets b
(left) and c (right), with data in 30-min bins shown in red. The phase-folded
best-fitting transit models for each planet are shown as the slightly thicker
blue and orange lines, with 25 random posterior draws overplotted as the
thinner lines. Right: The RV model. The top panel shows the RV time series
with Keck-HIRES data in black and the RV model in blue. Residuals are shown
below. Phase-folded RV curves for planets b and c are shown in the bottom
panels. Red points are data binned in 0.125 units of orbital phase. The phase-
folded best-fitting RV models for each planet are shown as the slightly
thicker blue and orange lines, with 25 random posterior draws overplotted as
the thinner lines.

###  10.1 HIP 8152 (TOI-266)

HIP 8152 is an inactive G dwarf for which we report the discovery of twin sub-
Neptunes, HIP 8152 b and HIP 8152 c. Figure 6 summarizes our joint analysis
and Table 10 summarizes the system properties.

As mentioned in §6.1.2, before CHEOPS follow-up, HIP 8152 c constituted a
duotransit scenario. TESS observed only two transits of HIP 8152 c, one in
Sector 3 and one in Sector 33, initially obfuscating the planet’s true period.
In 2022-May, a GLS periodogram analysis of the RVs with planet b removed
showed a clear peak at $P=19.6$ d. We built a joint model of the TESS
photometry and Keck-HIRES RVs with MonoTools (Osborn et al., 2022) which
statistically ruled out all aliases other than the 19.6 d period. CHEOPS then
confirmed the period of HIP 8152 c by observing an additional transit in
2022-Aug (via private communication with the CHEOPS GTO team, point of contact
H. Osborn).

We do not include a GP to model stellar activity in the RVs. HIP 8152 is
inactive according to Ca II H and K emission and the Keck-HIRES RVs and
$S_{\text{HK}}$ values are not correlated. Furthermore, there are no peaks
rising above the 0.1% FAP threshold in a GLS periodogram of the Keck-HIRES
$S_{\text{HK}}$ values. As discussed in §7, the only peak that rises above the
0.1% FAP threshold in the RV residuals is near 120 d and is likely related to
the RV window function. Figure 7 shows the GLS periodograms for the system.

HIP 8152 b and c sit near the peak of the sub-Neptune distribution in the
mass-radius plane and are on orbits slightly short of a 2:1 mean-motion
resonance (MMR; $P_{\mathrm{c}}$/$P_{\mathrm{b}}$ $=1.8$). With nearly
identical physical properties, these planets are attractive candidates for
comparative studies in planet composition. Interestingly, HIP 8152 b, though
closer to the G dwarf host, is slightly less dense than HIP 8152 c.

Figure 7: GLS periodograms for HIP 8152 are shown in black. The vertical blue
and orange lines mark the orbital periods of planet b and c, respectively. The
dashed horizontal line indicates the 0.1% FAP threshold (Baluev, 2008). We
compute the GLS periodogram of the RV window function with a minimum period
slightly longer than 1 d to avoid the strong nightly alias which otherwise
skews the y-axis scale.

###  10.2 HD 42813 (TOI-469)

HD 42813 is an inactive, metal-rich K0V dwarf (Houk & Smith-Moore, 1988). We
report the discovery of a transiting sub-Neptune in the system, HD 42813 b.
Our adopted model assumes a circular orbit for HD 42813 b since posterior
estimation for a model that included eccentricity was hindered by the planet’s
modest RV detection (see §9.1). In any case, a simpler, RV-only model that
includes eccentricity finds that the planet’s orbit is consistent with being
circular (see §11). Figure 13 summarizes our joint analysis and Table 11
summarizes the system properties.

While the Keck-HIRES RVs appear to be slightly correlated with the Keck-HIRES
$S_{\text{HK}}$ values (see Figure 4, left), we decided to leave out a GP fit
to the RVs and $S_{\text{HK}}$ values in our joint model. HD 42813 is
seemingly inactive according to Ca II H and K emission to begin with.
Furthermore, we do not find significant peaks in the GLS periodogram of the
residuals of the Keplerian-only RV model. There are no peaks rising above the
0.1% FAP threshold in the GLS periodogram of the Keck-HIRES $S_{\text{HK}}$
values save for broad peaks that are consistent with the yearly alias and its
harmonics. In addition, the GLS periodogram of the RV window function shows
significant power at the yearly and lunar monthly aliases, making it difficult
to disentangle a would-be activity signal in the RV residuals from artifacts
of the RV time sampling (Figure 14). We fit a joint model to the data that
included a GP and recovered a mass measurement for HD 42813 b that was
entirely consistent with the results of our Keplerian-only model. We therefore
adopt a joint model of the TESS photometry and Keck-HIRES RVs that does not
include a GP fit to the RVs and $S_{\text{HK}}$ values.

With its seemingly low density, HD 42813 b is a very attractive candidate for
atmospheric follow-up (TSM $>89$ at 98% confidence). However, continued RV
monitoring is first required to refine the planet’s mass measurement.

###  10.3 HD 25463 (TOI-554)

As an F8 dwarf (Cannon & Pickering, 1993), HD 25463 has the earliest spectral
type of all the hosts in our sample. We report the discovery of two transiting
planets orbiting HD 25463: the sub-Neptune HD 25463 b and the super-Earth HD
25463 c. Similar to HD 42813, we enforced circular orbits for planets b and c
to avoid divergences during the HMC posterior estimation, though an RV-only
model that included eccentricity found that the orbits were consistent with
being circular (see §11). Figure 15 summarizes our joint analysis and Table 12
summarizes the system properties.

HD 25463 is seemingly inactive according to Ca II H and K emission
($\log_{10}R^{\prime}_{\text{HK}}$ $=$ $-5.26\pm 0.05$). However, we were
prompted to explore a model that included a GP for activity mitigation because
of the star’s relatively rapid rotation. Using SpecMatch-Syn we find $v\sin
i_{*}$ $=$ $11.6\pm 1.0$ km/s. When combining this with our stellar radius
measurement and marginalizing over the inclination of the stellar spin axis,
the star’s projected rotation velocity implies
$P_{\mathrm{rot}}=4.3^{+1.7}_{-2.8}$ d. If this simplistic estimate is to be
trusted, the stellar rotation period is very close to the orbital periods of
the transiting planets ($P_{\mathrm{b}}$ $=7.0$ d and $P_{\mathrm{c}}$ $=3.0$
d). Vanderburg et al. (2016) highlight how stellar rotation can confuse the
search for the Doppler signals of planets when $P_{\mathrm{rot}}$ and its
first harmonic are in the neighborhood of the planets’ orbital periods.
Despite this concern, we do not find peaks above the 0.1% FAP threshold in a
GLS periodogram of the RVs after removing planets b and c. There are also no
peaks that rise above the 0.1% FAP threshold in a GLS periodogram of the Keck-
HIRES $S_{\text{HK}}$ values (see Figure 16). Nevertheless, we added a GP to
our joint model of the system with a Gaussian prior of $\mathcal{N}(4,1.5)$ d
on $P_{\mathrm{rot}}$, where 1.5 d is the Gaussian’s standard deviation. The
GP-enabled model finds best-fitting masses for planets b and c that are
consistent with the posterior estimates of the non-GP model, but the HMC
sampling had difficulty converging due to the large number of additional model
parameters introduced by the GP kernel. We therefore adopt the non-GP model,
whose HMC sampling does converge.

Given the brightness of the system ($V=6.9$ mag, $J=6.0$ mag), HD 25463 is
highly amenable to both ground- and space-based follow-up. With planets on
opposite sides of the radius valley (Fulton et al., 2017; Van Eylen et al.,
2018), HD 25463 represents an opportunity for comparative studies in
atmospheric mass loss. Though HD 25463 is too bright for single object
slitless spectroscopy (SOSS) with _JWST_ ’s Near Infrared Imager and Slitless
Spectrograph (NIRISS; brightness limit of $J=6.5$ mag), the system represents
an attractive target for HST. The mass measurement precision for the planets
must be improved, however, before they are subjected to detailed atmospheric
characterization (Batalha et al., 2019).

###  10.4 TOI-669

TOI-669 is an inactive G dwarf for which we report the discovery of a hot
transiting sub-Neptune, TOI-669 b. Figure 17 summarizes our joint analysis and
Table 13 summarizes the system properties.

We do not include a GP to model a stellar activity signal in the RVs and
$S_{\text{HK}}$ values. TOI-669 is inactive according to Ca II H and K
emission, and the Keck-HIRES RVs and $S_{\text{HK}}$ values are not
correlated. Furthermore, there are no peaks rising above the 0.1% FAP
threshold in a GLS periodogram of the Keck-HIRES $S_{\text{HK}}$ values or the
Keck-HIRES RV residuals (see Figure 18).

A less significant peak near 9.6 d is visible in the periodogram of the RV
residuals, but it is unclear whether the signal is planetary or related to the
RV window function, which has significant power near 180 d and 25 d. If the
$P=9.6$ d signal is in fact a planet, assuming a circular orbit, $b>1$ would
imply an orbital inclination of $i_{\mathrm{p}}$ $<86.8\degree$. For
reference, the orbit of TOI-669 b has $i_{\mathrm{b}}=$ $88.7\pm 1.0\degree$.
A two-planet fit to the RVs using RadVel (Fulton et al., 2018) does not result
in a significant detection of the $P=9.6$ d signal (it finds
$M_{\mathrm{p}}\sin i_{\mathrm{p}}=5.0\pm 2.6$ $M_{\oplus}$ for a Keplerian at
$P=9.61\pm 0.52$ d; the resulting mass of TOI-669 b in this two-planet fit is
consistent with our adopted joint model). The $\Delta$AIC between the one- and
two-planet RadVel models is $<1$, so there does not appear to be evidence for
including the $P=9.6$ d signal.

TOI-669 is a relatively bright G dwarf ($J=9.6$ mag) whose hot sub-Neptune
($T_{\mathrm{eq}}$ $=$ $1235\pm 37$ K, $S_{\mathrm{p}}$ $=$ $388.2\pm 47.1$
$S_{\oplus}$) lies just outside of the “sub-Neptune desert” (planets with
$2.2<R_{\mathrm{p}}<3.8$ $R_{\oplus}$ and $S_{\mathrm{p}}$ $>650$
$S_{\oplus}$; Lundkvist et al. 2016). TOI-669 b’s mass and radius measurements
place it at the mode of the sub-Neptune mass-radius distribution.

###  10.5 HD 135694 (TOI-1247)

HD 135694 is a K0 dwarf (Cannon & Pickering, 1993). We report the discovery of
a warm sub-Neptune in the system, HD 135694 b. Figure 19 summarizes our joint
analysis and Table 14 summarizes the system properties.

We do not include a GP to model stellar activity in the RV time series. HD
135694’s Ca II H and K emission ($\log_{10}R^{\prime}_{\text{HK}}$ $=$
$-4.99\pm 0.05$) indicates that the star is relatively inactive and the Keck-
HIRES RV residuals and $S_{\text{HK}}$ values do not appear to be correlated.
Although there are 14 sectors of TESS photometry available, a stellar rotation
period is not readily apparent in either the PDCSAP or SAP light curve. There
is a strong peak ($<0.1\%$ FAP) in the RV residuals at 45.6 d (see Figure 20),
but we attribute this power to the RV window function given that 45.6 is a
near-perfect divisor of 365.25. It is unclear why our mass measurement is so
imprecise (about 2.7$\sigma$) when we have nearly 200 RV measurements between
APF-Levy and Keck-HIRES. Below, we discuss scenarios that may be the cause of
model misspecification.

After the $P=45.6$ d peak, second-highest peak in the GLS periodogram of the
RV residuals rises above the 10% FAP threshold and is located at about
$P=32.5$ d. This period is not a clear harmonic of the yearly alias or the
lunar monthly alias, but it is just about twice the period of HD 135694 b
($32.5/15.9=2.0$). The highest peak short of $P=100$ d in the GLS periodogram
of the Keck-HIRES $S_{\text{HK}}$ values (beyond $P=100$ d the GLS power is
dominated by contributions from the RV window function) is at about $P=31$ d
and also rises above the 10% FAP threshold. While the $P=32.5$ d peak in the
RV residuals is not overwhelmingly significant, it could represent either a
nontransiting planet in a near 2:1 MMR with HD 135694 b, or the stellar
rotation period. The latter explanation seems slightly more preferable given
the $P\approx 31$ d peak in the GLS periodogram of the Keck-HIRES
$S_{\text{HK}}$ values. Furthermore, the activity-rotation relation from Noyes
et al. (1984) suggests that HD 135694 has $P_{\mathrm{rot}}\approx 29$ d,
making a rotation period of roughly 31–33 d seem reasonable for this K0 dwarf.

For completeness, we added a GP component to our joint model following the
methodology in §9.3. We placed a Gaussian prior on the GP rotation period at
32.5 d with a width of 1.5 d. The only other difference between this model and
our adopted model of the photometry and RVs is that we forced HD 135694 b’s
orbit to be circular so as to prevent the HMC sampling from diverging (which
it tended to do when allowing $e_{\mathrm{b}}$ and $\omega_{\mathrm{b}}$ to
float). This GP-enabled model returned $M_{\mathrm{b}}$ $=6.1\pm 2.1$
$M_{\oplus}$, in agreement with the results of our adopted model.

We also explored the idea that the signal near 32 d could be a nontransiting
planet in a near 2:1 MMR with planet b. We fit a joint model where the GP in
the model above was replaced with a nontransiting planet on a circular orbit.
This model finds $M_{\mathrm{b}}$ $=6.6\pm 2.1$ $M_{\oplus}$ for HD 135694 b
(which agrees with the results of our adopted, one-planet model) and
$M_{\mathrm{p}}\sin i_{\mathrm{p}}=10.6\pm 2.6$ $M_{\oplus}$ for the signal at
$P=32.5$ d. The AIC favors the one-planet plus GP model over the (adopted)
one-planet model, which itself is favored over the two-planet model (all at
the $\Delta$AIC $>10$ level). It should be noted, however, that GPs can be
susceptible to overfitting (e.g., Blunt et al., 2023), which can muddle the
interpretability of Bayesian model comparison statistics. In this context, the
AIC’s preference for the GP-enabled model is not entirely surprising. Setting
aside the AIC comparison, if there was a nontransiting planet at 32.5 d,
assuming a circular orbit, $b>1$ for the planet would imply $i_{\mathrm{p}}$
$<88.6\degree$. For reference, we find that planet b’s orbit has
$i_{\mathrm{b}}=$ $89.1\pm 0.5\degree$. Simulations similar to those conducted
by Lubin et al. (2022) could be used to place a lower limit on the inclination
of the potential nontransiting planet, but these are beyond the scope of this
work.

In summary, we cannot rule out the possibility that the $P=32.5$ d signal
represents either the stellar rotation period or a nontransiting planet.
However, the signal’s ambiguity, combined with the star’s lack of Ca II H and
K emission, encouraged us to adopt a model of the photometry and RVs that does
not use a GP for stellar activity mitigation and does not include
nontransiting planets. In any case, our experimentation with various models of
the data reassures us that the mass measurement of planet b is seemingly
insensitive to our choice of model. HD 135694 is a bright ($V=9.1$ mag,
$J=7.9$ mag) K0 dwarf that adds another planet to the mode of the sub-Neptune
mass-radius distribution. Continued Doppler monitoring is required to refine
the planet’s mass measurement precision.

###  10.6 HIP 9618 (TOI-1471)

####  10.6.1 Joint analysis from this work

HIP 9618 is a G5 dwarf (Cannon & Pickering, 1993). The system is host to two
warm sub-Neptunes, HIP 9618 b and HIP 9618 c. We robustly detect a linear RV
trend, indicating that there is also a distant, massive companion in the
system. Figure 21 summarizes our joint analysis and Table 15 summarizes the
system properties.

We exclude a GP fit to the RVs and Keck-HIRES $S_{\text{HK}}$ values from our
adopted model. HIP 9618 is nominally inactive based on its Ca II H and K
emission levels ($\log_{10}R^{\prime}_{\text{HK}}$ $=$ $-4.99\pm 0.05$) and
the Keck-HIRES RV residuals and $S_{\text{HK}}$ values are not correlated.
Stellar activity does not seem to be a concern for HIP 9618 according to our
GLS periodograms (Figure 22). There are no peaks that rise above the 0.1% FAP
threshold in the RV residuals of our joint model. In the GLS periodogram of
the Keck-HIRES $S_{\text{HK}}$ values, there is a peak just long of 30 d that
rises above the 0.1% FAP level, but it is unclear whether or not the power is
related to the window function—the periodogram also shows significant power
near 180 d and at $P>700$ d. Even if the signal near $P=30$ d in the
$S_{\text{HK}}$ values is related to the stellar rotation period for this late
G dwarf, the lack of power in the periodogram of the RV residuals indicates
that activity is not greatly impacting the planet mass measurements.

The nature of the massive companion causing the linear RV trend is uncertain.
We find that the companion must have $M\sin i_{\mathrm{p}}\gtrsim 4.7$
$M_{\mathrm{Jup}}$ and $a\gtrsim 5.0$ AU by making the following simplifying
assumptions: (1) the RV trend is caused by a single companion, (2) the
companion’s orbit is circular, (3) the companion’s orbital period is greater
than four times our APF-Levy and Keck-HIRES RV baseline (about 1042 d), and
(4) the RV semi-amplitude of the companion’s orbit is greater than the
$\Delta$RV caused by the trend over the baseline (about 62 m/s). A more
detailed investigation is beyond the scope of this work—in §10.6.2 we
summarize the results from O23, who conduct a thorough analysis of the trend
that includes constraints from RVs, astrometry, and direct imaging. Continued
RV monitoring is required to reveal the true nature of the distant companion.

HIP 9618 is perhaps the most exciting system in our sample for atmospheric
follow-up with _JWST_. As noted by O23, HIP 9618 is one of only five multi-
transiting systems with $K_{s}<8$ mag to host a planet with $P>50$ d. Our mass
constraints translate to TSM values of $161\pm 43$ and $>87$ for planet b and
c, respectively (with the lower limit for planet c reflecting 98% confidence).
These values place both planets above the Kempton et al. (2018) TSM cutoff
($>84$) for follow-up of planets with $2.75<$ $R_{\mathrm{p}}$ $<4$
$R_{\oplus}$. Moreover, HIP 9618 b’s TSM estimate places it in the top
quartile of all planets in its radius range (top quartile cutoff of TSM
$>146$).

Figure 8: The RV portion of a joint model for HIP 9618 that includes the CAFE,
HARPS-N, and SOPHIE RVs from O23 in addition to the APF-Levy and Keck-HIRES
data from this work. The figure description is the same as for Figure 6. This
RV model also includes a curvature term in addition to the linear trend, in
order to mimic the adopted solution from O23. The results of this model are
consistent with our joint model of the APF-Levy and Keck-HIRES RVs alone.

####  10.6.2 Comparison with the O23 results

O23 first reported the discovery and confirmation of HIP 9618 b and c using
space-based photometry from TESS and CHEOPS in combination with a total of 49
RVs from CAFE, HARPS-N, and SOPHIE. The authors report masses of
$M_{\mathrm{b}}$ $=10.0\pm 3.1$ $M_{\oplus}$ and $M_{\mathrm{c}}$ $<18$
$M_{\oplus}$ (at 3$\sigma$ confidence).

The only significant difference between the RV model presented in this work
and the adopted model from O23 is that, in addition to including a linear
trend in their model of the RVs, O23 also include a quadratic term
($\ddot{\gamma}$). O23 find $\dot{\gamma}=-0.067\pm 0.001$ m/s/d and
$\ddot{\gamma}=-4.99\times 10^{-6}\pm 6.4\times 10^{-7}$ m/s/d2. O23 use
orvara (Brandt et al., 2021) to translate their reported RV trend and
curvature, the lack of an astrometric detection with Hipparcos (Lindegren et
al., 1997) and _Gaia_ , constraints from HRI, the lack of secondary lines in
their high-resolution spectra of HIP 9618, and the assumed stability of the
inner transiting planet system into orbital separation and mass ratio
posteriors for a distant, massive, single companion. O23 suggest the companion
is either a brown dwarf or low-mass M dwarf with $0.08^{+0.12}_{-0.05}$
$M_{\odot}$ in an orbit at $26^{+19}_{-11}$ AU.

The curvature reported by O23 is driven by just three SOPHIE RVs acquired
between 2011-Oct and 2011-Dec—the next SOPHIE RV (which also happens to be the
next RV from any of their three RV instruments) was taken eight years later in
2019-Dec. If we exclude the three SOPHIE RVs collected in 2011 from our RV
analysis (either with our joint model or with RadVel), then the curvature
detection disappears. Perruchot et al. (2008) quote SOPHIE’s RV stability as
being $\approx 3$ m/s over several months. To explore the possibility that the