# Predictions of Planet Detections with Near-infrared Radial Velocities in the Upcoming SPIRou Legacy Survey-planet Search

Ryan Cloutier, Etienne Artigau, Xavier Delfosse, Lison Malo, Claire Moutou, René Doyon, Jean-François Donati, Andrew Cumming, Xavier Dumusque, Elodie E. Hebrard, et al.

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

To cite this version:

Ryan Cloutier, Etienne Artigau, Xavier Delfosse, Lison Malo, Claire Moutou, et al. Predictions of Planet Detections with Near-infrared Radial Velocities in the Upcoming SPIRou Legacy Survey-planet Search. The Astronomical Journal, 2018, 155 (2), pp.93. 10.3847/1538-3881/aaa54e hal-02109411

HAL Id: hal-02109411

https://hal.science/hal-02109411v1

Submitted on 20 Aug 2025

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

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

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

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

# Predictions of Planet Detections with Near-infrared Radial Velocities in the Upcoming SPIRou Legacy Survey-planet Search

Ryan Cloutier<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> Ⓓ, Étienne Artigau<fnref n="3" />, Xavier Delfosse<fnref n="4" />, Lison Malo<fnref n="3" />,<fnref n="5" />, Claire Moutou<fnref n="5" />,<fnref n="6" /> Ⓓ, René Doyon<fnref n="3" />, Jean-François Donati<fnref n="7" />,<fnref n="8" />, Andrew Cumming<fnref n="9" /> Ⓓ, Xavier Dumusque<fnref n="10" /> Ⓓ, Élodie Hébrard<fnref n="11" />, and Kristen Menou<fnref n="1" />,<fnref n="2" />

Received 2017 October 23; revised 2017 December 4; accepted 2017 December 15; published 2018 January 30

## Abstract

The SPIRou near-infrared spectropolarimeter is destined to begin science operations at the Canada-France-Hawaii Telescope in mid-2018. One of the instrument's primary science goals is to discover the closest exoplanets to the solar system by conducting a three- to five-year long radial velocity survey of nearby M dwarfs at an expected precision of ~1 m s⁻<fnref n="1" />, the SPIRou Legacy Survey-Planet Search (SLS-PS). In this study, we conduct a detailed Monte Carlo simulation of the SLS-PS using our current understanding of the occurrence rate of M dwarf planetary systems and physical models of stellar activity. From simultaneous modeling of planetary signals and activity, we predict the population of planets to be detected in the SLS-PS. With our fiducial survey strategy and expected instrument performance over a nominal survey length of ~3 years, we expect SPIRou to detect 85.3⁺²⁹.<fnref n="3" />⁻¹².<fnref n="4" /> planets including 20.0⁺¹⁶.<fnref n="8" />⁻<fnref n="7" />.<fnref n="2" /> habitable-zone planets and 8.1⁺<fnref n="7" />.<fnref n="6" />⁻<fnref n="3" />.<fnref n="2" /> Earth-like planets from a sample of 100 M1–M8.5 dwarfs out to 11 pc. By studying mid-to-late M dwarfs previously inaccessible to existing optical velocimeters, SPIRou will put meaningful constraints on the occurrence rate of planets around those stars including the value of η⊕ at an expected level of precision of ≲45%. We also predict that a subset of 46.7⁺¹⁶.⁰⁻<fnref n="6" />.⁰ planets may be accessible with dedicated high-contrast imagers on the next generation of extremely large telescopes including 4.9⁺<fnref n="4" />.<fnref n="7" />⁻<fnref n="2" />.⁰ potentially imagable Earth-like planets. Lastly, we compare the results of our fiducial survey strategy to other foreseeable survey versions to quantify which strategy is optimized to reach the SLS-PS science goals. The results of our simulations are made available to the community on GitHub (https://github.com/r-cloutier/SLSPS_Simulations).

**Key words:** methods: statistical – planets and satellites: detection – planets and satellites: general – stars: low-mass techniques: radial velocities

**Supporting material:** Interactive figures

## 1. Introduction

The radial velocity (RV) method of detecting exoplanets is one of the most successful methods of exoplanet detection and has been widely used since the first discovery of an exoplanet around a main-sequence star over two decades ago (Mayor & Queloz 1995). Since then, numerous international teams have successfully built and used precision velocimeters to grow the population of RV planets (e.g., HARPS, Mayor et al. 2003; HARPS-N, Cosentino et al. 2012; HIRES, Vogt et al. 1994). Up until recently, the majority of these precision velocimeters have operated in the visible wavelength regime where their sensitivity is maximized for the discovery of planets around Sun-like stars.

In recent years, much interest has been generated regarding the population of exoplanets around M dwarfs. M dwarfs, with effective temperatures ≲3800 K and masses ≲0.6 M⊕, outnumber Sun-like stars in the solar neighborhood (within ~10 pc) nearly 4:1. Furthermore, M dwarfs are known to frequently host multiple

small ($r_p \le 4 R_⊕$) planets (e.g., Dressing & Charbonneau 2015; Gaidos et al. 2016) including a large fraction of planets within the star's habitable zone (HZ), which itself spans shorter orbital periods than that around more luminous Sun-like stars (Kasting et al. 1993; Kopparapu et al. 2013). Lastly, the amplitude of the RV signal induced by a given planet is larger around M dwarfs than that around Sun-like stars owing to their smaller masses. These favorable qualities have many astronomers committed to uncovering the M dwarf exoplanet population with purpose-built transit (e.g., MEarth, Irwin et al. 2015; ExTrA, Bonfils et al. 2015; TRAPPIST, Gillon et al. 2011; SPECULOOS, Gillon et al. 2013) and RV instrumentation (e.g., SPIRou, Delfosse et al. 2013b, Artigau et al. 2014; NIRPS, Bouchy et al. 2017; CARMENES, Quirrenbach et al. 2014; HPF, Mahadevan et al. 2012; IRD, Tamura et al. 2012).

One particular precision velocimeter optimized for the detection of exoplanets around M dwarfs in RV is SPIRou (Un Spectro-Polarimètre Infra-Rouge; Delfosse et al. 2013b; Artigau et al. 2014). SPIRou is a high-resolution near-infrared

<footnote n="1">¹Department of Astronomy & Astrophysics, University of Toronto. 50 St. George Street, Toronto, Ontario, M5S 3H4, Canada</footnote>

<footnote n="2">²Centre for Planetary Sciences, Department of Physical & Environmental Sciences, University of Toronto Scarborough.</footnote>

<footnote n="1265">1265 Military Trail, Toronto, Ontario, M1C 1A4, Canada</footnote>

<footnote n="3">³Institut de recherche sur les exoplanètes, Département de physique, Université de Montréal.
2900 boul. douard-Montpetit, Montréal, Quebec, H3T 1J4, Canada</footnote>

<footnote n="4">⁴Université Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="5">⁵CFHT Corporation, 65-1238 Mamalahoa Hwy, Kamuela, HI 96743, USA</footnote>

<footnote n="6">⁶Aix Marseille Université, CNRS, LAM, Laboratoire d'Astrophysique de Marseille, Marseille, France</footnote>

<footnote n="7">⁷Université de Toulouse, UPS-OMP, IRAP, 14 avenue E. Belin, Toulouse, F-31400 France</footnote>

<footnote n="8">⁸CNRS, IRAP/UMR 5277, Toulouse, 14 avenue E. Belin, F-31400, France</footnote>

<footnote n="9">⁹Department of Physics and McGill Space Institute, McGill University, 3600 rue des Universités, Montréal, Quebec H3A 2T8, Canada</footnote>

<footnote n="10">¹⁰Observatoire Astronomique de l'Université de Genève, 51 Chemin des Maillettes, 1290 Versoix, Switzerland</footnote>

<footnote n="11">¹¹Department of Physics and Astronomy, York University, Toronto, Ontario L3T 3R1, Canada</footnote>
velocimeter whose first light is scheduled on the Canada–France–Hawaii Telescope (CFHT) in 2018. A significant fraction of SPIRou’s allocated time will be spent surveying nearby M dwarfs searching for new exoplanets in a campaign known as the SPIRou Legacy Survey-Planet Search (SLS-PS). With the increased sensitivity to cool M dwarfs enabled by near-IR detectors, the planet detections resulting from the SLS-PS will be able to constrain the occurrence rate of planets around stars later than ~M4.5 and find the closest exoplanetary systems beyond Proxima Centauri (1.3 pc; Anglada-Escudé et al. 2016), which may be amenable to direct imaging with the next generation of imagers on an extremely large telescope (ELT).

In this study, we present a comprehensive simulation of the SLS-PS to estimate its planet yield as well as the bulk properties of the detected SPIRou planet population. These simulations were performed for a variety of survey strategies which enabled the SPIRou science team to establish an experimental setup that optimizes both the detection sensitivity and survey yield given the nominal time allocation of the SLS-PS. The main results of this study are based on the survey version deemed to be optimal for meeting the science goals of SPIRou. However, the results of the various surveys are summarized in Section 13. Comparison of the various survey versions may be useful to inform other upcoming RV planet searches similar to the SLS-PS.

The paper is organized as follows:

1. Section 2 gives an overview of the important aspects of the SPIRou spectropolarimeter.

2. Section 3 describes the stellar input catalog.

3. Section 4 describes the simulated SLS-PS.

4. Section 5 describes the population of simulated planetary systems.

5. Sections 6 and 7 describe how we mitigate the effects of stellar activity and detect planets.

6. Sections 8 and 9 describe the results of the survey.

7. Section 10 considers the effect of an increased planet frequency on the survey results.

8. Section 11 describes how well we can measure planet occurrence rates based on the results of the SLS-PS.

9. Section 12 discusses the potential for targeting SPIRou planets in direct imaging campaigns with ELTs,

10. Section 13 compares the merits of various potential versions of SLS-PS with the fiducial version presented throughout this paper.

## 2. Un Spectro-polarimètre Infra-rouge

SPIRou is an upcoming near-IR échelle spectropolarimeter and high-precision velocimeter whose first light is scheduled for 2018 on the Canada–France–Hawaii Telescope (CFHT) on Maunakea. The instrument is optimized to observe exoplanets via the RV technique around low-mass stars and to study the magnetic fields of young embedded protostars (Delfosse et al. 2013b). The design of SPIRou is intended to address its main science goals of detecting and characterizing M dwarf exoplanetary systems and to investigate the role that magnetic fields have in star and planet formation processes. SPIRou can be considered a heritage instrument built upon the success of the previous generation of optical spectropolarimeters and high-precision velocimeters such as the ESPaDOnS spectropolarimeter (Donati et al. 2006) as well as the RV

spectrographs SOPHIE (Bouchy & Sophie Team 2006) and HARPS (Mayor et al. 2003).

Here we provide a brief overview of the main instrument specifications as they pertain to the detection of new exoplanetary systems around nearby M dwarfs. Details of the optical and mechanical design of the instrument can be found in Artigau et al. (2014). The instrument itself is a fiber-fed, bench-mounted, double-pass, cross-dispersed spectropolarimeter that is cryogenically cooled to an operation temperature of 80 K and provides simultaneous spectroscopic and polarimetric observations. The optical fiber link connecting the Cassegrain unit—used for polarimetric analysis and guiding—to the calibration module and spectrograph is made from purified fluoride; a special optical material featuring improved transmission at wavelengths >2.0 µm enables the inclusion of the K band. The inclusion of the K band in the SPIRou spectral coverage is unique among most near-IR velocimeters and is highly desirable for the velocimetry of late M dwarfs as a large fraction of the RV information content is contained in the K band (É. Artigau et al. 2018, in preparation). With the broad continuous spectral coverage of SPIRou spanning the near-IR *YJHK* bands (0.98–2.35 µm), SPIRou will pioneer infrared planet searches by targeting low-mass stars with flux peaks in the near-IR wavelength domain. In order to detect Earth-size planets, SPIRou is required to achieve a long-term RV precision of 1 m s⁻¹ while simultaneously monitoring the star’s intrinsic activity, which is enabled by its high spectral resolution (λ/Δλ = 70,000).

Thanks to its spectropolarimetric capabilities, SPIRou is also optimized for characterizing stellar activity—an obvious asset for studying M dwarfs—and particularly late-type M dwarfs, which are known for their significant levels of magnetic activity (West et al. 2015). This will allow users to (i) minimize the impact of activity on RV curves and ease planet detections and (ii) to characterize the impact of stellar activity on the close-in HZ planets that SPIRou will detect.

## 3. Stellar Input Catalog

### 3.1. Stellar Sample

The SPIRou input catalog used in our simulated SLS-PS (see Section 4) contains 100 stars visible from CFHT on Maunakea ($\delta \gtrsim -30^{\circ}$). We note that the stellar sample used in these simulations is intended to be an approximation to the true SPIRou input catalog which has yet to be formalized exactly. The stars chosen were selected based on their high scores in the SPIRou merit function. The merit function is based on being able to detect the RV semi-amplitude from the gravitational pull of a 3 $M_{\oplus}$ planet at an equilibrium temperature of 250 K (L. Malo et al. 2018, in preparation). Of the optimal selection of 120 stars, approximately 20 were deemed to result in poor detection sensitivities based on the measured fractions of simulated planets detected around those stars in preliminary simulations of the SLS-PS. The primary culprit for the rejection of these stars was their large projected rotation velocities $v \sin i_s$, which have a large detrimental effect on the RV measurement precision and hence on our ability to detect planets.

The properties of stars in the SPIRou input catalog are presented in Figure 1. For comparison purposes, in Figure 1 we include six nearby M dwarf planetary systems with at least one known planet in or near the HZ: Proxima Centauri (Anglada-Escudé et al.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Scatter plots and histograms depicting the distribution of the SPIRou input catalog *J*-band magnitudes, stellar masses, distances, and effective temperatures for our fiducial version of the SPIRou Legacy Survey-Planet Search. Histograms are in linear units. Nearby M dwarf planetary systems with at least one known HZ planet are depicted with diamonds. K2-18 at 34 pc does not appear in the scatter plot in right panel.</figcaption></figure>

2016), Ross 128 (Bonfils et al. 2017), GJ 273 (Astudillo-Defru et al. 2017), LHS 1140 (Dittmann et al. 2017), TRAPPIST-1 (Gillon et al. 2017), and K2-18 (Montet et al. 2015; Cloutier et al. 2017a). The global properties of the 100 stars were defined from the Boyajian et al. (2012) relation (effective temperature and radii, when [Fe/H] is fixed to the solar value) and the Delfosse et al. (2000) relation (estimated stellar mass based on absolute *J* magnitude). Projected rotational velocities and rotation periods were taken from an exhaustive literature search or derived from the CFHT-CoolSnap program (Moutou et al. 2017). In our sample, we consider stars with masses between 0.08 and 0.57 *M*⊙ with *J*-band magnitudes of 4.2–10.3 in a range of distances spanning ~1.8–11 pc.

Rotational information for stars in our sample is important for RV modeling as stellar rotation strongly affects the activity arising from the rotationally modulated active regions (ARs) observed in RVs (Saar & Donahue 1997; Meunier et al. 2010; Aigrain et al. 2012; Dumusque et al. 2014). Rotation timescales also restrict the periodicities at which we can detect planets due to difficulties in detecting RV planets with orbital periods close to the stellar rotation period or its harmonics (Vanderburg et al. 2016). However, the rotational information for our stellar sample is incomplete. For stars with no known available rotation measurements (8/100 stars) or with a $v \sin i_s$ upper limit only (2/100 stars), we sample $P_{\text{rot}}$ from a modified empirical distribution of M dwarf rotation periods from ground-based photometry as a function of stellar mass (Newton et al. 2016). The nature of this “modification” is discussed in the subsequent paragraph. The corresponding $v \sin i_s$ is then computed from the sampled value of $P_{\text{rot}}$, the known stellar radius, and the inclination of the stellar spin axis to the line-of-sight $i_s$, which we draw from a geometrical distribution (i.e., uniform in $\cos i_s$). For stars that only have a measured upper limit on $v \sin i_s$, the modified empirical distribution from which $P_{\text{rot}}$ is sampled is truncated at the minimum $P_{\text{rot}}$ corresponding to the upper limit on $v \sin i_s$.

The necessary modification to the empirical $P_{\text{rot}}$ distribution arises from an observational bias in the Newton et al. (2016) sample, which favors rapidly rotating stars. The detection of short photometric rotation periods (i.e., rapid rotators) is attained more easily than long rotation periods because full phase coverage is more readily obtained over many rotation cycles. Furthermore, in the case of early M dwarfs ($M_s > 0.25 \ M_{\odot}$), there is evidence for a positive correlation between the star’s rotation rate and the amplitude of its photometric variability (Newton et al. 2016). Rapid rotators therefore tend to exhibit larger amplitudes of variability, thus making the signal more easily detectable. The raw empirical distribution therefore does not represent the true underlying distribution of M dwarf rotation periods in the solar neighborhood. We attempt to account for this bias in a simplified way by modifying the empirical $P_{\text{rot}}$ distribution by insisting that only ~25% of sampled rotation periods can be < 10 days, as estimated from the volume-limited sample of field M dwarfs with measured rotation velocities from Delfosse et al. (1998). This constraint reduces the fraction of fast rotators with $P_{\text{rot}}$ < 10 days by a factor of ~2.

The modification to the empirical $P_{\text{rot}}$ distribution is visualized in Figure 2. Here we compare the empirical $P_{\text{rot}}$ distribution, based on the full stellar sample with detected $P_{\text{rot}}$ from Newton et al. (2016), with the modified distributions of stars in the SPIRou input catalog with either a $v \sin i_s$ upper limit only or no available rotation data. The latter two distributions are nearly equivalent as the majority of $v \sin i_s$ upper limits do not provide substantial new information regarding the star’s $P_{\text{rot}}$, and in both cases we insist that only ~25% of sampled $P_{\text{rot}}$ can be < 10 days. We impose this condition by noting that ~60% of stars in the empirical distribution have $P_{\text{rot}}$ < 10 days and resample a particular fraction of those stars from the subset of the empirical distribution restricted to $P_{\text{rot}} \ge 10$ days. The fraction of stars with $P_{\text{rot}}$ < 10 days that get resampled is 1 - 0.25/0.6 ≈ 0.58.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Cumulative distribution functions of the empirical distribution of M dwarf rotation periods from Newton et al. (2016), the modified empirical distribution of SPIRou stars with a v sin i_s upper limit measured, the modified empirical distribution of SPIRou stars with no rotation data available, and the full SLS-PS sample.</figcaption></figure>

## 3.2. Physical Models of Stellar RV Activity

### 3.2.1. SOAP 2.0: Activity Simulations

ARs in the stellar photosphere (e.g., star spots and faculae) and in the hot chromosphere (e.g., plages) are expected to be present in M dwarfs. These surface inhomogeneities have characteristic temperatures that differ from the star's effective temperature and therefore disrupt the symmetry of the visible stellar disk as they rotate in and out of view at the stellar rotation period. One resulting source of RV activity from ARs, known as the flux effect (Dumusque et al. 2014), results in an anomalous RV signal as the ARs block a fraction of Doppler-shifted photons from the rotating stellar limbs. The strong local magnetic fields associated with ARs at the stellar photospheric boundary also act to inhibit the upward flow of hot convective material in an effect known as the suppression of convective blueshift (Dravins et al. 1981).

In the limit of simple distributions of ARs, the observed RV structure from the two aforementioned effects is dependent on the fractional coverage of the visible stellar disk by the ARs and its first time derivative (Aigrain et al. 2012). For each simulated RV time series, we sample the relevant physical parameters of the ARs (i.e., AR sizes and spatial distribution) and simulate the corresponding RV activity, FWHM, bi-sector inverse slope (BIS), and photometric time series arising from both the flux effect and from the suppression of convective blueshift using the SOAP 2.0 code (Dumusque et al. 2014). The FWHM and BIS time series are shape parameters of the cross-correlation function between the observed stellar spectra and the template spectrum used to measure the stellar RVs in an observing campaign. These ancillary time series are sensitive to the presence of ARs but not to planets, making them useful diagnostics for distinguishing activity-induced RV signals from planetary signals. In particular, the FWHM time series will be used in Section 6 to train our RV activity model and disentangle RV activity signals from planetary signals.

The SOAP 2.0 code outputs time series that are phase-folded to the input stellar rotation period. These time series are initially treated as strictly periodic and interpolated to the

epochs of observation. In this way we ignore any contribution from differential rotation, the amplitude of which has been shown to decrease with decreasing stellar mass (Donati et al. 2008; Morin et al. 2008; Kitchatinov & Olemskoy 2011), before invoking rigid-body rotation in fully convective M dwarfs ($M_s \lesssim 0.2 M_\odot$). In Section 3.2.3, the strictly periodic condition is relaxed to account for the finite lifetimes of ARs and the existence of long-term magnetic activity cycles.

Unfortunately, very little is presently known about the physical nature of ARs on M dwarfs but we do observe quasi-periodic photometric variability which arises from evolving ARs (O'Neal et al. 2005). We use the empirical distribution of photometric amplitudes as a function of stellar mass from Newton et al. (2016) to sample the average photometric variability amplitudes<fnref n="12" /> A for each star in the simulated SLS-PS. The sampled value of A is interpreted as an average value because the phase in the star's magnetic activity cycle at the time of the Newton et al. (2016) observations is unknown. We then use A to constrain the size of up to four ARs, with each AR being treated as either a cool spot or a bright plage. The spatial distribution of ARs is determined from random draws in latitude and longitude as Doppler imaging provides evidence for more uniformly distributed ARs on M dwarfs than on Sun-like (FGK) stars (Barnes & Collier Cameron 2001; Barnes et al. 2004), the ARs of which tend to be more localized along the stellar equator. We note, however, that such Doppler imaging observations are limited to rapid rotators that are largely avoided in the selection of the SPIRou input catalog. The temperature contrast between the ARs and the stellar effective temperature is fixed to 200 K in all realizations (Berdyugina 2005). The SOAP 2.0 code is designed to model Sun-like stars at an optical wavelength of $\lambda \sim 529$ nm. The resulting activity is then scaled from the default SOAP 2.0 wavelength to the approximate central H-band wavelength of $\lambda' \sim 1.6 \mu$m via the ratio of blackbody emission at $\lambda$ to $\lambda'$ with a characteristic temperature of $T_{eff}$. This scaling decreases the amplitude of the flux and convective blueshift effects by a typical factor of a few in the near-IR compared to that at optical wavelengths (Martín et al. 2006; Huélamo et al. 2008; Prato et al. 2008; Reiners et al. 2010; Mahmud et al. 2011).

### 3.2.2. Zeeman Broadening

Unlike the flux and convective blueshift effects, the RV activity due to Zeeman broadening tends to increase toward the near-IR. The Zeeman broadening of spectral features in unpolarized light occurs in the presence of strong magnetic fields that cause Zeeman splitting, an effect that grows with wavelength. Reiners et al. (2013) and Hébrard et al. (2014) used polarized radiative transfer at near-IR wavelengths to compute the effect of Zeeman splitting from both atomic and molecular sources on stellar line profiles. Reiners et al. (2013) report the following simplified model for the RV signal resulting from Zeeman broadening in M dwarfs ($T_{eff} \in [2800, 3700]$ K):

$$RV_z(t) = 300 \text{ m s}^{-1} f(t) \left(\frac{B}{1 \text{ kG}}\right)^2 \left(\frac{\lambda}{1 \mu\text{m}}\right)^a, \quad (1)$$

<footnote n="12">¹² For photometric variations measured in the near-IR with MEarth over a custom passband spanning ~0.7–1μm known as the i + z band (Nutzman & Charbonneau 2008).</footnote>
where *f* is the filling factor or fraction of the visible stellar disk that is spanned by ARs, *B* is the local magnetic field strength within the AR, and $\lambda$ is the wavelength of observation. The power-law index $a \in [0, 2]$ describes the increase of RVz with $\lambda$ and is variable as a result of the apparent distribution of molecular Landé g-values in cool stars. However, only the FeH and CO bands are considered in the stellar atmospheric model from which Equation (1) is derived (Reiners et al. 2013).

Computing RVz to add to our complete physical RV activity model requires knowledge of the local *B* field strength within ARs. The empirical distribution of this quantity in M dwarfs is incomplete despite contributions from various observing campaigns (e.g., Reiners & Basri 2007; Shulyak et al. 2014; Hébrard et al. 2016; Moutou et al. 2017; Shulyak et al. 2017). Instead of sampling *Bf* from an empirical distribution, we use an ad hoc method of sampling *Bf* that exploits what is known about small-scale *Bf* fields in M dwarfs as a function of spectral type and rotation, namely, the fraction of M dwarfs that are magnetically active as a function of rotation period differs between early-type and late-type M dwarfs (M5–M8) as later M dwarfs are able to remain magnetically active late into their lives even after considerable spin-down (West et al. 2015). For each star, we assign an activity flag indicating whether it is an active or inactive star, where the probability of being flagged as an active star is equal to the measured activity fraction of M dwarfs from West et al. (2015) and is dependent on the star’s spectral type and $P_{\text{rot}}$. If magnetically inactive, we sample the localized magnetic field strength from $B \sim U(1, 1)$ kG (Moutou et al. 2017). If magnetically active, we instead draw from $B \sim U(1, 3.1)$ kG (Moutou et al. 2017). We then calculate the value of *a* based on the star’s sampled *B* and spectral type before evaluating the Zeeman broadening model (Equation (1)) at the central *H*-band wavelength of 1.6 µm as a function of the time-evolving filling fraction, which is known from our SOAP 2.0 simulations.

### 3.2.3. Active Region Lifetimes

The ARs giving rise to stellar activity in our simulations are short lived compared to the baseline of our observations. RV observations of M dwarfs have suggested that the lifetimes of individual ARs may persist from one to a few stellar rotations and up to $\geq 10$ (e.g., Bonfils et al. 2007; Forveille et al. 2009; Hébrard et al. 2016). RV observations have also elucidated that M dwarfs undergo long-term magnetic activity cycles similarly to the Sun (e.g., Gomes da Silva et al. 2012; Route 2016). Following the prescription of Dumusque (2016) for Sun-like stars, we proceed in deriving the temporal variation of AR sizes by scaling the total RV activity signal according to each AR’s appearance rate $\lambda(t)$. Dumusque (2016) also included the time-dependent latitude of ARs which we ignore here due to the more uniform distribution of ARs observed on M dwarfs (Barnes & Collier Cameron 2001; Barnes et al. 2004). Furthermore, we assume that the appearance rate for both star spots and bright plages are consistent.

The probability that an AR appears at a time *t* is governed by the Poisson distribution

$$P(t) = \frac{e^{-\lambda(t)\tau}(\lambda(t)\tau)^k}{k!}, \quad (2)$$

where $\tau$ is the time step in days and $k = 0, 1, 2, 3$ as we only consider a maximum of four ARs. Next, as a function of time,

we draw from the probability distribution in Equation (2), which dictates at which epochs an AR is formed. For each newly formed AR, we insist that it spends the first third of its lifetime evolving linearly to its maximum size before shrinking toward zero over the remaining two-thirds (Dumusque 2016). Each AR’s lifetime is sampled from a truncated Gaussian distribution with mean $3P_{\text{rot}}$ and standard deviation $P_{\text{rot}}$. The Gaussian distribution is truncated at $P_{\text{rot}}$ such that all sampled ARs persist for a minimum of one stellar rotation (e.g., Bonfils et al. 2007; Forveille et al. 2009; Hébrard et al. 2016). Also recall that the maximum size of the AR is determined by the star’s sampled amplitude of photometric variability.

The AR appearance rate per unit time is

$$\lambda(t) = (\lambda_{\text{max,act}} - 0.5) \left[ -0.5 \cos \left( \frac{2\pi t}{P_{\text{cycle}}} + \phi \right) + 0.5 \right] + 0.5, \quad (3)$$

where $\lambda_{\text{max,act}}$ describes the maximum appearance rate during the maximum of the star’s magnetic activity cycle with period $P_{\text{cycle}}$. We set $\lambda_{\text{max,act}} = 10$ ARs per day and sample $P_{\text{cycle}}$ from $U(6, 10)$ years (Suárez Mascareño et al. 2016; Wargelin et al. 2017), which is a factor of two or more greater than the baseline of the observations. The added term of 0.5 ARs per day to Equation (3) ensures that we maintain a low but non-zero probability of forming an AR close to the minimum of the stellar activity cycle.

Recall that the sampled amplitude of photometric variability sets the size of ARs in our simulations. Furthermore, because the phase within a star’s magnetic activity cycle at the time of photometric observations is unknown, we treat the observed photometric variability amplitude as an average value. To account for this approximately, we rescale our derived AR lifetime scaling to the interval 0.1–2—instead of 0–1—to account for the varying levels of stellar activity up to a factor of two greater than the maximum value and down to a minimum value slightly greater than zero. We then use this scaling to rescale the injected activity in both RVs and the ancillary time series.

In general, our rescalings were optimized such that the root mean square (rms) of the injected RV activity is roughly consistent with the upper envelope of the RV activity rms observed with HARPS. Figure 3 depicts the distribution of RV activity rms in our simulated stellar sample as a function of $P_{\text{rot}}$ and directly compares it to the values of M dwarfs observed with HARPS. The HARPS M dwarfs with measured RV activity rms in Figure 3 include Proxima Centauri (Anglada-Escudé et al. 2016), GJ 3293, GJ 3341, GJ 3542 (Astudillo-Defru et al. 2015), GJ 1132 (Berta-Thompson et al. 2015), GJ 876 (Correia et al. 2010), GJ 674 (Bonfils et al. 2007), Gl 205, Gl 358, Gl 388, Gl 479, Gl 526, Gl 846 (Bonfils et al. 2013a), GJ 163 (Bonfils et al. 2013b), Gl 433, Gl 667C (Delfosse et al. 2013a), Gl 176 (Forveille et al. 2009), GJ 205, GJ 358, GJ 410, GJ 479, GJ 846 (Hébrard et al. 2016), and GJ 436 (Lanotte et al. 2014).

We note that the comparison depicted in Figure 3 is not one to one as RV activity is an intrinsically chromatic effect, and our simulated time series are computed at near-IR wavelengths whereas the HARPS measurements are taken in the visible. The effect of temperature contrast between ARs and the stellar photosphere on observed RV activity is known to
<figure><img src="image_6.png" /><figcaption>**Figure 3.** RV activity rms as a function of stellar rotation period for stars in the simulated SLS-PS compared to a population of M dwarfs observed with the HARPS spectrograph. The dashed white curve depicts the mean RV activity rms in each $P_{\text{rot}}$ bin.</figcaption></figure>

decrease toward longer wavelengths. Conversely, activity from Zeeman broadening increases toward longer wavelengths and is known to be an important source of activity in the near-IR (Hébrard et al. 2014; Moutou et al. 2017). Our adopted activity scaling is chosen to best match the HARPS observations at the most frequently sampled rotation periods in our stellar sample of ~50–120 days. In Figure 3, it is clear that the mean RV activity rms of our sample closely matches that of the HARPS stars, albeit with a large dispersion. However, at smaller $P_{\text{rot}}$ (≤5 days), the mean RV activity rms in our simulated sample becomes slightly underestimated relative to the small sample of HARPS stars at those rotation periods. A larger sample of M dwarfs with measured RV activity rms at small $P_{\text{rot}}$ will be required to precisely characterize the correlation between the activity rms and $P_{\text{rot}}$.

### 3.2.4. Ignored Sources of RV Activity

Various other minor sources of stellar activity in M dwarfs are ignored. These include radial pulsations, which do not appear to persist with any significant amplitude within the interiors of M dwarfs (Rodríguez-López et al. 2015). Second, we disregard the effect of granulation, the amplitude of which scales with the velocity of convective cells which itself decreases toward later spectral types (Dumusque et al. 2011; Meunier et al. 2017). Lastly, we ignore flares in our activity model because their distinctive spectral signature allows them to be flagged and removed from the subsequent analysis (Schmidt et al. 2012; Anglada-Escudé et al. 2016). However, the high occurrence rate of flaring events on many M dwarfs may prove costly when attempting to construct large time series of uncontaminated RVs.

## 4. Simulated Survey

In this study, we conduct a detailed Monte Carlo (MC) simulation of the SLS-PS for the purpose of predicting the SPIRou planet detection yield. To do so, we must construct a statistically significant number of unique RV time series for each of the 100 stars that we target throughout the simulated SLS-PS. In practice, we simulate 100 unique RV time series for each star, totalling $10^4$ realizations in the full simulated SLS-PS. Each RV time series contains unique signals from planets, which are sampled from their known occurrence rates around early M dwarfs (see Section 5), from physical models of stellar activity (as described in Section 3.2), and from instrumental noise. Each contribution is sampled in time using a unique window function spanning ~300 nights over ~3 years, which represents the typical time baseline of observations for a single star in accordance with the expected subset of SPIRou's time allocation that will be dedicated to the discovery of new exoplanetary systems in the SLS-PS. We note, however, that at this time, neither the exact number of available nights dedicated to the SLS-PS nor the duration of the full SLS-PS have been established absolutely. In addition to the RVs, we also derive various spectroscopic activity indicators arising from stellar activity that are contemporaneous with the RVs. One such ancillary time series is the FWHM of the cross-correlation function which will be used to train non-parametric Gaussian process models of RV activity based on its common covariance structure with the FWHM (see Section 6). SPIRou is also unique in that it simultaneously operates its spectroscopic and polarimetry modes, thus providing a contemporaneous diagnostic of the star's magnetic topology. Such time series may also be used to model RV signals from magnetically ARs (Hébrard et al. 2016), although we do not consider such time series in this study.

In each simulated RV time series, we attempt to recover the injected planets to form an estimate of the expected planet population that will be discovered with SPIRou. This is facilitated by the joint modeling of stellar activity and planetary signals. This allows for the self-consistent characterization of each RV signal and the detection of a subset of planets that are nominally hidden by their host star's intrinsic RV activity.

### 4.1. Window Functions

In each MC realization, the unique window function $t$ is the vector of length $n_{\text{obs}}$ containing the epochs of observation in barycentric Julian dates. The window function describes the time sampling of our time series. To derive a set of window functions for each star in the SPIRou input catalog, we run a separate MC simulation of stellar observing sequences for all targets, taking into account when each star is visible from CFHT on Maunakea with an airmass of < 2.5 based on its celestial coordinates. During every available night, all visible stars are observed up to two times, each with an integration time required to achieve a signal-to-noise ratio (S/N) per resolution element of 150—at the central J-band wavelength of 1.25 µm—with a minimum integration time of 15 minutes. The imposed lower limit on the integration time may also be necessary to mitigate the effects of granulation (Lovis et al. 2005), which is expected to be low on M dwarfs. Integration times required to achieve at least the target S/N are computed for each star based on using their *YJHK* magnitudes. An overhead of 5 minutes is added to each integration for
<figure><img src="image_7.png" /><figcaption>Figure 4. Lomb–Scargle periodogram of an example window function from the simulated SLS-PS. The periodogram exhibits strong aliases from our time sampling at the periodicities highlighted by the vertical lines. The power corresponding to a 1% false-alarm probability in the periodogram is highlighted by the horizontal dashed line.</figcaption></figure>

guiding and setup purposes. The output from these simulations is a set of window functions, each pertaining to a star in the SPIRou input catalog. Multiple MC simulations are run for various observing sequences and thus provide unique window functions for the simulated SLS-PS. We then sample from these window functions for the purpose of investigating the sensitivity of our planet detection results to the exact form of the window function.

The available nights for observation with SPIRou are limited by two important considerations, the first being the effect of stochastic weather, which limits the number of epochs in our derived window functions according to the CFHT observatory's weather statistics. The second effect is somewhat unique to CFHT as the telescope hosts a suite of instruments that do not operate simultaneously. In particular, the telescope's wide-field optical imager requires dark time to conduct its observations whereas the SPIRou spectropolarimeter does not. As such, we proceed with constructing SPIRou window functions that only include non-dark-time observation and thus correspond to higher levels of lunar contamination. This represents a worst case scenario for SPIRou as aliases from the window function will undoubtedly arise at periodicities close to the cadence of the non-dark-time observing sequences. This cadence evolves with a period close to the period of the lunar cycle at ~30 days. Because the SPIRou time sampling occurs in windows separated by ~30 days, at that period and its first harmonic at 15 days, significant aliases in the Lomb–Scargle (LS) periodogram (Scargle 1982) of the window function can arise as is shown in Figure 4. The effect of these aliases are detrimental to the detection of periodic planetary signals at these periods because in the LS periodogram of the RVs, one cannot distinguish a priori between these periodic signals as a planet or as an alias of the time sampling (Dawson & Fabrycky 2010). This effect has already been shown to mimic planetary signals (e.g., Rajpaul et al. 2016) and is particularly detrimental to finding HZ planets around ~M2–M4 dwarfs with HZ spanning ~30 day orbital periods. We quantify the magnitude of this aliasing effect on the SPIRou planet detection sensitivity in Section 8.

For the purpose of performing accurate statistical inference of the SPIRou planet population following the full SLS-PS (see Section 11.1), we attempt to split the total observing time

evenly among the stars in the SPIRou input catalog. Doing so, while taking into account stochastic weather and non-dark-time restrictions, limits the size of our window functions over an a priori three-year-long survey to an average of 198.1 RV measurements per star for 100 stars over ~300 nights. The variance in the number of RV measurements per star is relatively small and ranges from 181 to 212. In the other versions of the SLS-PS containing either more or fewer stars than in our fiducial survey version (see Section 13), our MC calculations of the window functions consequently contain fewer and more RV measurements per star, respectively.

From preliminary simulations of the SLS-PS using 10 unique window functions per sampled planetary system per star, we found that the net planet detection results are largely independent of the exact window function used. Note that all sampled window functions contained the same restrictions discussed previously but do not include the exact same epochs of observation for each star. Explicitly, the SPIRou planet yield was found to vary by only ~1% across the various window functions used. This dispersion is much less than the uncertainties on the resulting estimates of the SPIRou planet yield due solely to uncertainties in the input planet occurrence rates. Therefore, we conclude that using multiple window functions for each planetary system in our simulated SLS-PS is an unnecessary computational expense that can be mitigated by considering a single window function per planetary system and will not significantly affect the results of our study. However, each unique window function considered was still assumed to be restricted to non-dark-time observations, thus preserving the aliasing effect on planets with orbital periods of ~30 days.

## 4.2. RV Time Series Construction

In each MC realization, the RV contribution from $N_p$ injected planets in the simulated planetary system (see Section 5) is calculated via the superposition of $N_p$ Keplerian orbital solutions. Each Keplerian RV$_{kep}$ is parameterized by the planet's orbital period $P$, time of inferior conjunction $T_0$, orbital eccentricity $e$, argument of periastron $\omega$, and RV semi-amplitude $K$ according to

$$RV_{\text{kep}}(t) = K[\cos(\nu(t) + \omega) + e \cos \omega], \quad (4)$$

where $\nu(t)$ is the true anomaly and is computed by solving Kepler's equation and the eccentric anomaly as a function of time $t$ contained in the window function $t$. The RV semi-amplitude $K$ is computed from the planet's minimum mass $m_p \sin i$, the stellar mass $M_s$, $P$, and $e$ using the standard formula

$$K = 1.05 \text{ m s}^{-1} \left( \frac{m_p \sin i}{2 M_{\oplus}} \right) \left( \frac{P}{20 \text{ days}} \right)^{-1/3} \times \left( \frac{M_s}{0.3 M_{\odot}} \right)^{-2/3} \frac{1}{\sqrt{1 - e^2}}. \quad (5)$$

The Keplerian RV approximation is valid in all single-planet systems and the majority of multiplanet systems considered and is best motivated by its ability to negate the need to perform costly numerical integrations of each sampled planetary system. However, we note that this approximation naturally excludes certain dynamical effects in multiplanet systems such as planet–planet interactions and mean-motion resonances,
which can affect the planet-induced periodicities within the RV time series and therefore also affect our ability to detect those planets. Although the former effect is not accounted for in our simulated time series, its amplitude in real systems is typically small compared to the RV measurement uncertainty as a result of dynamical restrictions on multiplanetary systems, making tightly packed systems less stable over long timescales and therefore rarely seen in nature. The significance of the latter effect is also expected to be small given the dearth of multiplanet systems at low-order period ratio commensurabilities (Lissauer et al. 2011; Fabrycky et al. 2014).

All RV components are evaluated at *t* and contain the additive (i) RV activity $RV_{act}$ derived from physical models (see Section 3.2), (ii) $N_p$ Keplerian models, and (iii) a white noise term with an rms equal to the median RV measurement uncertainty $\sigma_{RV}$ expected for the host star. Each star's $\sigma_{RV}$ value is calculated using the near-IR RV information content calculations from Figueira et al. (2016) corrected using the empirical spectra from É. Artigau et al. (2018, in preparation). We use the results from condition 3 in Figueira et al. (2016) to estimate $\sigma_{RV}$: the photon-noise contribution to the spectrum being amplified by the limited spectral window due to atmospheric telluric transmission. The information content calculations for each star are computed using SPIRou's spectral resolution, spectral coverage, and the expected S/N per resolution element obtained during an integration of the star.

The RV activity model contains additive contributions from the flux effect, the suppression of convective blueshift, and Zeeman broadening from ARs. The complete RV model is therefore

$$RV_{\text{model}}(t) = RV_{\text{act}}(t) + \sum_{i=1}^{N_p} RV_{\text{kep},i}(t) + \mathcal{N}(0, \sigma_{RV}). \quad (6)$$

The RV measurement uncertainties pertaining to each measured RV is contained in the vector $\sigma_{RV}$ and is modified from the scalar value of $\sigma_{RV}$ based on the variable absorption by terrestrial water vapor. The water vapor correction at each epoch in *t* is the product of the airmass and the zenith water column from the well-documented CFHT observing condition statistics throughout the calendar year. The distribution of the median $\sigma_{RV}$ in our MC simulation for our full stellar sample is shown in Figure 5 as a function of the stellar J-band magnitude and $v \sin i_s$. A noise floor is imposed at the expected long-term RV precision limit of SPIRou at $1 \text{ m s}^{-1}$. This results in a median $\sigma_{RV} = 1.52 \text{ m s}^{-1}$. The maximum $\sigma_{RV}$ is 6.58 m s$^{-1}$. Not depicted in Figure 5 is the dependence of $\sigma_{RV}$ on spectral type as the RV measurement uncertainty tends to decrease toward later spectral types due to the increased number of available spectral features and corresponding increase in RV information content.

Similarly to the construction of the RV activity time series at the epochs in *t*, the FWHM and BIS ancillary time series are constructed from our SOAP 2.0 simulations (see Section 3.2.1) at the same epochs as the RVs. These contemporaneous time series will be used in Section 6 to model the RV activity and help to detect underlying planetary signals.

<figure><img src="image_8.png" /><figcaption>Figure 5. Median RV measurement uncertainty in each simulated RV time series $\sigma_{RV}$ as a function of the star's J-band magnitude and projected stellar rotation velocity $v \sin i_s$. A noise floor is imposed at $1 \text{ m s}^{-1}$. The histogram shown in the right panel is in linear units.</figcaption></figure>

# 5. Planet Sample

## 5.1. Planet Sample from Known Occurrence Rates

In each MC realization, we populate the simulated M dwarf planetary system with planets according to their known occurrence rates from the Kepler transit survey (Dressing & Charbonneau 2015). In practice, we sample from a grid in orbital period *P* and planetary radius $r_p$ according to the measured planetary occurrence rate at each point in the parameter space *f*(*P*, $r_p$). We assume that each point in the parameter space is uncorrelated aside from the dynamical stability constraints that we will apply to closely packed systems (see Section 5.2).

The occurrence rates $f(P, r_p)$ (or equivalently, $f$) from which we sample planets have associated asymmetric uncertainties which are derived from data-driven posterior distributions rather than from analytically defined distributions which are more easily sampled from. Dressing & Charbonneau (2015) only present the modes and $1\sigma$ dispersions of their *f* posteriors, which is insufficient information to reconstruct their *f* posteriors in each *P* and $r_p$ bin. To account for their reported uncertainties on $f$, we must assume a functional form of the *f* posteriors. For each bin in *P* and $r_p$, we sample planets with a probability drawn from $N(\mu, \bar{\sigma})$, where $\mu$ is the most likely value of $f$ and $\bar{\sigma} = \text{mean}(\sigma_{\text{upper}}, \sigma_{\text{lower}})$, where $\sigma_{\text{upper}}$ and $\sigma_{\text{lower}}$ represent the asymmetric uncertainties on $f$. The values of $\mu$, $\sigma_{\text{upper}}$, and $\sigma_{\text{lower}}$ are provided in Dressing & Charbonneau (2015; see their Table 4) for all *P* and $r_p$ bins for which $f > 0$. For bins in which $f$ is consistent with zero, we set $\mu = 0$ and $\bar{\sigma}$ equal to the $1\sigma$ upper limit on $f$. For bins in which $f$ is completely unconstrained as a result of a low detection sensitivity, we assume that $f(P, r_p)$ evolves smoothly such that we can extrapolate the values of $f$ from surrounding bins with constraints in order to estimate the values of $f$ where it is unconstrained by the data (e.g., $r_p \le 1 R_\oplus$ and $P \ge 18.2$ days). The extrapolated values are treated as $1\sigma$ upper limits, and $\mu$ is set to zero such that when integrating the most likely values of $f(P, r_p)$ over the range of periods and
planet radii considered by Dressing & Charbonneau (2015; $P \in [0.5, 200]$ days, $r_p \in [0.5, 4]$ $R_⊕$) we recover their cumulative planet occurrence rate of 2.5 ± 0.2 planets per M dwarf. However, our dynamical stability constraints will reject a subset of sampled planetary systems, thus effectively reducing the cumulative planet occurrence rate to < 2.5.

The dispersion in mutual inclinations among planets in multiplanet systems is related to the dispersion in eccentricities by $\langle i^2 \rangle \sim \langle e^2 \rangle/4$ (Stewart & Ida 2000; Quillen et al. 2007) with the mean orbital inclination to the plane of the sky $i$ being set to the value of $i_s$ obtained for the host star (i.e., from a geometrical distribution). Although the distribution of spin-orbit angles for small planets around M dwarfs has yet to be established, the observed low dispersion in orbital eccentricities (Van Eylen & Albrecht 2015) and mutual inclinations (Figueira et al. 2012; Fabrycky et al. 2014) in these types of planetary systems suggests that they are dynamically cold. If this is indeed the case, then the normal vector to the mean planetary orbital plane is expected to be close to parallel to the stellar spin axis, $i \approx i_s$. In each planetary system, we sample each planet's orbital eccentricity $e$ from the $\beta$ probability distribution describing the high detection significance sample of RV planets reported in Cloutier et al. (2015; see also Kipping 2013).

## 5.2. Modifications to the Planet Sample

Our adopted approach for sampling planetary systems is accompanied by four important caveats.

Converting planet radii to masses. First, the Kepler-derived $f$ directly samples planetary radii whereas RV surveys are only sensitive to the planetary minimum masses. Therefore, an assumption must be made regarding the planetary mass-radius relation required to convert the sampled planetary radii into masses. We opt for the following empirical mass-radius relation derived in Weiss & Marcy (2014) from an unbiased sample of known transiting planets:

$$ \frac{m_p}{M_\oplus} = \begin{cases} 0.440 \left(\frac{r_p}{R_\oplus}\right)^3 + 0.614 \left(\frac{r_p}{R_\oplus}\right)^4, & r_p < 1.5 R_\oplus \\ 2.69 \left(\frac{r_p}{R_\oplus}\right)^{0.93}, & r_p \ge 1.5 R_\oplus. \end{cases} \quad (7) $$

This piecewise mass-radius relation distinguishes between small rocky planets with $r_p < 1.5 R_⊕$ and larger gaseous planets (e.g., Dressing et al. 2015; Rogers 2015; Fulton et al. 2017). The intrinsic dispersion about the mean mass-radius relation in Equation (7) has characteristic rms values of

$$ \sigma_{\text{rms}} = \begin{cases} 2.7 M_{\oplus}, & r_p < 1.5 R_{\oplus} \\ 4.7 M_{\oplus}, & r_p \ge 1.5 R_{\oplus}. \end{cases} \quad (8) $$

To include this intrinsic dispersion in our sampled planet population, for each sampled $r_p$, we compute the mean $m_p$ using Equation (7) and add an additional offset drawn from $N(0, \sigma_{\text{rms}})$, where the value of $\sigma_{\text{rms}}$ is given by Equation (8). We reject planets with unphysical negative masses, which naturally biases our sample to larger planet masses than are predicted by the mean mass-radius relation. We then apply a unique correction factor to each radius bin $r_p$ with a width of 0.2 $R_⊕$ such that we recover the mean mass-radius relation given in Equation (7) in our final planet sample. Our Kepler-derived planet sample in the mass-radius plane is shown in

<figure><img src="image_9.png" /><figcaption>**Figure 6.** 2D histogram depicting the masses and radii of the sampled small planet population with $r_p \le 4 R_⊕$. Planetary radii are drawn from the planet occurrence rates derived from Kepler and are mapped to planetary masses using the mass-radius relation and Gaussian dispersion given in Equations (7) and (8). The mean mass-radius relation is overplotted as the dashed curves. The mean planet mass and standard deviation of the mean in each planetary radius bin are plotted as black circles to show their close correspondence with the mean mass-radius relation.</figcaption></figure>

Figure 6 along with the mean mass-radius relation. Overplotted are the mean planet masses in each $r_p$ bin demonstrating the accuracy of our unique correction factors used to recover the mean mass-radius relation. Note the large uncertainties in the mean planet masses resulting from the large dispersion observed in the empirical mass-radius relation (Weiss & Marcy 2014).

Restricted range of $P$ and $r_p$. The second caveat with the Kepler-derived $f$ is that we are restricted in our range of sampled orbital periods and planetary radii as imposed by the size of the parameter space for which $f$ is robustly measured with Kepler. Specifically, $f(P, r_p)$ from Dressing & Charbonneau (2015) is restricted to the orbital period domain of $P \in [0.5, 200]$ days and the planetary radius domain of $r_p \in [0.5, 4]$ $R_⊕$, which includes small Mars-sized planets up to about the size of Neptune. However, given that the duration of the SLS-PS is longer than 200 days, it is important to also consider planets with orbital periods > 200 days. To inform the inclusion of such long-period planets, we use the RV M dwarf survey detection of the non-zero frequency of giant planets ($m_p \sin i \ge 100 M_⊕$) with orbital periods of > 200 days albeit with a low $f(\lesssim 5\%$; Bonfils et al. 2013a). We therefore supplement the Kepler occurrence rates with draws from the HARPS occurrence rates of giant planets (Bonfils et al. 2013a). Long-period planet sampling is carried out similarly to how the Kepler occurrence rates are sampled. Here we draw the minimum planetary masses rather than planetary radii, $f \to f(P, m_p \sin i)$. We also note that in Bonfils et al. (2013a), the $f$ for giant planets is computed over a coarse grid in minimum mass spanning an order of magnitude ($10^2 \le m_p \sin i/M_⊕ \le 10^3$) despite the most massive HARPS
detection having only $m_p \sin i \sim 112 M_\oplus$. Therefore, it is
possible that giant planets with minimum masses in excess of
$\sim 112 M_\oplus$ (0.3 $M_{Jup}$) do not exist in nature around M dwarfs. We
note that throughout this study, we will primarily focus on the
population of small planets from the *Kepler*-derived $f(P, r_p)$
because of their much higher frequency compared to giant
planets around M dwarfs.

Applicability of occurrence rates to the full input catalog.
Robust statistics regarding *f* were derived from the sample of
small *Kepler* stars, which almost exclusively contained early-
to-mid M dwarfs and late K-dwarfs. The assumption that the
resulting *f* extends to later spectral types is still largely
uncertain, but the early discovery of seven transiting Earth-
sized planets around the ultracool dwarf TRAPPIST-1 (Gillon
et al. 2017), from a small sample of observed stars, hints that
small planets around late M dwarfs might be as common and
potentially more common than around the early M dwarfs
observed with *Kepler*. Preliminary estimates of *f* around late M
dwarfs with K2 suggest a potential lack of super-Earth-sized
planets on close-in orbits but have been insufficient to probe
the population of Earth-sized planets (Demory et al. 2016).
Theoretically, planet formation scenarios have also predicted
the existence of many such small planets on close-in orbits
around late M dwarfs (e.g., Alibert et al. 2013; Alibert &
Benz 2017). In Section 10, we will investigate the effect of
increasing *f* on our planet detections in the SLS-PS.

Dynamical considerations of multiplanet systems. The final
caveat arises from $f(P, r_p)$ being derived in uncorrelated bins
whereas dynamical constraints will prevent certain types of
planetary systems from existing in nature, for example, close
pairs of massive planets. To ensure that sampled multiplanetary
systems in our simulations are dynamically stable, we impose
two priors on each system with multiplicity >1. The first
constraint is the analytic assessment of Lagrange stability
(Barnes & Greenberg 2006), which only depends on the masses
of the central star and planets, the planets' semimajor axes, and
eccentricities. For adjacent planet pairs that are Lagrange
stable, by definition their ordering remains fixed, both planets
remain bound to the central star, and the criterion limits
permissible changes to planets' semimajor axes. Lagrange
stability can be thought of as a more stringent extension of Hill
stability (Gladman 1993). However, we note that the analytic
treatment of Lagrange stability is only applicable to a three-
body system. In planetary systems from our simulations with
>2 planets, we apply the Lagrange stability criterion to every
adjacent planet pair. We then supplement the Lagrange stability
criterion with the heuristic criterion from Fabrycky et al.
(2012), which is applicable to systems with multiplicity >2.
The resulting criterion for stability<fnref n="13" /> is derived from the
stability analysis of a set of numerical integrations. In our MC
simulations, we only include multiplanetary systems, which
satisfy both aforementioned stability criteria.

As mentioned previously, our dynamical considerations
cause the injected planet population to not exactly match the
adopted *Kepler* occurrence rates. The modified *Kepler* planet
occurrence rates, or equivalently our injected planetary
population, are shown in Figure 7 as a function of the input
variables *P* and *r*<sub>*p*</sub>. We then report the same population after

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Top: histograms of the injected small planet population in the simulated SLS-PS as a function of orbital period and planetary radius from Dressing & Charbonneau (2015). Bottom: the same planet population as above converted to planetary mass using Equations (7) and (8). For clarity, each *r*<sub>*p*</sub> bin, and the approximately corresponding *m*<sub>*p*</sub> bin can be viewed individually: *r*<sub>*p,min*</sub> = 0.5 *R*<sub>⊕</sub>, *r*<sub>*p,min*</sub> = 1 *R*<sub>⊕</sub>, *r*<sub>*p,min*</sub> = 1.5 *R*<sub>⊕</sub>, *r*<sub>*p,min*</sub> = 2 *R*<sub>⊕</sub>, *r*<sub>*p,min*</sub> = 3 *R*<sub>⊕</sub>.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Probability density function of planet multiplicities in the simulated SLS-PS. The fraction of simulated systems with a given multiplicity—in percentage—is annotated for each integer multiplicity. We find an average planet multiplicity of ~2.4, which is highlighted by the vertical dashed line and is slightly less than the cumulative injected multiplicity of 2.5 planets with $P \le 200$ days from Dressing & Charbonneau (2015).</figcaption></figure>

converting the planetary radii into masses using Equations (7)
and (8).

The resulting distribution of planetary system multiplicities
is shown in Figure 8. Simulated planetary systems that obey

<footnote n="13">¹³ Long-term stability requires $a_{\text{out}} - a_{\text{in}} \gtrsim 3.5 R_{\text{Hill}}$, where $a_{\text{in}}$ and $a_{\text{out}}$ are the semimajor axes of the inner and outer planets, respectively, and $R_{\text{Hill}}$ is their mutual Hill radius.</footnote>
our dynamical stability criteria contain zero to seven planets, although only 0.04% of simulated planetary systems can survive with seven planets. Similarly, ~4% of the simulated planetary systems contain no planets at all. The most common planet multiplicity is 2, with ~32% of the simulated planetary systems containing two planets. The resulting average planet multiplicity is ~2.4, which is slightly less than the cumulative injected multiplicity of 2.5 ± 0.2. When recovering the planet yield of our simulated survey, we will have to correct for this small discrepancy between the cumulative planet multiplicity of our injected population and the true multiplicity of 2.5 for planets with $P \in [0.5, 200]$ days and $r_p \in [0.5, 4]$ $R_⊕$ (see Section 9).

# 6. Activity Mitigation

## 6.1. Overview of the Gaussian Process Formalism

RV activity signals present in M dwarfs (e.g., Gl 176, Forveille et al. 2009; Gl674, Bonfils et al. 2007; Proxima Centauri, Robertson et al. 2016) will deter our ability to detect new exoplanets in the SLS-PS. A number of correction techniques have been developed by various groups within the field to mitigate these effects and detect planets in RV (see reviews in Fischer et al. 2016; Dumusque et al. 2017). One particularly promising method is the use of Gaussian process regression to model correlated noise (i.e., activity) in RV time series.

Gaussian processes (GPs) belong to a class of non-parametric regression models from the field of machine learning (Rasmussen & Williams 2005). Being non-parametric, the functional form of the activity model is unspecified and instead is determined by the data itself. In this way, the GP activity model does not assume any physical model of the underlying processes responsible for the observed activity. Instead, the GP is used to model the covariance properties of an input time series according to a user-defined covariance function which itself is described by a small number of hyperparameters.<fnref n="14" /> Because of its non-parametric nature, GP regression is an attractive method for modeling the stochastic processes that give rise to observable stellar activity.

The GP prior distribution, which is described by the aforementioned covariance function, is a multivariate Gaussian distribution of functions specified by a mean function $\mu(t)$, evaluated at the epochs contained in the window function $t$, and a covariance matrix,

$$K_{ij} = \sigma_i^2 \delta_{ij} + k(t_i, t_j), \quad (9)$$

which is computed from the vector of measurement uncertainties $\sigma(t)$, the Kronecker delta function $\delta_{ij}$, and the user-defined covariance function $k(t_i, t_j)$ describing the covariance between two measurements taken at times $t_i$ and $t_j$ for $i, j \in [1, n_{\text{obs}}]$. Obtaining the GP model of an arbitrary input time series $y$ and uncertainty vector $\sigma$ is done by maximizing the Gaussian logarithmic likelihood function

$$\ln \mathcal{L} = -\frac{1}{2} (y - \mu)^T K^{-1} (y - \mu) + \ln \det K + n_{\text{obs}} \ln 2\pi), \quad (10)$$

to obtain the “best-fit” values of the GP hyperparameters describing the covariance of the observations $y$ through the covariance function $k(t_i, t_j)$. In the case of M dwarfs, activity predominantly arises from rotationally modulated ARs, which also evolve in time due to their varying lifetimes, spatial distribution, and contrast. Because of this temporal evolution, we adopt a quasi-periodic covariance function of the form

$$k(t_i, t_j) = a^2 \exp \left[ -\frac{|t_i - t_j|^2}{2\lambda^2} - \Gamma^2 \sin^2 \left( \frac{\pi|t_i - t_j|}{P_{GP}} \right) \right], \quad (11)$$

described by four hyperparameters: $a$ the amplitude of the covariance in units of $y$, $\lambda$ the exponential decay of correlations, $\Gamma$ the coherence scale of correlations, and $P_{GP}$ the periodic timescale. From identical reasoning to our own, quasi-periodic covariance functions are commonly adopted in a number of related astrophysical applications such as the recovery of stellar rotation periods (i.e., setting $P_{GP} = P_{\text{rot}}$; e.g., Angus et al. 2017) and the modeling of RV activity from spectroscopic activity diagnostics (e.g., Haywood et al. 2014; Rajpaul et al. 2015; Cloutier et al. 2017a), photometry (e.g., Cloutier et al. 2017a, 2017b), or the raw RVs themselves (e.g., Faria et al. 2016; Donati et al. 2017; Yu et al. 2017), and is applicable when the baseline of the observations spans at least a few rotation periods (Pont et al. 2013).

Following the optimization of the GP hyperparameters, one obtains a unique GP prior distribution. Conditioning the GP prior on the data set $y(t)$ results in the GP predictive distribution, the mean function and posterior variance of which can be evaluated at previously unseen epochs $t^*$ via

$$\mu(t^*) = K(t^*, t) \cdot K(t, t)^{-1} \cdot y(t) \quad (12)$$

and

$$C(t^*) = K(t^*, t) - K(t^*, t) \cdot K(t, t)^{-1} \cdot K(t^*, t)^T. \quad (13)$$

Here, the covariance matrix $K(t^*, t)$ is evaluated at unseen epochs such that $\sigma(t^*)$ must be set to zero in Equation (9).

## 6.2. Modeling RV Activity with GP Regression

As in nature, the RV time series from our MC simulations contain contributions from both planetary companions as well as activity from ARs. These sources of RV activity signals also have manifestations in other spectroscopic time series, albeit not strictly at the same rotation period (Hébrard et al. 2016). One example is the FWHM of the CCF which we simulate along with the RV signals (see Section 3.2). As is often the case when searching for small planets in RV, one cannot distinguish a priori an RV signal from a planet or from activity. However, the covariance properties in the FWHM time series will be related to the covariance properties in RV residuals after the removal of planetary sources. We therefore use the FWHM time series to train the hyperparameters of our quasi-periodic GP prior distribution. In principle, one could have chosen an alternative activity indicator such as the BIS or the contrast of the CCF; however, recent HARPS observations have demonstrated that the strongest correlations between the RVs and an activity indicator often exist between the RV and FWHM time series (e.g., Astudillo-Defru et al. 2017; X. Bonfils et al. 2018, in preparation). However, this is not universal. Similarly, the contrast or depth of the CCF is not used in place of the FWHM

<footnote n="14">¹⁴ “Small” compared to the size of the input time series, which would be required to fully describe the covariance properties of the input data set.</footnote>
<table><caption>**Table 1**
Gaussian Process Hyperparameter Priors Used in Training</caption><thead><tr><th>Hyperparameter</th><th>Prior</th></tr></thead><tbody><tr><td>Covariance amplitude, <i>a</i> [m s<sup>-1</sup>]</td><td>J(10<sup>-2</sup>, 10<sup>2</sup>) · max|FWHM - &lt;FWHM&gt;|)</td></tr><tr><td>Exponential decay timescale, λ [days]</td><td>J(1, 10 · (max<i>t</i> - min<i>t</i>))</td></tr><tr><td>Coherence scale, Γ</td><td>J(10<sup>-2</sup>, 10<sup>2</sup>)</td></tr><tr><td>Periodic timescale, <i>P</i><sub>GP</sub> [days]</td><td>J(0.1, 300)<sup>a</sup><br>U(0.9, 1.1) · <i>P</i><sub>tot</sub><sup>b</sup></td></tr></tbody></table>

**Notes.**

<sup>a</sup> If $P_{\text{tot}}$ is not detected in the FWHM time series.

<sup>b</sup> If $P_{\text{tot}}$ is detected in the FWHM time series.

because although the two time series are highly correlated, the CCF contrast is known to be more strongly affected by instrumental noise in practice. Using a spectroscopic training set like the FWHM also has the benefit of being obtained contemporaneously with the RV measurements therefore probing the star's activity at the same epochs in which we are searching for planets. In this way, our GP activity model is also sensitive to variations in the RV activity occurring on long timescales such as that from magnetic activity cycles.

To reduce the computational wall time of our full MC simulation, we do not compute a GP activity model in each MC realization of our simulated SLS-PS. The GP activity model is only computed when the rms of the injected RV activity exceeds the median RV measurement uncertainty. In such cases, we run a Markov Chain Monte Carlo (MCMC) on the FWHM time series to obtain the marginalized posterior probability density functions (PDF) of the four GP hyperparameters: $a$, $\lambda$, $\Gamma$, and $P_{\text{GP}}$. To sample the posterior PDFs, we use the emcee affine invariant MCMC ensemble sampler (Foreman-Mackey et al. 2013) coupled with the fast GP package george (Ambikasaran et al. 2015) to evaluate the likelihood function in Equation (10). Specifically, we initialize 100 walkers with an effective chain length of ~10 autocorrelation times to ensure the convergence of the chains. Samples of the posterior PDFs are only saved following a burn-in phase of ~10 autocorrelation times. From a set of supervised preliminary tests of the MCMC procedure, walkers are initialized in the parameter space within Gaussian balls whose variance is chosen such that the mean acceptance fraction among the walkers is 20%–60%.

We select broad non-informative priors for all GP hyperparameters with the exception of $P_{\text{GP}}$, which is constrained to the narrow uniform range of $P_{\text{GP}}/P_{\text{rot}} \in U(0.9, 1.1)$ either when the photometric rotation period is known for a particular star or if $P_{\text{rot}}$ is detected in an LS periodogram of the FWHM with a false-alarm probability (FAP) $\le 1\%$. The prescriptions for the GP hyperparameter priors are reported in Table 1. Throughout this paper, we calculate FAPs via bootstrapping with replacement using $10^4$ iterations, and each LS periodogram is normalized by its standard deviation. From the marginalized posterior PDFs of the GP hyperparameters, we adopt the maximum a posteriori (MAP) values of each hyperparameter to construct a unique covariance matrix $K$ and thus our mean GP model and its $1\sigma$ confidence interval from Equations (12) and (13). Examples of resulting GP models of the FWHM time series for a rapidly and a slowly rotating star are shown in the upper panels in Figure 9.

Following the training phase on the FWHM time series, we proceed with modeling the RVs simultaneously with a trained GP activity model plus Keplerian planetary signals. The marginalized posterior PDFs of the GP hyperparameters $\lambda$, $\Gamma$, and $P_{\text{GP}}$ from training are used as informative priors in the joint RV analysis, which treats the GP amplitude $a$ as a free parameter. For each assumed mean function $\mu$, containing between zero and three Keplerian solutions, we compute the MAP GP activity model from the hyperparameter values sampled using MCMC. Assuming a zero-mean function, the resulting mean GP activity models for the two stars shown in Figure 9 are shown in the middle panels of the figure. The residuals following the removal of the mean GP activity model is also shown. In each case, the stellar rotation period is detected in the LS periodogram of the FWHM time series and therefore is used to constrain $P_{\text{GP}}$ during training. For the rapid rotator with $P_{\text{rot}} = 2.8$ days, the rms of the injected activity is 4.43 m s⁻¹ and is reduced to 1.59 m s⁻¹ after removing the mean GP activity model ($\chi^2_{\text{red}} = 1.4$ for four GP hyperparameters). This resulting RV rms is more comparable to the median RV measurement uncertainty $\sigma_{\text{RV}} = 1.35$ m s⁻¹. In the slow rotator case ($P_{\text{rot}} = 83.1$ days), the rms of the injected activity is reduced from 3.60 to 1.53 m s⁻¹ compared to $\sigma_{\text{RV}} = 1.63$ m s⁻¹ ($\chi^2_{\text{red}} = 0.9$). In both test cases considered in Figure 9, a planet is detected in the RV residuals and fit simultaneously with the activity, assuming a new mean model containing a single Keplerian solution. Details of our planet detection algorithm are discussed in Section 7.

### 6.3. GP Activity Model Performance

As an overview of the performance of our GP activity modeling, we can compare the rms of the known injected RV activity with the rms of the residuals following the removal of our mean GP activity model. The residual rms should never exceed the rms of the injected activity, otherwise our GP formalism would be adding additional noise into the RV time series rather than modeling and reducing it, which is its intended purpose. Similarly, optimal GP fits will result in a residual rms that is close to the median RV measurement uncertainty of the time series.

Figure 10 compares the rms of the injected RV activity to the rms of the residual RV activity after removing the mean GP activity model computed assuming a zero planet model. Recall that we only compute a GP activity model when the injected rms is $> \sigma_{\text{RV}}$ such that the injected activity rms in units of $\sigma_{\text{RV}}$ is always greater than or equal to unity in Figure 10. We note that the residual activity rms never exceeds the injected activity rms as expected, i.e., the residual rms always lies beneath the line $y = x$ in Figure 10. Second, there appears to be a positive correlation between $P_{\text{rot}}$ and the relative reduction of the activity rms, which is analogous to GP performance. That is within our sample of RV time series, which are modeled with GP activity model; the activity rms is maximally reduced in systems with the longest rotation periods. Conversely, the GP performance when applied to rapid rotators ($P_{\text{rot}} \gtrsim 2$ days) is often marginal in comparison. The exact cause of these effects above may be related to the poor time sampling of our observations compared to the short stellar rotation period but is ultimately beyond the scope of this paper and is reserved for a future study (R. Cloutier et al. 2018, in preparation).

We note that in many slow rotator cases, the dimensionless residual rms is often less than unity. This suggests that in such
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Examples of our GP formalism used to mitigate RV activity and detect underlying planetary signals for two simulated systems with either a rapidly ($P_{\text{rot}}$ = 2.8 days; left column) or slowly ($P_{\text{rot}}$ = 83.1 days; right column) rotating host star. Top panels: the FWHM time series used for training along with the mean GP regression model and its $1\sigma$ confidence interval shown in red. Middle panels: the raw RV time series along with the mean GP activity model (in the absence of planets) and its $1\sigma$ confidence interval. Bottom panels: the RV residuals after removal of the mean GP activity model shown in the middle panels. Each system contains three planets that contribute to the O-C residuals albeit with semi-amplitudes that are all $\lesssim\sigma_{\text{RV}}$. The two upper rows report the MAP values of GP hyperparameters along with their 16th and 84th percentiles.</figcaption></figure>

cases, the non-parametric GP is actually modeling the noise and not just the stellar activity signal. Unfortunately, when the GP activity model overfits the RV noise, planetary signals can be absorbed into the mean GP model and thus avoid detection. Furthermore, this apparent overfitting is a common feature in many of our simulations and does not appear to depend on whether or not the stellar rotation period is detected a priori in the FWHM time series. However, it is true that on average, the GP modeling outperforms cases in which $P_{\text{rot}}$ is known from the FWHM compared to cases in which $P_{\text{rot}}$ remains undetected, thus providing very weak constraints on $P_{\text{GP}}$. This highlights the importance having a priori knowledge of $P_{\text{rot}}$ from any of the ancillary spectroscopic time series, polarimetric time series (Hébrard et al. 2016), or from previously obtained long-baseline photometry (e.g., Newton et al. 2016).

## 7. Automated Planet Detection

Due to the large number of planetary systems in our simulated SLS-PS, we must detect planets in an automated way. The steps in our automated planet detection algorithm represent computationally tractable calculations given the large number of planetary systems for which each step must be performed. We note however that other—potentially more robust—planet detection algorithms may be adopted in the real SLS-PS, which are likely to include more human intervention than the automated techniques described in the following subsections.

### 7.1. Establishing Putative Planetary Detections

We proceed by searching for planetary periodicities in an iterative manner using the LS periodograms of the RVs following the removal of various periodic signals. An example of this iterative process is visualized in Figure 11 for a three-planet system ($P_b, P_c, P_d$ = 1.85, 7.17, 14.66 days; $K_b, K_c, K_d$ = 4.7, 2.3, 0.07 m s⁻¹, respectively) around a moderately active star with a measured photometric rotation period $P_{\text{rot}} \sim$ 8.8 days and an RV rms in the absence of planets of $\sim$4.2 m s⁻¹.

For MC realizations in which the rms of the injected RV activity exceeds the median RV measurement uncertainty—as is the case for the system shown in Figure 11—we compute two versions of our initial periodogram: one of the raw RVs *only* and a second of the RV residuals after fitting the data with the trained GP activity model and zero-mean function (i.e., no-planet model). This GP model predicts the RV activity in the absence of planetary signals. For the remaining MC realizations containing quiet stars, we only compute the LS periodogram of the raw RVs, thus ignoring any modeling of correlated RV residuals. For cases in which we use a trained GP to model RV activity, it is beneficial for the periodic term of the assumed quasi-periodic covariance kernel to be constrained by our training set. For the example shown in Figure 11, the top panel shows the LS periodogram of the FWHM in which $P_{\text{rot}}$ is detected and is subsequently used in our RV modeling (see Section 6.2).

The second panel in Figure 11 depicts the LS periodogram of the raw RVs only whereas the third panel depicts the raw RVs corrected by the trained GP activity model with $P_{\text{GP}} = P_{\text{rot}} =$
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Comparison of the rms of the injected RV activity to the rms of the residual RV activity following the removal of the mean GP activity model, which is computed in the absence of planets. Each rms value is normalized by the RV measurement uncertainty $\sigma_{RV}$ of its time series. The dashed line $y = x$ is indicative of cases wherein the GP activity model does not reduce the rms due to activity. Time series in which $P_{rot}$ is detected in the FWHM time series are depicted as circles, whereas time series in which $P_{rot}$ is not detected in the FWHM time series are depicted as triangles.</figcaption></figure>

8.8 days. Comparing these two periodograms, it is apparent that some power at $P_{rot}$ is diminished when removing the mean GP activity model along with power at long periods due to the large exponential timescale found during training. More importantly, the strongest periodic signal in each periodogram is at the orbital period of the innermost planet at ~1.85 days, suggesting the presence of a planet. In these periodograms, and in all subsequent iterations, we claim a putative planet detection for periodicities with (i) FAP $\le$ 1%, (ii) is within 2% of an injected planet's orbital period, (iii) is not associated with significant periodic signals seen in the FWHM time series, (iv) is not an alias of the time sampling (see example in Figure 4), and (v) is >2% away from the stellar rotation period and any of its first four harmonics. The signals at ~1.85 days in the second and third panels of Figure 11 obey these criteria and therefore constitute a putative planet detection. We proceed by referring to this putatively detected planet as “b.”

We note that the second condition for a putative planet detection above cannot be utilized in a real survey because any planetary periodicities are not known a priori. Instead, we invoke this condition to accelerate our automated planet detection algorithm relative to the steps that must be taken in a real survey to secure planet detections (e.g., Bayesian model comparison; Ford & Gregory 2007). Without a priori knowledge of planet orbital periods, LS periodicities at high significance that obey all of the remaining aforementioned putative planet criteria may represent false positive signals if not carefully modeled. In our simulations, we find that on average, our time series generate ~0.5 such false positives in their LS periodograms following the removal of planet models,

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Top to bottom: Lomb-Scargle periodograms of the FWHM of the cross-correlation function, raw radial velocities (RVs), RVs corrected for activity (RVs-GP), RVs corrected for activity and a planet (RVs-GP-b), and the RVs corrected for activity and two planets (RVs-GP-b-c). Periodicities equal to the stellar rotation period, its first harmonic, and the orbital periods of the three planets in the system are highlighted with vertical lines. The power corresponding to a 1% false-alarm probability in each periodogram is highlighted by a horizontal dashed line in each panel.</figcaption></figure>

although a fraction of these false positive periodicities are likely to be aliases of each other, making the above estimate an upper limit. Ideally, the determination of these signals as false positives or as true planetary signals would be solved via a formal model comparison in a subsequent analysis. However, such calculations are not guaranteed to converge to the correct solution (see Dumusque et al. 2017).

The second periodogram iteration requires that we fit the putative planetary signal and search for additional signals in the LS periodogram of the residuals. To model planet “b,” we adopt the maximum-likelihood Keplerian model parameters at the detected periodicity and recompute the GP activity model, which is modified due to the new mean function (i.e., a one-planet model rather than the previously assumed no-planet model). Note that the GP activity model is recomputed after each iteration, due to the changing mean function. The fourth panel in Figure 11 shows the resulting periodogram after correcting the raw RVs with the new GP activity model and a
Keplerian solution for planet “b.” The strongest residual periodic signal is at the orbital period of the middle planet at ~7.17 days, which we claim as a second putative planet detection because the periodicity obeys the aforementioned criteria. This second putative planet is referred to as “c.”

The bottom panel of Figure 11 depicts the LS periodogram of the RVs after being corrected for activity and the superposition of the two maximum-likelihood planet models. Now, the strongest residual periodicity has an FAP > 1% implying that our automated planet detection algorithm has ceased to detect planetary signals. Therefore, in the example shown in Figure 11, the third injected planet at ~14.66 days remains undetected—an unsurprising result given the small RV semi-amplitude of the planet ($K_d = 7$ cm s⁻¹) compared to the time series’ median RV measurement uncertainty $\sigma_{RV} = 1.8$ m s⁻¹. In each MC realization, we perform this iterative procedure until no putative planets are detected and up to a maximum of three planets, despite many planetary systems having >3 injected planets (see Figure 8). In this way, we are at least sensitive to the expected number of planets per M dwarf (2.5 ± 0.2; Dressing & Charbonneau 2015) and limit the computational expense of detecting planets dominated by repeatedly computing GP activity models and LS periodogram FAPs.

## 7.2. Model Selection

In Section 7.1, we established putative planet detections based on low FAP LS periodogram periodicities. However, the robust detection of a planet with a particular set of Keplerian model parameters must be favored over competing models that lack such a planet. The proper diagnostic for model selection is the ratio of Bayesian model evidences, which are notoriously difficult and time consuming to calculate (Ford & Gregory 2007). As an alternative model selection technique, we turn to time series cross-validation (CV). This technique is a specialized version of the general K-fold CV and is suitable for data featuring strong correlations in time as is the case with RV time series (Arlot & Celisse 2010).

For MC realizations featuring at least one putative planetary detection, we perform time series CV on models that contain an increasing number of planets, including the null hypothesis, i.e., no planets. The latter model has zero Keplerian parameters, whereas a model containing $N_p \ge 1$ planets contains $3N_p$ model parameters where the three parameters per planet are its orbital period, time of inferior conjunction, and RV semi-amplitude. For the purpose of model selection, we will assume circular orbits for all planets to limit the size of the parameter space.

The CV algorithm proceeds by first splitting the RV time series $y_1, \dots, y_{n_{obs}}$ into training and testing sets. For some $t > 1$, each competing model is fit to the training set $y_1, \dots, y_t$ using a Levenberg-Marquardt optimization routine. The optimized model is then evaluated at the next epoch $t+1$ ($\mu_{t+1}$), and the ln likelihood of the testing set $y_{t+1}$ given the optimized model is computed using Equation (10). When computing the ln likelihood, we adopt a white covariance matrix for systems wherein GP analysis is not used but otherwise assumes the MAP GP hyperparameters from the iterative procedure in Section 7.1. These steps are repeated for $t = N_{min}, \dots, N-1$, where the minimum size of the training set $N_{min}$ is set to 20. The favored model is determined by which of the competing models has the largest median ln likelihood per measurement

among the $N - N_{min}$ CV iterations. In cases wherein two models are consistent within their median absolute deviations, the model containing fewer planets is accepted as an imposition of Occam’s razor.

## 7.3. Vetting of Planet Detections

A consequence of our automated planet detection methodology is various non-deterministic effects, which can result in planet detections that are highly unlikely to be favored by model comparison in the real SLS-PS yet are marginally detected in our simulations. Such planets are commonly those whose RV semi-amplitude is close to the rms of the RV time series. These planets would likely be rejected by any human vetting, which we do not conduct in our simulations. We therefore undergo a vetting procedure in an attempt to restrict the detected planet population to be maximally realistic. Our adopted vetting procedure is based on the methods of Cumming et al. (2008) from the Keck Planet Search. For vetting, we define the condition that a bona fide planet detection must satisfy $K/\sigma_K \ge 3$, i.e., that a true planet detection is one in which the planet’s semi-amplitude K is detected with a minimum expected detection significance of $3\sigma$.

To estimate the expected uncertainty in the RV semi-amplitude $\sigma_K$, we compute the Fisher information matrix B, which quantifies the information content in an observable time series $y(t)$ regarding unknown model parameters $\theta$. The model parameter covariance matrix is related to the Fisher information via $C = B^{-1}$. We can therefore use the Fisher information matrix to analytically predict the measurement uncertainty of the RV model parameters of interest given an input time series with $n_{obs}$ measurements $y = (y_1, \dots, y_{n_{obs}})$ obtained at the epochs $t = (t_1, \dots, t_N)$ and with measurement uncertainties $\sigma = (\sigma_1, \dots, \sigma_N)$.

The Fisher information matrix is a Hessian matrix of the ln likelihood of a single Keplerian model with respect to its model parameters $\theta = \{\hat{P}, T_0, K\}$:

$$B_{i,j} = - \frac{\partial^2 \ln \mathcal{L}}{\partial \theta_i \partial \theta_j} \qquad (14)$$

The Fisher information matrix is symmetric and is $3 \times 3$ in our case as we assume circular orbits.

In order to simplify the calculation of B, we consider planets individually and account for the residual RV signal from additional planets through an “effective” RV uncertainty $\sigma_{eff}$ in place of the RV measurement uncertainties $\sigma_i$ in Equation (9). The effective RV uncertainty is the rms of the RVs after removal of the Keplerian signal from the planet being considered. It therefore contains contributions from any additional planets, stellar activity, and systematic errors. Because we do not fit for each planet’s orbital eccentricity, the Keplerian model simplifies to $\mu(t_k) = -K \sin \phi_k$, where $\phi_k = \frac{2\pi}{p}(t_k - T_0)$. Using this mean model in the ln likelihood (Equation (10)) along with the $K_{ij}$ approximated with a white covariance matrix $K_{ij} = \sigma_{eff} \delta_{ij}$, we can compute each element of B analytically (see the Appendix). The Fisher information matrix is then inverted to obtain the covariance matrix of the model parameters C. The measurement uncertainty of the semi-amplitude is then $\sigma_K = \sqrt{C_{K,K}}$.
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Number of planets detected as a function of $m_p \sin i$ after each of the three steps in our automated planet detection algorithm. The number of detected planets in each $m_p \sin i$ bin is annotated above the bin to help quantify the decrease in planet detections following each step. The distribution of detected planet minimum masses following each step can be viewed independently for clarity (1, 2, 3).</figcaption></figure>

## 7.4. Summary of the Automated Planet Detection Algorithm

To recapitulate our process of claiming planet detections in our simulated SLS-PS, we recall the three steps discussed throughout this section. First, we search for putative planetary signals in the LS periodogram of either the raw RVs or the RVs corrected for activity using a zero-mean GP activity model. We are careful to ensure that significant periodicities have a planetary origin and are not associated with stellar activity signals or our time sampling. Second, we compute model in likelihoods using time series CV, and compare models with and without the putative planet. We only retain planets favored by this model selection technique. The third and final step consists of vetting our planet detections by insisting that they must have an RV semi-amplitude detection significance greater than $3\sigma$, where the detection significance is estimated from the planet's known semi-amplitude *K* and an analytical estimate of the measurement uncertainty on *K* from the Fisher information. Planet detections that pass our vetting procedure are treated as bona fide detections.

The distributions of detected planet minimum masses after each step in our automated planet detection algorithm are shown in Figure 12. Here we only include planet detections with $0.1 \le m_p \sin i/M_\oplus \le 20$. In this SLS-PS, there are a small number of giant planets with $m_p \sin i > 20 M_\oplus$ detected, thus resulting in an underestimated total planet yield annotated in Figure 12 (see Section 9 for a full description of the detected planet population).

It is clear from Figure 12 that each step in our automated planet detection algorithm reduces the number of detected planets somewhat. In each $m_p \sin i$ bin other than the most massive bin ($9.4 \le m_p \sin i/M_\oplus \le 20$), our CV model selection technique rejects between ~1 and 3 planets, or ~13% of all putative planet detections in the range of minimum planet masses considered. Similarly, our vetting procedure reduces the number of detected planets in each $m_p \sin i$ bin by ~1–2 planets for intermediate minimum masses, $0.7 \le m_p \sin i/M_\oplus \le 6$. Our vetting procedure therefore rejects a significant number of

detected planets neither at the lowest masses ($m_p \sin i \lesssim 0.7 M_\oplus$) nor at the highest ($m_p \sin i \gtrsim 6 M_\oplus$), the former being the result of the small number of putative low-mass planets detected and the latter being due to the large RV semi-amplitude of the most massive planets thus resulting in a typically large detection significance. Vetting rejects ~7% of planet detections favored by CV. Therefore, ~80% of putative planet detections from the periodogram analysis materialize into bona fide planet detections.

## 8. SLS-PS Sensitivity

The detection sensitivity is defined as the recovery fraction of injected planets in our simulated SLS-PS. Because we have a priori knowledge of the injected planet population, we can compute the detection sensitivity for each star in our sample by simply dividing the number of detected planets by the number of injected planets over any desired range of planet properties. We perform this calculation over the discretized parameter space in *P*, $m_p \sin i$, and in *S*, $m_p \sin i$. Here, $S = L_s/4\pi a^2$ is the insolation received by the planet, where $L_s$ is the stellar luminosity<fnref n="15" /> and *a* is the planet's semimajor axis. We focus on the following ranges of parameter values, which encompass the vast majority of the injected planet population: $P \in [0.5, 200]$ days, $S \in [0.01, 100]$ $S_\oplus$, and $m_p \sin i \in [0.4, 15]$ $M_\oplus$.

We note that in this study the recovery fraction is uniquely determined by the performance of our automated planet detection algorithm (see Section 7). Conversely, the actual SLS-PS will have a much higher degree of human intervention on the data analysis effort. This is afforded by the relatively small size of the RV data sets compared to large surveys (e.g., *Kepler* and *TESS*), which benefit greatly from automated detection algorithms. Therefore, the detection sensitivity in the actual SLS-PS may not correspond exactly to what is presented here although the automated algorithm used in this study is designed to closely mimic the analysis that will be conducted on the actual SLS-PS data.

The detection sensitivity to planets varies from star to star due to their changing stellar properties, which can affect our ability to detect planets in RV (e.g., apparent magnitude, stellar mass, stellar rotation, etc.). Computing the detection sensitivity for each star individually is necessary for calculating planet occurrence rates (see Section 11). To improve the detection statistics across the full range of planetary parameters considered, we augment the MC realizations for each star with an additional set of planetary systems with logarithmic *P*, *S*, and $m_p \sin i$ sampled uniformly rather than from the planet occurrence rates. The individual detection sensitivity maps for each star can then be combined to obtain the average sensitivity maps for the full SLS-PS as a function of *P*, *S*, and $m_p \sin i$ as shown in Figure 13. In this way, we marginalize over the stellar properties of our sample stars including the aforementioned parameters, which are known to influence the detection sensitivity for each individual star. Hence, our sensitivity results might be scaled to various stellar samples provided that their global properties are consistent with those of our current sample. In Figure 13, the uncertainties in the detection sensitivity within each grid cell comes from counting or Poisson statistics, and therefore benefit from a large number of simulated planetary systems.

<footnote n="15">¹⁵ The stellar luminosity is computed from the evolutionary models of Baraffe et al. (1998) based on the stellar mass on the main sequence at 2 Gyr.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Binned maps of the detection sensitivity of the full SPIRou input catalog derived using our automated planet detection algorithm as a function of minimum planet mass and orbital period (left) or insolation (right). The dashed vertical lines in the insolation panel indicate the approximate “water-loss” and “maximum-greenhouse” insolation limits of the HZ from Kopparapu et al. (2013). The dashed-dotted vertical lines indicate the more conservative “recent-Venus” and “early-Mars” HZ limits (Kopparapu et al. 2013) Shaded bins highlight regions of the parameter space wherein our detection sensitivity is < 20%. The solid lines highlight the curve with $K = 1 \text{ m s}^{-1}$ for a star with mass equal to our sample’s median stellar mass of 0.25 $M_\odot$.</figcaption></figure>

Given the binning in Figure 13, nowhere do we achieve a 100% detection sensitivity. This is true even for the most massive close-in planets whose RV semi-amplitudes are typically much greater than the characteristic RV measurement uncertainty so long as the system is not orientated close to face-on (i.e., $\sin i \sim 0^{\circ}$). The geometric effect as well as the potential aliasing of periodic signals arising from stellar rotation or the window function also prevent the detection of certain types of planets. For example, in Figure 2 we see that in our stellar sample there is a dearth of $P_{\text{rot}} \sim 3–10$ days. Thus, we run into minimal aliasing from $P_{\text{rot}}$ at those orbital periods and achieve an increased detection sensitivity relative to planets with equivalent $m_p \sin i$ but with smaller $P$ (left panel Figure 13).

Similarly, we can see a steep decrease in detection sensitivity at periods of ~30 days. Recall that in our window functions we are restricted by the telescope’s observing schedule to only observe outside of dark time, whose cycle follows the lunar cycle with a cadence of ~30 days. Consequently, the origin of the decrease in detection sensitivity at ~30 days can be traced back to aliasing from our time sampling at that period, an unfortunate circumstance as the range of orbital periods corresponding to the HZ around M2-4 dwarfs spans 30 days.

Due to the aforementioned aliasing effects from stellar rotation and our window functions, the SPIRou detection sensitivity is maximized for the most massive planets with $3 \lesssim P \lesssim 10$ days. Another pertinent effect at orbital periods close to one day is the potential for signal aliasing due to the Earth’s rotation. This—in part—is responsible for the low detection sensitivity at orbital periods close to a day in Figure 13, another unfortunate circumstance given the general interest in quantifying the occurrence rate of close-in planets (Mulders et al. 2015).

The average detection sensitivity across the full range of $P$ and $m_p \sin i$ in Figure 13 is 44.8% ± 0.5%. This average detection sensitivity is sufficiently high such that the resulting planet detections from the SLS-PS will be able to place strong

constraints on the cumulative occurrence rate of planets around SPIRou stars (see Section 11).

When considering the detection sensitivity as a function of insolation, we see the same qualitative structure seen as a function of orbital period (right panel of Figure 13). For example, the increased sensitivity between ~3 and 10 days has a broad manifestation at $S \in [1, 10] S_\oplus$. Similarly, in both cases, it is unsurprising to see that the detection sensitivity increases toward more massive planets. In both cases, the SPIRou detection sensitivity reaches its lowest values for the least massive planets on wide orbits. Notably, we achieve a detection sensitivity of $\lesssim 20\%$ for all planets with $m_p \sin i \lesssim 1 M_\oplus$, thus making it difficult to detect Earth-mass planets and smaller in the SLS-PS. The average detection sensitivity across the full range of $S$ and $m_p \sin i$ considered here is 47.6% ± 0.5%.

## 8.1. Detection Sensitivity to HZ Planets

It is also important to consider our detection sensitivity to HZ planets as these targets are often flagged for various observational follow-up campaigns. In this study, we adopt the “water loss” and “maximum greenhouse” definitions as our fiducial HZ limits from Kopparapu et al. (2013), which are derived from a 1D radiative-convective climate model in the absence of clouds. Following Kasting et al. (1993), the inner edge of the HZ is defined by the “water-loss” limit, which arises from the photolysis of water in the upper atmosphere and subsequent hydrogen escape. The outer edge of the HZ is determined by the “maximum-greenhouse” limit, wherein an increase in atmospheric $CO_2$ levels will result in a net cooling as the increased albedo from Rayleigh scattering begins to dominate over the increasing greenhouse effect. The insolation levels approximately corresponding to our adopted HZ definition are $S \in [0.22, 0.90] S_\oplus$ for our stellar sample. For comparison, we also consider the more conservative ($S \in [0.20, 1.48] S_\oplus$; Kopparapu et al. 2013) “recent-Venus” and “early-Mars” HZ limits, which assume that both Venus and Mars were habitable early on in the lifetime of the solar system.
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Coarsely binned maps of the predicted planet yield from the SLS-PS as a function of minimum planet mass and orbital period (left) or insolation (right). The dashed vertical lines in the insolation panel indicate the approximate “water-loss” and “maximum-greenhouse” insolation limits of the HZ from Kopparapu et al. (2013). The dashed-dotted vertical lines indicate the more conservative “recent-Venus” and “early-Mars” HZ limits (Kopparapu et al. 2013). The annotated numbers in each bin report the predicted number of detected planets and the uncertainties on the prediction, which are dominated by uncertainties in the input planet occurrence rates.</figcaption></figure>

From the fiducial definition of the HZ, we find an average detection sensitivity to HZ planets of 47.3% ± 1.0%, which is consistent with the average detection sensitivity over the full $S, m_p \sin i$ grid in the right panel of Figure 13 (46.2% ± 0.5%). Adopting the more generous HZ limits, we find a comparable average detection sensitivity of 48.7% ± 0.9%. However, the average detection sensitivity to Earth-like planets ($m_p \in [1, 5] M_⊕$) in the HZ is significantly reduced to 36.7% ± 1.2%. Here, we have defined Earth-like planets in terms of their absolute mass, where the mass upper limit is approximately equal to the planet mass obtained when evaluating the mean mass-radius relation (Equation (7)) at the proposed maximum radius of a rocky planet, ~1.5–1.8 $R_⊕$ (Weiss & Marcy 2014; Rogers 2015; Fulton et al. 2017).

## 9. SLS-PS Predicted Yield

For each star in the SPIRou input catalog, we compute the number of planets detected as a function of $P, S$, and $m_p \sin i$ given the input occurrence rates and calculations of each star's detection sensitivity. Here we must correct the reduced injected cumulative planet occurrence rate of 2.4 planets per star over our grid of $P \in [0.5, 200]$ days and $m_p \sin i \in [0.4, 15] M_⊕$ to be equal to the intended 2.5 planets per star over the Dressing & Charbonneau (2015) grid ($P \in [0.5, 200]$ days, $r_p \in [0.5, 4] M_⊕$). Then after dividing out the number of simulated planetary systems per star, we compute the total planet detection yield predicted by our simulated survey. The predictions are shown in Figure 14 over a more coarsely binned map than in Figure 13, due to the small number of planets detected in each bin.

The cumulative planet yield over the orbital period domain and minimum mass range considered in Figure 14 is 85.7$^{+29.3}_{-12.5}$ out of ~180 injected planets. Of these, ~53.7 planets (62.6%) are the only planets detected in the system while ~26.8 planets (31.3%) are detected in a two-planet system. The remaining ~5.2 planets (6.1%) are found in systems with three detected

planets, i.e., we expect to detect one to two three-planet systems in the SLS-PS.

The number of simulated planetary systems per SPIRou star is large in our simulations. The result is that the uncertainties in the predicted yield are dominated by uncertainties in the input planet occurrence rates from *Kepler*. Due to *Kepler's* low detection sensitivity to small planets on wide orbits, the planet occurrence rate and hence the predicted SPIRou yield is poorly constrained at large orbital periods/low insolation levels. Approximately half of our planet detections are super-Earths with minimum masses $m_p \sin i \in [3, 7] M_⊕$ owing to their assumed frequency, which is intrinsically high, and the moderately high detection sensitivity achieved across the range of orbital periods considered (~30%–85%). Considering the “water-loss” and “maximum-greenhouse” definitions of the HZ limits from Kopparapu et al. (2013), we detect $22.0^{+18.4}_{-7.9}$ out of ~47 injected HZ planets in the SLS-PS over the range of minimum masses considered in Figure 14. These include $9.0^{+8.5}_{-3.6}$ Earth-like HZ planets out of ~25 injected planets with $m_p \in [1, 5] M_⊕$. When adopting the more generous “recent-Venus” and “early-Mars” HZ limits, these numbers get inflated to $31.5^{+26.8}_{-11.4}$ and $13.4^{+12.7}_{-5.3}$, respectively.

The population of SPIRou planets can be visualized slightly differently in the insolation/minimum planet mass plane as shown in Figure 15. Here we present a random subset of all simulated planetary systems. The size of this subset is chosen such that the integer number of detected planets in the subset is consistent with the predicted planet yield of 89.9 planets. The resulting total number of planets shown in Figure 15 is 250 with 90 planets detected. Because the subset of planets shown in Figure 15 is random, it does not preserve the number of detected planets in various subsets of the planet parameter space. For example, in the full simulated SLS-PS, we detected nine Earth-like planets whereas the random subset shown in Figure 15 only depicts four Earth-like planet detections.

The stellar parameters, planetary parameters, and simulated time series for each MC realization are made available to the
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Random subset of the simulated SPIRou planets representative of the underlying planet population investigated in the SLS-PS in the insolation/minimum planet mass plane. SPIRou planet detections are marked by solid circles whereas open circles represent injected planets that remain undetected by SPIRou. We detect 90 planets around 100 stars in the subset of simulated planets shown here. The size of each planet's marker is proportional to its radius. The inner shaded region highlights the approximate "water-loss" and "maximum-greenhouse" limits of the HZ, whereas the outer shaded region highlights the "recent-Venus" and "early-Mars" HZ limits (Kopparapu et al. 2013).</figcaption></figure>

community on GitHub<sup>16</sup> in the form of Python pickles. We also provide the combined results of all realizations for interested users to analyze the full SPIRou input catalog and each star's suite of simulated planetary systems. These data may be used, for example, to reconstruct most of the figures shown in this paper. Instructions and examples of how to read and access these data are also provided along with all the required Python scripts.

## 9.1. Giant Planet Detections

Although the vast majority of our predicted SPIRou planet population are derived from the *Kepler* occurrence rates of small planets ($r_p \le 4 R_⊕$ or $m_p \sin i \lesssim 15 M_⊕$), the non-zero frequency of giant planets from RV surveys results in some giant planet detections. From our simulations, we find that 4.1 SPIRou detections are giant planets with $m_p \sin i \ge 20 M_⊕$. Thus, the total SPIRou planet yield including all planets becomes $89.9_{-13.0}^{+30.7}$. Furthermore, we find that the average multiplicity of simulated planetary systems containing a giant planet is 2.2 despite our imposed dynamical stability criteria which disfavor giant planets in multiplanet systems. This suggests that M dwarf systems containing both small and giant planets could be discovered with SPIRou. Such systems— should they exist—would be crucial for informing planet formation scenarios around M dwarfs.

## 10. The Effect of an Increased Planet Frequency around Late M Dwarfs

Recall that thus far we have assumed that the *Kepler* occurrence rates, computed using stars with $T_{\text{eff}} \gtrsim 3200$ K (i.e., spectral types of M4.5 and earlier; Luhman et al. 2003), are also applicable to the later M dwarfs in our stellar sample. However, there are some lines of evidence that suggest that the cumulative planet occurrence rate will increase toward later M dwarfs. For example, planet formation models around very low-mass stars predict many small planets ($r_p \sim 1 R_⊕$) at short orbital periods (e.g., Alibert et al. 2013; Alibert & Benz 2017), and the detection of seven Earth-sized planets around the ultracool dwarf TRAPPIST-1 from a small sample of ultracool dwarfs (Gillon et al. 2017; Luger et al. 2017) are suggestive of an increased cumulative planet occurrence rate around late M dwarfs compared to early M dwarfs. If this is true, then two possible outcomes on the SPIRou SLS-PS planet yield are imaginable. Either the predicted planet yield will increase as a result of the greater number of potential planets to detect, or the predicted planet yield will decrease because adding additional planets will contribute to the observed RV rms, thus deterring our ability to detect individual planets. Because the latter effect will modify our detection sensitivity, we cannot simply estimate the resulting planet yield from the product of a scaled-up planet occurrence rate with our nominal detection sensitivity from Figure 13. Instead, to address this caveat in a simplified way, we simulate a new version of the SLS-PS in which we artificially increase the planet occurrence rates by a

<footnote n="16"><sup>16</sup> https://github.com/r-cloutier/SLSPS_Simulations</footnote>
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Top: detection sensitivity as a function of minimum planet mass for the nominal fiducial version of the SLS-PS containing 2.39 planets per star and for a modified version of the SLS-PS that is nearly identical except for the increased cumulative planet occurrence rate of 3.56 planets per star. Bottom: predicted planet yield as a function of minimum planet mass for the two aforementioned survey versions. Over the range of $m_p \sin i \in [0.4, 15] M_\oplus$, a factor of 1.3 more planets are detected in the survey with the increased cumulative planet occurrence rate.</figcaption></figure>

scaling factor and recompute the SLS-PS detection sensitivity and planet yield. All inputs in this simulation other than the planet occurrence rates are identical to the fiducial survey presented throughout this paper.

To increase the planet occurrence rates, we will use a simple scaling factor. That is because in practice, the planet occurrence rates around late M dwarfs are not well-characterized so we adopt the same *Kepler* occurrence rates $f(P, r_p)$ but increase it by a factor of 2 (i.e., a revised cumulative planet occurrence rate of $5 \pm 0.4$ planets per M dwarf). We then conduct the simulated SLS-PS identically as before. However, the dynamical stability considerations (see Section 5.2) that restrict the types of simulated multiplanet systems can have a more pertinent effect when the number of sampled planets per planetary system is doubled. We find that in practice these considerations cause the resulting planet occurrence rates to only be increased by a factor of $\sim 1.5$ instead of 2, such that the resulting cumulative planet occurrence rate is 3.56 rather than 5 planets per star.

Figure 16 compares the detection sensitivity and planet yield recovered by the fiducial version of the SLS-PS—with 2.4 planets per star—with those obtained after increasing the cumulative planet occurrence rate to 3.56 planets per star. The detection sensitivity as a function of $m_p \sin i$ is only slightly modified by the addition of, on average, 1.2 planets per planetary system. That is, the fiducial version of the survey presented in this paper only has a slightly greater detection

sensitivity to planets with $m_p \sin i \in [0.4, 20] M_\oplus$ because the fewer planets on average per planetary system contribute to a lower RV rms than when the planet occurrence is increased. Explicitly, the average detection sensitivity over the range of $m_p \sin i$ considered in Figure 16 is 47.5% for our fiducial survey version compared to 41.9% when the cumulative planet occurrence rate is increased.

The net effect in the new version of the SLS-PS of having more planets per planetary system and a comparable average detection sensitivity is that more planets are detected overall. In each $m_p \sin i$ bin shown in Figure 16, more planets are detected when the cumulative planet occurrence rate is increased, resulting in 1.3 times more SPIRou planet detections ($110.7_{-12.6}^{+28.5}$ planets) than when assuming the nominal *Kepler* planet occurrence rates.

Despite there being evidence for an increased planet occurrence rate $f$ around late M dwarfs than around the early-to-mid M dwarfs observed by *Kepler*, consideration of an SLS-PS in which the cumulative planet occurrence rate is increased is not intended to be interpreted as a new estimate of the SLS-PS planet yield. This is because our simplified methodology for scaling-up $f$ is not intended to replicate the true $f$ around late M dwarfs. We have adopted a simple scaling of $f$ as measurements of $f$ around late M dwarfs are currently poorly constrained. Furthermore, there exist some lines of empirical and theoretical evidence (Gillon et al. 2017; Pan & Schlichting 2017) where the fraction of multiplanetary systems that form low-order mean-motion resonant chains is larger around late M dwarfs than is seen in the *Kepler* sample (Lissauer et al. 2011; Fabrycky et al. 2014). In this case, the Keplerian RV modeling used in this study would not be applicable, which may affect the resulting planet yield. However, this exercise is intended to demonstrate the nature of the effect on our detection sensitivity and ultimately on the SPIRou planet yield if the cumulative planet occurrence rate increases toward later M dwarfs. Fortunately, the effect seems to be a positive one for the detection of small planets around late M dwarfs with SPIRou.

# 11. Measuring the RV Planet Frequency

## 11.1. Recovering Planet Frequency

With consistent RV monitoring of the target stars in the SLS-PS, the resulting detections from the survey will be able to provide independent constraints on the occurrence rate of RV planets around M dwarfs. Of course, this calculation has previously been done using either RV or transit survey data (e.g., Bonfils et al. 2013a; Dressing & Charbonneau 2015), although those studies were limited to early M dwarfs, whereas SPIRou is uniquely designed to study mid-to-late M dwarfs whose planet occurrence rates are less certain (Demory et al. 2016). Furthermore, SPIRou is expected to uncover a larger set of planet detections than the HARPS M dwarf sample presented in Bonfils et al. (2013a), and thus provide stronger constraints on the full M dwarf planet occurrence rate as a function of minimum planetary mass.

Here we wish to estimate the precision with which we expect to measure the planet frequency based on the expected results of the SLS-PS. The planet frequency $f$ differs somewhat from the planetary occurrence rate in that the planetary frequency does not represent the number of a particular type of planet per host star but instead is the fraction of stars that host a particular type of planet and is therefore only defined on the closed
interval $f \in [0, 1]$. To compute the planet frequency as a function of $P$ (or similarly $S$) and $m_p \sin i$, we will adopt the formalism from Carson et al. (2006) and Lafrenière et al. (2007a). In our survey of $N$ stars denoted by $j = 1, \dots, N$, we wish to compute the fraction of M dwarfs that host a planet within a particular range of orbital periods and minimum masses, $f \to f(P, m_p \sin i)$. First, we must estimate the probability that such a planet will be detected around the $j$th star $p_j$, based on the star's known detection sensitivity (see Section 8). The probability that a particular planet will be detected orbiting the $j$th star is then $fp_j$, whereas the probability of a non-detection is $1 - fp_j$. Given the resulting planet yield from our simulated SLS-PS, we can identify around which stars a particular planet is detected. For a given range of $P$ and $m_p \sin i$, we denote planet detections within that range around the $j$th star by $d_j$, which equals 1 if such a planet is detected and 0 otherwise. Now we can write down the likelihood of our planet detections given $f$ as

$$ \mathcal{L}(d_j|f) = \prod_{j=1}^{N} (1 - fp_j)^{1-d_j}(fp_j)^{d_j}. \quad (15) $$

In order to compute the value of $f(P, m_p \sin i)$ in various $P$, $m_p \sin i$ bins, we will invoke Bayes' theorem,

$$ P(f|d_j) = \frac{\mathcal{L}(d_j|f)p(f)}{\int_0^1 \mathcal{L}(d_j|f)p(f)df}, \quad (16) $$

where $p(f)$ is the prior probability of measuring $f$ and $P(f|d_j)$ is the posterior PDF of measuring a frequency $f$ given the observations $d_j$. To compute $f(P, m_p \sin i)$ over the full grid of $P$ and $m_p \sin i$, the above formalism is applied independently to each logarithmically spaced bin in $P$ and $m_p \sin i$. The fraction of M dwarfs with a particular planet is defined as the MAP value of the posterior PDF with its uncertainties characterized by the 16th and 84th percentiles of the distribution.

Before computing $f(P, m_p \sin i)$ from the results of our simulated SLS-PS, we must first assign planet detections around each star in the target sample to integer values rather than the statistical averages used to present the results of the full survey. To do so, we round the number of detected planets in each $P$, $m_p \sin i$ bin to the nearest integer for each star individually and for the full distribution of detected planets from the SLS-PS (see Figure 14). Rounding to integers alters the total number of detected planets so we use a small multiplicative correction factor on each star's detected planet population such that we recover the correct total planet yield of 85 planets. By combining the results of all realizations, the number of planets detected in each $P$, $m_p \sin i$ bin over the full simulated survey always exceeds the rounded average value shown in Figure 14. To account for this when assigning the 85 planet detections, we assign planet detections to stars based on their sensitivity to that particular type of planet. Specifically, in each $P$, $m_p \sin i$ bin in which we detect at least one planet, we identify the subset of stars which have at least one planet detection in that bin and sort those stars by their detection sensitivities within that bin. We then select the $n$ stars with the highest detection sensitivity, where $n$ is the rounded total number of planets detected within that bin. Those stars have $d_j$ set to unity within that bin, whereas $d_j = 0$ for all remaining stars. We note that by this routine stars with the highest detection sensitivities are frequently chosen so we limit the

number of planets that can be detected around a single star to be $\le 3$, the maximum number of planets detected by our automated detection algorithm in the simulated SLS-PS. We also note that by selecting stars with the highest detection sensitivity when assigning planet detections, we are maximizing the likelihood (Equation (15)) and consequentially computing the maximally constrained planet frequency values.

The planet frequency derived using the predicted SPIRou planet detections spanning $P \in [0.5, 200]$ days, $S \in [0.01, 100]$ $S_⊕$, and $m_p \sin i \in [0.4, 15]$ $M_⊕$ are shown in Figure 17. For bins in which the MAP value of the $f$ posterior PDF is non-zero, we report the MAP value and its 16th and 84th percentiles in Figure 17. For bins in which we detect only a small number of planets or have a low detection sensitivity, we find an MAP $f = 0$ and report the 68th percentile as an upper limit.

By considering the SPIRou planet detections and our detection sensitivity within various ranges of planet parameters, we can estimate the frequency of such planets and the precision of the measurement. We remind the reader that it is the precision of the measurement that is meaningful here as the MAP values of $f$ are simply the result of combining the known Kepler occurrence rates, an empirical mass-radius relation, and our detection sensitivity, and therefore does not provide any new information regarding the planet frequencies themselves. However, our adopted formalism (see Equations (15) and (16)) computes the planet frequency over a specified range of planetary parameters. Recall that the planet frequency is not equivalent to the planet occurrence rate when the occurrence rate is greater than unity. Therefore, over the full range of $P$, $S$, and $m_p \sin i$ presented in Figure 7, we cannot use this formalism to compute the cumulative SPIRou planet occurrence rate, which is $> 1$ at $\sim 2.5$ planets per star. Instead, we must take the ratio of the SPIRou planet detection map shown in Figure 14—rounded to integer values—with the sensitivity map shown in Figure 13. The resulting map depicts the SPIRou-derived planet frequency as a function of $P$, $S$, and $m_p \sin i$, which we can then integrate over to estimate the cumulative planet occurrence rate of $1.8 \pm 0.2$ planets per M dwarf.

## 11.2. Measuring $\eta_⊕$

Using the formalism from Section 11.1 to compute the planet frequency, we can estimate the frequency of any subset of planets. Of particular interest is the frequency of potentially habitable planets around M dwarfs, $\eta_⊕$. We will define potentially habitable—Earth-like—planets as those with $m_p \in [1, 5]$ $M_⊕$ and within the fiducial HZ period limits defined by the equations in Kopparapu et al. (2013) for “water loss” and “maximum greenhouse.” Here, the absolute planet mass is inferred from the detected population of minimum planet masses by correcting for the geometrical effect of randomly oriented orbits following a geometrical distribution. The upper limit on $m_p \sim 5 M_⊕$ approximately corresponds to the expected mass of a $1.5 R_⊕$ planet, which marks the approximate radius boundary between rocky/Earth-like and gaseous planets (e.g., Valencia et al. 2013; Lopez & Fortney 2014; Fulton et al. 2017). By our definition, we expect to detect nine potentially habitable planets in the full SLS-PS of 100 stars. Based on our derived detection sensitivity to such planets, we measure $\eta_⊕ = 0.28_{-0.07}^{+0.12}$, derived from its posterior PDF shown in Figure 18.
<figure><img src="image_20.png" /></figure>

**Figure 17.** Coarsely binned maps of the RV planet frequency derived from the SLS-PS as a function of minimum planet mass and orbital period (left) and insolation (right). The dashed vertical lines in the insolation panel indicate the approximate “water-loss” and “maximum-greenhouse” insolation limits of the HZ from Kopparapu et al. (2013). The dashed-dotted vertical lines indicate the more conservative “recent-Venus” and “early-Mars” HZ limits (Kopparapu et al. 2013). Overplotted in each bin are the planet frequency posterior PDFs derived from the SLS-PS planet detections (Figure 14) and sensitivity (Figure 13) using Equation (16). When the planet frequency PDF is consistent with 0, we report the 68th percentile as an upper limit. When a non-zero planet frequency is detected, we report the MAP value along with the 16th and 84th percentiles.

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

**Figure 18.** Probability density function of the RV value of $\eta_{\oplus}$ derived from the simulated SLS-PS. Here, Earth-like HZ planets are defined according to the “water-loss” and “maximum-greenhouse” HZ limits defined in Kopparapu et al. (2013, K13) and have absolute planet masses $\le 5 M_{\oplus}$. The inner shaded region corresponds to 16th and 84th percentiles whereas the outer regions mark the 2nd and 98th percentiles.

Focusing solely on late M dwarfs (M5–M9) then, we expect to detect $3.4_{-1.4}^{+3.2}$ potentially habitable planets. As a result we measure a new $\eta_{\oplus} = 0.29_{-0.10}^{+0.24}$ whose measurement uncertainty is effectively double relative to the uncertainty on the value of $\eta_{\oplus}$ derived from the full SPIRou input catalog. This is because only 37/100 stars from the input catalog are classified as late M dwarfs with spectral type later than M5. The result of the decreased sample size and the, on average, lower detection sensitivity around dim late M dwarfs is a greater uncertainty on $\eta_{\oplus}$ around late M dwarfs. However, recall that existing estimates of $\eta_{\oplus}$ around M dwarfs have been limited to early-

to-mid M dwarfs, making the SPIRou estimate of $\eta_{\oplus}$ around late M dwarfs potentially the first of its kind.

## 12. Direct Imaging of Nearby Planetary Systems

Assessing the habitability of exoplanets relies heavily on probing the planet's atmospheric conditions. One potential avenue for studying the atmospheres of non-transiting exoplanets is via high-contrast imaging wherein photons from the spatially resolved planet are directly detected following the suppression of quasi-static speckles associated with the bright host star. Various observational techniques such as adaptive optics/coronagraphy and post-processing techniques (e.g., ADI, Marois et al. 2006; LOCI, Lafrenière et al. 2007b; and KLIP, Soummer et al. 2012) have enabled the direct detection of a number of exoplanets via high-contrast imaging. However, the planet population for which this technique is currently amenable is limited to self-luminous substellar objects at large angular separations from their host star, e.g., young gas giants on wide orbits with planet-to-star contrasts of $O(10^{-4})$. However, the large apertures on board the upcoming generation of ELTs will offer sufficiently high spatial and spectral resolution to reach the small planet-to-star contrasts in the near-IR ($\sim O(10^{-7})$) required to directly image a small number of small HZ planets around the closest stars. Time-resolved rotational color variations of small planets may even permit the detection of the number, reflectance spectra, sizes, and longitudinal positions of major surface features such as continents and/or liquid oceans, and the measurement of cloud properties (Ford et al. 2001; Cowan et al. 2009; Fujii et al. 2010, 2011; Cowan & Strait 2013).

To obtain a set of direct images of HZ exoplanets, we must first find the closest habitable worlds. Of particular focus are M dwarfs in the solar neighborhood because of their abundance and the favorable contrasts of their HZ planets compared to HZ planets around Sun-like stars (Crossfield 2013). The closest HZ exoplanet has already been discovered around Proxima Centauri (1.3 pc; Anglada-Escudé et al. 2016), but many more
M dwarf HZ planets likely remain undetected within ≲10 pc. The SLS-PS will uncover many of these systems and potentially with a lesser total observation time per planetary system than that which is required to detect M dwarf planetary systems using optical velocimeters. To identify the subset of SPIRou detections amenable to direct imaging, we adopt the near-IR contrast performance expected for a number of dedicated imagers and compare the detected SPIRou planet population to these contrast curves in the planet's projected separation/contrast space.

With the upcoming ELT imagers optimized at near-IR wavelengths, targeted planets are observed in reflected light such that the expected planet-to-star contrast is

$$C = 1.81 \times 10^{-7} A \left( \frac{r_p}{1 R_{\oplus}} \right)^2 \left( \frac{a}{0.1 \text{ au}} \right)^{-2}, \quad (17)$$

where $r_p$ is the planet's radius, $a$ is the semimajor axis, and $A$ is the geometric albedo, which we assume is 0.3 for all detected planets in the simulated SLS-PS. At such low planet-to-star contrasts, no existing high-contrast imager is presently capable of imaging HZ planets. However, there are proposed techniques to achieve such small near-IR contrasts, which involve coupling high-contrast imaging capabilities to spatially resolve the targeted planet followed by the use of high-dispersion spectroscopy to filter out the stellar component as a result of the differentially Doppler-shifted planet and stellar spectra (e.g., Snellen et al. 2015; Lovis et al. 2017). This technique was recently used with CRIRES on the VLT to measure carbon monoxide (Snellen et al. 2010; Brogi et al. 2012; de Kok et al. 2013) and water (Birkby et al. 2013) in the atmospheres of hot Jupiters. This technique has also become considerably more topical since the discovery of an HZ planet in the closest exoplanetary system—Proxima Centauri—and the prospect of detecting the potential biosignature $O_2$ in this system using the VLT (Lovis et al. 2017).

Similar promising techniques have led to a suite of contrast curve predictions which indicate the types of planets that may be detected at a particular detection significance given the capabilities of the imager and the expected S/N achievable based on the target properties and observational strategy. Of these, we will consider the expected geometric mean of the E-ELT EPICS<fnref n="17" /> (Kasper et al. 2010) and TMT PFI<fnref n="18" /> (Macintosh et al. 2006) $5\sigma$ H-band contrast curves. We also consider the space-based hybrid Lyot coronagraph on board WFIRST (Trauger et al. 2015).

A random subset of planets from the simulated SLS-PS are shown in Figure 19. Comparing this population to the expected performance of various ELT imaging instruments, we expect $46.7_{-6.0}^{+16.0}$ SPIRou planets to be imaged.<fnref n="19" /> Here, we have defined an imagable planet as one whose expected projected angular separation and contrast lie above the geometric mean of the EPICS and PFI $5\sigma$ H-band contrast curves in Figure 19. Note that this definition does not impose

a minimum inner working angle although such a cut at, say, $3\lambda/D$, would decrease the total number of imagable SPIRou planets by a factor of ~2.7. The subset of imagable planets represents ~55% of all SPIRou planets. Among the imagable planets are $13.7_{-4.9}^{+11.5}$ HZ planets and $4.9_{-2.0}^{+4.7}$ HZ planets with $m_p \in [1, 5]$ $M_{\oplus}$, the so-called Earth-like planets. These SPIRou planets along with the known Proxima Centauri b, Ross 128b (Bonfils et al. 2017), and GJ 273b (Astudillo-Defru et al. 2017) will represent the best potential targets for the imaging of small HZ exoplanets with ELTs.

## 13. Comparison of Different Versions of the SLS-PS

Table 2 summarizes the main results of six simulated versions of the SLS-PS including the fiducial version presented in this paper, which we now refer to as the optimized version in Table 2. Brief descriptions and motivations for each additional version of the SLS-PS are given below.

1. **Optimized:** the SLS-PS version used throughout this study containing 100 stars in the input catalog. This version approximately represents the optimal compromise between maximizing detection sensitivity and producing a satisfactory number of planet detections, including a sizable fraction of planets that may be amenable to direct imaging with ELTs.

2. **Closest:** contains the 50 closest stars ($d \lesssim 6.8$ pc) from the optimized input catalog. Here we target the closest M dwarfs to the solar system in an effort to focus observational resources on a small number of target stars, thus maximizing detection sensitivity and the number of detections of small planets that may be imagable with ELTs.

3. **Large:** contains 360 stars in the input catalog, where 360 is the number of targeted stars in the original SPIRou Science Case<fnref n="20" /> proposal from 2013. Targeting a large sample of M dwarfs may result in the greatest planet yield, which is desirable for putting tight constraints on the cumulative planet occurrence rates around M dwarfs and $\eta_{\oplus}$.

4. **Short:** has an identical input catalog to optimized with 100 stars but includes only half of total available time on sky; we obtain half as many RVs per star as in optimized. This experimental setup is undesirable but the simulation is used to demonstrate by how much the SLS-PS detection sensitivity and yield suffer given fewer measurements.

5. **Degraded:** has an identical input catalog and window functions to optimized with 100 stars but imposes a degraded noise floor on the RV measurement precision of $2 \text{ m s}^{-1}$ compared to the $1 \text{ m s}^{-1}$ noise floor assumed in optimized. This survey version is used to assess the impact of a degradation in the RV measurement precision on the detection sensitivity and planet yield.

6. **Dark:** is identical to the optimized survey with 100 stars but whose window functions are not restricted to non-dark time only. Here we include dark-time observations, thus alleviating the strong aliasing of planets with $P = 15$ or 30 days as in optimized.

<footnote n="17">¹⁷ European-ELT-Exoplanet Imaging Camera and Spectrograph, since renamed the Planetary Camera and Spectrograph (PCS).</footnote>

<footnote n="18">¹⁸ Thirty Meter Telescope-Planet Formation Imager.</footnote>

<footnote n="19">¹⁹ Here, the term “imagable” need not correspond exactly to an imaging observation as achieving the small planet/star contrasts exhibited by SPIRou planets will likely require high-contrast imaging coupled with high-dispersion spectroscopy (e.g., Snellen et al. 2015).</footnote>

<footnote n="20">²⁰ http://spirou.irap.omp.eu/content/download/232/1348/file/SPIRou-sciencecase-last.pdf</footnote>
<figure><img src="image_22.png" /><figcaption>**Figure 19.** Same random subset of the simulated SPIRou planets from Figure 15 in the projected angular separation/reflected light contrast plane. The geometric albedo A is set to 0.3 for all planets. Yellow circles highlight detected Earth-like planets ($m_p \in [1, 5] M_⊕$), orange circles highlight the remaining detected HZ planets, and red circles highlight all non-HZ detected planets. Open circles represent planets that remain undetected by SPIRou. Yellow diamonds labelled 1–5 depict known, likely rocky planets at or near the HZ. The size of each planet’s marker is proportional to its radius. The planet population is compared to the geometric mean of the predicted $5\sigma$ H-band contrast curves for EPICS and PFI (solid curve) and the predicted performance for the WFIRST hybrid Lyot Coronagraph (dotted curve). The projected separation is also depicted in units of $\lambda/D$ in the H band ($\lambda = 1.66$ µm) for a D = 30 m telescope.</figcaption></figure>

## 13.1. Optimized: the Optimal Survey Strategy

The *optimized* survey version represents an ideal compromise between (i) achieving sufficient planet detection sensitivity to put meaningful constraints on the planet occurrence rates —including $\eta_⊕$— and (ii) producing a large planet yield including a set of imagable planets and, in particular, Earth-like imagable planets from the SLS-PS. For other RV planet search campaigns with similar science goals to those aforementioned, we advocate for a similar survey strategy to *optimized*. With that said, at the time of writing, the SPIRou input catalog has not been absolutely defined and will likely be altered between the time of these simulations and the beginning of the actual SLS-PS. The results presented in this paper are therefore intended as a guideline to inform how many stars should be included in the SLS-PS.

Here, we detect ~8 Earth-like planets with a sensitivity of ~33.5%. This will result in a constraint on $\eta_⊕$ at a level of precision of ≲45. We also detect ~5 Earth-like planets that may be amenable to direct imaging with ELTs. This is the largest number of imagable Earth-like planets detected with SPIRou compared to any other simulated survey version other than *dark*, whose idealized experimental setup is not feasible with the suite of instruments on the CFHT.

## 13.2. Closest: the Closest M Dwarfs

The *closest* survey version reduces the size of the SPIRou input catalog to 50 stars but maintains the same volume of nights. Specifically, the 50 closest stars from the *optimized* SPIRou input catalog are retained. By focusing on fewer stars with the same amount of total survey time as in the *optimized* survey, more observing time can be dedicated to each star. The result is a higher detection sensitivity to all planets around each star compared to *optimized*. In this way, we are able to detect ~5 Earth-like planets, ~4 of which will be amenable to direct imaging with ELTs due to their close proximity to the solar system.

Although the sensitivity to any given planet is maximized in the *closest* survey version, the stellar input catalog is too small to result in a larger total planet yield than in *optimized*. The small number of planet detections also has a detrimental effect on measuring the cumulative planet occurrence rate around M dwarfs and in particular late M dwarfs.

## 13.3. Large: Lots of Stars

The *large* survey version was originally considered as the tentative survey strategy for the SLS-PS. Its aim is to discover the greatest number of exoplanets by surveying many more stars (360 stars) than any other considered survey version. Targeting
<table><caption>**Table 2**
Overview of SLS-PS Versions</caption><thead><tr><th></th><th>Optimized</th><th>Closest</th><th>Large</th><th>Short</th><th>Degraded</th><th>Dark</th></tr></thead><tbody><tr><td>Fraction of total available time on sky<sup>a</sup></td><td>1</td><td>1</td><td>1</td><td>0.5</td><td>1</td><td>1</td></tr><tr><td>Number of target stars</td><td>100</td><td>50</td><td>360</td><td>100</td><td>100</td><td>100</td></tr><tr><td>Average number of RVs per star</td><td>198.1</td><td>396.0</td><td>62.0</td><td>99.5</td><td>198.1</td><td>198.1</td></tr><tr><td>Median σ<sub>RV</sub> [m s<sup>-1</sup>]</td><td>1.33</td><td>1.21</td><td>1.88</td><td>1.33</td><td>2.52</td><td>1.33</td></tr><tr><td>Median expected K measurement uncertainty [m s<sup>-1</sup>]</td><td>0.19</td><td>0.14</td><td>0.40</td><td>0.25</td><td>0.26</td><td>0.20</td></tr><tr><td>Average detection sensitivity<sup>b</sup> [%]</td><td>44.9 ± 0.5</td><td>53.6 ± 1.1 (1.2)<sup>c</sup></td><td>19.5 ± 0.2 (0.4)</td><td>34.4 ± 0.4 (0.8)</td><td>34.6 ± 0.6 (0.8)</td><td>47.2 ± 0.7 (1.1)</td></tr><tr><td>Average detection sensitivity to HZ planets<sup>d</sup> [%]</td><td>43.1 ± 1.0</td><td>51.2 ± 2.1 (1.2)</td><td>19.6 ± 0.3 (0.5)</td><td>34.8 ± 0.9 (0.8)</td><td>34.0 ± 1.2 (0.8)</td><td>48.1 ± 1.3 (1.1)</td></tr><tr><td>Average detection sensitivity to Earth-like planets<sup>e</sup> [%]</td><td>33.5 ± 1.2</td><td>45.4 ± 2.9 (1.4)</td><td>8.1 ± 0.3 (0.2)</td><td>21.4 ± 1.0 (0.6)</td><td>20.7 ± 1.3 (0.6)</td><td>35.2 ± 1.5 (1.1)</td></tr><tr><td>Average detection sensitivity to imagerable planets [%]</td><td>47.7 ± 0.7</td><td>52.8 ± 1.3 (1.1)</td><td>24.5 ± 0.4 (0.5)</td><td>37.1 ± 0.6 (0.8)</td><td>36.0 ± 0.9 (0.8)</td><td>49.7 ± 1.0 (1.0)</td></tr><tr><td>Average detection sensitivity to imagerable HZ planets [%]</td><td>46.1 ± 1.3</td><td>50.3 ± 2.3 (1.1)</td><td>27.9 ± 1.0 (0.6)</td><td>37.8 ± 1.2 (0.8)</td><td>35.2 ± 1.6 (0.8)</td><td>51.5 ± 1.8 (1.1)</td></tr><tr><td>Average detection sensitivity to imagerable Earth-like planets [%]</td><td>33.7 ± 1.6</td><td>42.5 ± 3.1 (1.3)</td><td>11.4 ± 1.0 (0.3)</td><td>20.4 ± 1.2 (0.6)</td><td>18.1 ± 1.7 (0.5)</td><td>35.7 ± 2.1 (1.1)</td></tr><tr><td>Total planet yield</td><td>85.3<sup>+29.3</sup><sub>-12.4</sub></td><td>50.6<sup>+15.7</sup><sub>-6.5</sub> (0.6)</td><td>142.7<sup>+77.7</sup><sub>-29.5</sub> (1.7)</td><td>65.7<sup>+27.8</sup><sub>-11.1</sub> (0.8)</td><td>65.2<sup>+27.1</sup><sub>-11.0</sub> (0.8)</td><td>88.7<sup>+30.6</sup><sub>-13.0</sub> (1.0)</td></tr><tr><td>Total yield of HZ planets</td><td>20.0<sup>+16.8</sup><sub>-7.2</sub></td><td>12.1<sup>+10.1</sup><sub>-4.3</sub> (0.6)</td><td>35.3<sup>+29.5</sup><sub>-12.6</sub> (1.8)</td><td>16.4<sup>+13.7</sup><sub>-5.9</sub> (0.8)</td><td>15.9<sup>+13.3</sup><sub>-5.7</sub> (0.8)</td><td>22.9<sup>+19.1</sup><sub>-8.2</sub> (1.1)</td></tr><tr><td>Total yield of Earth-like planets</td><td>8.1<sup>+7.6</sup><sub>-3.2</sub></td><td>5.4<sup>+5.0</sup><sub>-2.1</sub> (0.7)</td><td>7.6<sup>+7.2</sup><sub>-3.0</sub> (0.9)</td><td>5.2<sup>+4.9</sup><sub>-2.1</sub> (0.6)</td><td>5.0<sup>+4.8</sup><sub>-2.0</sub> (0.6)</td><td>8.6<sup>+8.1</sup><sub>-3.4</sub> (1.1)</td></tr><tr><td>Total yield of imagerable planets</td><td>46.7<sup>+16.0</sup><sub>-6.0</sub></td><td>33.5<sup>+10.4</sup><sub>-3.8</sub> (0.7)</td><td>38.3<sup>+20.9</sup><sub>-7.0</sub> (0.8)</td><td>36.5<sup>+15.4</sup><sub>-5.4</sub> (0.8)</td><td>34.4<sup>+14.3</sup><sub>-5.1</sub> (0.7)</td><td>48.3<sup>+16.7</sup><sub>-6.2</sub> (1.0)</td></tr><tr><td>Total yield of imagerable HZ planets</td><td>13.7<sup>+11.5</sup><sub>-4.9</sub></td><td>9.7<sup>+8.2</sup><sub>-3.5</sub> (0.7)</td><td>8.2<sup>+6.9</sup><sub>-2.9</sub> (0.6)</td><td>11.3<sup>+9.5</sup><sub>-4.0</sub> (0.8)</td><td>10.1<sup>+8.5</sup><sub>-3.6</sub> (0.7)</td><td>15.5<sup>+12.9</sup><sub>-5.5</sub> (1.1)</td></tr><tr><td>Total yield of imagerable Earth-like planets</td><td>4.9<sup>+4.7</sup><sub>-2.0</sub></td><td>4.0<sup>+3.7</sup><sub>-1.6</sub> (0.8)</td><td>1.5<sup>+1.4</sup><sub>-0.6</sub> (0.3)</td><td>3.0<sup>+2.8</sup><sub>-1.2</sub> (0.6)</td><td>2.5<sup>+2.5</sup><sub>-1.0</sub> (0.5)</td><td>5.3<sup>+5.0</sup><sub>-2.1</sub> (1.1)</td></tr><tr><td>Cumulative planet occurrence rate<sup>b</sup> [planets per star]</td><td>1.8 ± 0.2</td><td>1.8 ± 0.3</td><td>1.8 ± 0.2</td><td>1.8 ± 0.2</td><td>1.8 ± 0.2</td><td>1.8 ± 0.2</td></tr><tr><td>Frequency of Earth-like planets, η<sub>⊕</sub></td><td>0.29<sup>+0.13</sup><sub>-0.07</sub></td><td>0.26<sup>+0.16</sup><sub>-0.08</sub></td><td>0.32<sup>+0.18</sup><sub>-0.09</sub></td><td>0.30<sup>+0.18</sup><sub>-0.09</sub></td><td>0.26<sup>+0.18</sup><sub>-0.09</sub></td><td>0.27<sup>+0.11</sup><sub>-0.07</sub></td></tr></tbody></table>

**Notes.**

a For a notional survey duration of ~300 nights over ~3 years.

b Over the range of $P \in [0.5, 200]$ days and $m_p \sin i \in [0.4, 20]$ $M_\odot$.

c Numbers in parentheses indicate the fractional value of the corresponding quantity relative to the optimized survey version.

d Based on the “water-loss” and “maximum greenhouse” limits of the HZ from Kopparapu et al. (2013).

e Earth-like planets are defined by the “water-loss” and “maximum greenhouse” HZ limits (Kopparapu et al. 2013) and have $m_p \in [1, 5]$ $M_\odot$.

so many stars comes at the expense of a reduced detection sensitivity per star and a particularly low detection sensitivity to Earth-like planets, which seem to require ≳150 RV measurements to detect (e.g., Astudillo-Defru et al. 2017). Although the overall planet yield in the *large* survey version is high, many of the most interesting systems—Earth-like planets that may be imagerable with ELTs—will remain largely undetected. For example, we only detect ≲2 imagerable Earth-like planets in *large* compared to the ~5 in *optimized*. Furthermore, the increase in precision with which the cumulative planet occurrence rate can be measured in *large* is only marginal, and in our opinion, not worth the small yield of imagerable Earth-like planets.

## 13.4. Short: the Effect of Fewer Observations

The short survey version features an identical input catalog to the optimized version but with half as many RV observations per star. This survey strategy is useful to characterize the loss in detection sensitivity and planet yield if fewer than the total number of possible measurements are obtained throughout the campaign. The loss in detection sensitivity—and hence in planet yield—compared to optimized evolves approximately as $\sqrt{n_{\text{obs}}}$ for all but the smallest, Earth-like planets. The loss in sensitivity for Earth-like planets is worsened by the small number of measured RVs in small. The rough scaling of

detection sensitivity with $\sqrt{n_{\text{obs}}}$ is the result of the use of our GP regression activity modeling to model non-white noise in active time series, and that modeling performs well on the majority of applicable systems (see Figure 10).

## 13.5. Degraded: a Degraded RV Measurement Precision

In the optimized version of the survey, we had assumed that SPIRou will operate with a long-term RV precision of $\sigma_{RV} = 1$ m s⁻¹. This was imposed as an RV noise floor on all stars for which we are likely to be able to achieve a photon-noise-limited RV uncertainty of < 1 m s⁻¹. However, given that the long-term RV stability of SPIRou has yet to be tested on sky, it is conceivable that a degraded level of RV precision, >1 m s⁻¹, may instead be realized. Here, we repeat the simulation of the optimized survey but increase the RV noise floor to 2 m s⁻¹.

Similarly to that in short, the degradation in detection sensitivity approximately scales as $1/\sqrt{\sigma_{RV}}$. The results for degraded are therefore closely related to the results for short.

## 13.6. Dark: No Window Function Restrictions

The dark survey version relaxes the assumption used in optimized that SPIRou observations may only be obtained during non-dark time. Although such a scenario is unlikely to be realized for SPIRou on CFHT—given the scheduling of
other instruments on the telescope—window functions that do include dark-time observations may be obtained with other near-IR velocimeters like NIRPS on the ESO 3.6 m telescope at La Silla (Bouchy et al. 2017). This experimental setup represents a best-case scenario for SPIRou and illustrates how many more planets can be detected with dark-time observations including a number of HZ planets around M2-4 dwarfs.

We can detect ~3 more planets in total compared to *optimized*. Most of these planets lie within the HZ. The number of imagable planets detected is only increased by ~2 new planets with fewer than one additional imagable Earth-like planet in *dark* compared to *optimized*. These modest increases in the planet yield result in only slight improvements to the measured cumulative planet occurrence rate and $\eta_{\oplus}$.

R.C. thanks the Canadian Institute for Theoretical Astrophysics for use of the Sunnyvale computing cluster throughout this work. R.C. is partially supported in this work by the National Science and Engineering Research Council of Canada. J.F.D. thanks the IDEX initiative of Université Fédérale Toulouse Midi-Pyrénées (UFTMiP) for generous fund allocation without which SPIRou could not have been built. J.F.D. also thanks the ERC for funding the NewWorlds project focused on SPIRou-related science. X. Dumusque acknowledges the Society in Science—The Branco Weiss Fellowship for its financial support.

## Appendix

### Computing $\sigma_K$ from the Fisher Information Matrix

As discussed in Section 7.3, we vet putative planet detections by insisting that all bona fide planet detections have an RV semi-amplitude detection significance of at least $3\sigma$; $K/\sigma_K \ge 3$. In order to estimate the precision with which a planet's semi-amplitude $K$ can be measured in a given RV time series $y_k = y(t_k)$, we must estimate the semi-amplitude measurement uncertainty $\sigma_K$ from the Fisher information matrix. The Fisher information matrix $B$ encodes the amount of information about $K$ contained within the data set $y_k$ and is computed analytically from the ln likelihood given in Equation (10) under some simplifying assumptions. Namely, we treat each detected planet in a given system individually, thus restricting the Keplerian model $\mu_k = \mu(t_k)$ to a fixed number of parameters describing a single Keplerian orbital solution. Second, we adopt the simplifying assumption that the planet is on a circular orbit such that only three model parameters (i.e., $\theta = \{P, T_0, K\}$) need be constrained by the time series. When considering planets individually, we must absorb all noise sources from unmodeled RV activity and additional planets into an effective RV uncertainty $\sigma_{\text{eff}}$ equal to the rms of the RVs after the removal of the planet's Keplerian model. Lastly, we assume that the noise properties of the time series are Gaussian distributed such that we can write the ln likelihood as

$$ \ln \mathcal{L} = - \frac{1}{2\sigma_{\text{eff}}^2} \sum_{k=1}^{n_{\text{obs}}} (y_k - \mu_k)^2 + c, \quad (18) $$

where $c$ is a constant that is independent of the model parameters in $\theta$ and $n_{\text{obs}}$ is the number of observations in the time series. Recall that for a circularized planet, the Keplerian

orbital solution reduces to a simple sinusoid which we write as

$$ \mu_k = -K \sin \phi_k, \quad (19) $$

where $\phi_k = 2\pi(t_k - T_0)/P$ and $T_0$ represents the epoch of inferior conjunction for a circularized planet. The elements of the matrix $B$ are then computed via

$$ B_{ij} = - \frac{\partial^2 \ln \mathcal{L}}{\partial \theta_i \partial \theta_j}. \quad (20) $$

Taking $\theta = \{P, T_0, K\}$, the symmetric Fisher matrix takes the form

$$ B = \begin{bmatrix} B_{P,P} & B_{P,T_0} & B_{P,K} \\ B_{T_0,P} & B_{T_0,T_0} & B_{T_0,K} \\ B_{K,P} & B_{K,T_0} & B_{K,K} \end{bmatrix} \quad (21) $$

and contains six independent terms. Here, we explicitly compute the analytical forms of each independent element of $B$ calculated using Equations (18)–(20). The first partials of the ln likelihood with respect to each of the model parameters in $\theta$ are

$$ \frac{\partial \ln \mathcal{L}}{\partial P} = - \left( \frac{\pi K}{P^2 \sigma_{\text{eff}}^2} \right) \sum_{k=1}^{N} (K \sin 2\phi_k - 2y_k \cos \phi_k)(t_k - T_0), \quad (22) $$

$$ \frac{\partial \ln \mathcal{L}}{\partial T_0} = - \left( \frac{\pi K}{P \sigma_{\text{eff}}^2} \right) \sum_{k=1}^{N} (K \sin 2\phi_k - 2y_k \cos \phi_k), \quad (23) $$

$$ \frac{\partial \ln \mathcal{L}}{\partial K} = \left( \frac{1}{\sigma_{\text{eff}}^2} \right) \sum_{k=1}^{N} (K \sin^2 \phi_k - y_k \sin \phi_k). \quad (24) $$

We are now in a position to compute the six independent elements of the Fisher matrix using Equation (20) and the first partials given in Equations (22)–(24):

$$ \begin{aligned} B_{P,P} &= -\frac{\partial}{\partial P} \left( \frac{\partial \ln \mathcal{L}}{\partial P} \right) \\ &= \frac{2\pi K}{P^3 \sigma_{\text{eff}}^2} \left[ \sum_{k=1}^{N} \left( 2y_k \cos \phi_k - \frac{2\pi K(t_k - T_0)}{P} \cos 2\phi_k - \frac{2\pi y_k(t_k - T_0)}{P} \sin \phi_k - K \sin 2\phi_k \right) (t_k - T_0) \right], \end{aligned} \quad (25) $$

$$ \begin{aligned} B_{T_0,P} &= -\frac{\partial}{\partial T_0} \left( \frac{\partial \ln \mathcal{L}}{\partial P} \right) \\ &= \frac{\pi K}{P^2 \sigma_{\text{eff}}^2} \left[ \sum_{k=1}^{N} \left( 2y_k \cos \phi_k - \frac{4\pi K(t_k - T_0)}{P} \cos 2\phi_k - \frac{4\pi y_k(t_k - T_0)}{P} \sin \phi_k - K \sin 2\phi_k \right) \right], \end{aligned} \quad (26) $$

$$ \begin{aligned} B_{K,P} &= -\frac{\partial}{\partial K} \left( \frac{\partial \ln \mathcal{L}}{\partial P} \right) = \frac{2\pi}{P^2 \sigma_{\text{eff}}^2} \\ &\quad \times \left[ \sum_{k=1}^{N} (K \sin 2\phi_k - y_k \cos \phi_k)(t_k - T_0) \right], \end{aligned} \quad (27) $$
$$B_{T_0, T_0} = -\frac{\partial}{\partial T_0} \left( \frac{\partial \ln \mathcal{L}}{\partial T_0} \right) = \frac{4\pi^2 K}{P^2 \sigma_{\text{eff}}^2} \times \left[ \sum_{k=1}^{N} (-K \cos 2\phi_k - y_k \sin \phi_k) \right], \quad (28)$$

$$B_{K,T_0} = -\frac{\partial}{\partial K} \left( \frac{\partial \ln \mathcal{L}}{\partial T_0} \right) = \frac{2\pi}{P \sigma_{\text{eff}}^2} \times \left[ \sum_{k=1}^{N} (K \sin 2\phi_k - y_k \cos \phi_k) \right], \quad (29)$$

$$B_{K,K} = -\frac{\partial}{\partial K} \left( \frac{\partial \ln \mathcal{L}}{\partial K} \right) = -\frac{1}{\sigma_{\text{eff}}^2} \left[ \sum_{k=1}^{N} \sin^2 \phi_k \right]. \quad (30)$$

Using the above expressions to compute the elements of $B$, we can then compute the covariance matrix $C$ of the model parameters in $\theta$ via $C = |B^{-1}|$. The diagonal elements of the $3 \times 3$ matrix $C$ are the estimated measurement variances of the three model parameters $\theta$. Therefore, the measurement uncertainty of the planet's semi-amplitude is $\sigma_K = \sqrt{C_{K,K}}$.

## ORCID iDs

Ryan Cloutier https://orcid.org/0000-0001-5383-9393

Claire Moutou https://orcid.org/0000-0002-2842-3924

Andrew Cumming https://orcid.org/0000-0002-6335-0169

Xavier Dumosque https://orcid.org/0000-0002-9332-2011

## References

* Aigrain, S., Pont, F., & Zucker, S. 2012, MNRAS, 419, 3147

* Alibert, Y., & Benz, W. 2017, A&A, 598, L5

* Alibert, Y., Carron, F., Fortier, A., et al. 2013, A&A, 558, A109

* Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., & O'Neil, M. 2015, ITPAM, 38 arXiv:1403.6015

* Anglada-Escudé, G., Amado, P. J., Barnes, J., et al. 2016, Natur, 536, 437

* Angus, R., Morton, T., Aigrain, S., Foreman-Mackey, D., & Rajpaul, V. 2017, arXiv:1706.05459

* Arlot, S., & Celisse, A. 2010, Statist. Surv., 4, 40

* Artigau, É., Kouach, D., Donati, J.-F., et al. 2014, Proc. SPIE, 9147, 15

* Astudillo-Defru, N., Bonfils, X., Delfosse, X., et al. 2015, A&A, 575, A119

* Astudillo-Defru, N., Forveille, T., Bonfils, X., et al. 2017, A&A, 602, A88

* Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 1998, A&A, 337, 403

* Barnes, J. R., & Collier Cameron, A. 2001, MNRAS, 326, 950

* Barnes, J. R., James, D. J., & Collier Cameron, A. 2004, MNRAS, 352, 589

* Barnes, R., & Greenberg, R. 2006, ApJL, 647, L163

* Berdyugina, S. V. 2005, LRSP, 2, 8

* Berta-Thompson, Z. K., Irwin, J., Charbonneau, D., et al. 2015, Natur, 527, 204

* Birkby, J. L., de Kok, R. J., Brogi, M., et al. 2013, MNRAS, 436, L35

* Bonfils, X., Almenara, J. M., Jocou, L., et al. 2015, Proc. SPIE, 9605, 96051L

* Bonfils, X., Astudillo-Defru, N., Díaz, R., et al. 2017, arXiv:1711.06177

* Bonfils, X., Delfosse, X., Udry, S., et al. 2013a, A&A, 549, A109

* Bonfils, X., Lo Curto, G., Correia, A. C. M., et al. 2013b, A&A, 556, A110

* Bonfils, X., Mayor, M., Delfosse, X., et al. 2007, A&A, 474, 293

* Bouchy, F. & Sophie Team 2006, in 10th Anniversary of 51 Peg-b: Status of and Prospects for hot Jupiter studies, ed. L. Arnold, F. Bouchy, & C. Moutou (Paris: Frontier Group), 319

* Bouchy, F., Doyon, R., Cloutier, R., & Valencia, D. 2017, Msngr, 169, 21

* Boyajian, T. S., von Braun, K., van Belle, G., et al. 2012, ApJ, 757, 112

* Brogi, M., Snellen, I. A. G., de Kok, R. J., et al. 2012, Natur, 486, 502

* Carson, J. C., Eikenberry, S. S., Smith, J. J., & Cordes, J. M. 2006, AJ, 132, 1146

* Cloutier, R., Astudillo-Defru, N., Doyon, R., et al. 2017a, A&A, 608, A35

* Cloutier, R., Doyon, R., Menou, K., et al. 2017b, AJ, 153, 9

* Cloutier, R., Tamayo, D., & Valencia, D. 2015, ApJ, 813, 8

* Correia, A. C. M., Couetdic, J., Laskar, J., et al. 2010, A&A, 511, A21

* Cosentino, R., Lovis, C., Pepe, F., et al. 2012, Proc. SPIE, 8446, 84461V

* Cowan, N. B., Agol, E., Meadows, V. S., et al. 2009, ApJ, 700, 915

* Cowan, N. B., & Strait, T. E. 2013, ApJL, 765, L17

* Crossfield, I. J. M. 2013, A&A, 551, A99

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

* Dawson, R. I., & Fabrycky, D. C. 2010, ApJ, 722, 937

* de Kok, R. J., Brogi, M., Snellen, I. A. G., et al. 2013, A&A, 554, A82

* Delfosse, X., Bonfils, X., Forveille, T., et al. 2013a, A&A, 553, A8

* Delfosse, X., Donati, J.-F., Kouach, D., et al. 2013b, in SF2A-2013: Proc. Annual Meeting of the French Society of Astronomy and Astrophysics, ed. L. Cambresy et al., 497

* Delfosse, X., Forveille, T., Perrier, C., & Mayor, M. 1998, A&A, 331, 581

* Delfosse, X., Forveille, T., Ségransan, D., et al. 2000, A&A, 364, 217

* Demory, B.-O., Queloz, D., Alibert, Y., Gillen, E., & Gillon, M. 2016, ApJL, 825, L25

* Dittmann, J. A., Irwin, J. M., Charbonneau, D., et al. 2017, Natur, 544, 333

* Donati, J.-F., Catala, C., Landstreet, J. D., & Petit, P. 2006, in ASP Conf. Ser. 358, ed. R. Casini & B. W. Lites (San Francisco, CA: ASP), 362

* Hébrard, É. M., Donati, J.-F., Delfosse, X., et al. 2014, MNRAS, 443, 2599

* Hébrard, É. M., Donati, J.-F., Delfosse, X., et al. 2016, MNRAS, 461, 1465

* Donati, J.-F., Morin, J., Petit, P., et al. 2008, MNRAS, 390, 545

* Donati, J.-F., Yu, L., Moutou, C., et al. 2017, MNRAS, 465, 3343

* Dravins, D., Lindegren, L., & Nordlund, A. 1981, A&A, 96, 345

* Dressing, C. D., & Charbonneau, D. 2015, ApJ, 807, 45

* Dressing, C. D., Charbonneau, D., Dumusque, X., et al. 2015, ApJ, 800, 135

* Dumusque, X. 2016, A&A, 593, A5

* Dumusque, X., Boisse, I., & Santos, N. C. 2014, ApJ, 796, 132

* Dumusque, X., Borsa, F., Damasso, M., et al. 2017, A&A, 598, A133

* Dumusque, X., Udry, S., Lovis, C., Santos, N. C., & Monteiro, M. J. P. F. G. 2011, A&A, 525, A140

* Fabrycky, D. C., Ford, E. B., Steffen, J. H., et al. 2012, ApJ, 750, 114

* Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ, 790, 146

* Faria, J. P., Haywood, R. D., Brewer, B. J., et al. 2016, A&A, 588, A31

* Figueira, P., Adibekyan, V. Z., Oshagh, M., et al. 2016, A&A, 586, A101

* Figueira, P., Marmier, M., Boué, G., et al. 2012, A&A, 541, A139

* Fischer, D. A., Anglada-Escude, G., Arriagada, P., et al. 2016, PASP, 128, 066001

* Ford, E. B., & Gregory, P. C. 2007, in ASP Conf. Ser. 371, Statistical Challenges in Modern Astronomy IV, ed. G. J. Babu & E. D. Feigelson (San Francisco, CA: ASP), 189

* Ford, E. B., Seager, S., & Turner, E. L. 2001, Natur, 412, 885

* Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

* Forveille, T., Bonfils, X., Delfosse, X., et al. 2009, A&A, 493, 645

* Fujii, Y., Kawahara, H., Suto, Y., et al. 2010, ApJ, 715, 866

* Fujii, Y., Kawahara, H., Suto, Y., et al. 2011, ApJ, 738, 184

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

* Gaidos, E., Mann, A. W., Kraus, A. L., & Ireland, M. 2016, MNRAS, 457, 2877

* Gillon, M., Jehin, E., Delrez, L., et al. 2013, Protostars and Planets VI Posters

* Gillon, M., Jehin, E., Magain, P., et al. 2011, in European Physical Journal Web of Conf. 11, Detection and Dynamics of Transiting Exoplanets, ed. F. Bouchy, R. Díaz, & C. Moutou (Les Ulis: EDP Sciences), 06002

* Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Natur, 542, 456

* Gladman, B. 1993, Icar, 106, 247

* Gomes da Silva, J., Santos, N. C., Bonfils, X., et al. 2012, A&A, 541, A9

* Haywood, R. D., Collier Cameron, A., Queloz, D., et al. 2014, MNRAS, 443, 2517

* Huélamo, N., Figueira, P., Bonfils, X., et al. 2008, A&A, 489, L9

* Irwin, J. M., Berta-Thompson, Z. K., Charbonneau, D., et al. 2015, in 18th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, ed. G. T. van Belle & H. C. Harris (Cambridge: Cambridge Univ. Press), 767

* Kasper, M., Beuzit, J.-L., Verinaud, C., et al. 2010, Proc. SPIE, 7735, 77352E

* Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icar, 101, 108

* Kipping, D. M. 2013, MNRAS, 434, L51

* Kitchatinov, L. L., & Olemskoy, S. V. 2011, MNRAS, 411, 1059

* Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, ApJ, 765, 131

* Lafrenière, D., Doyon, R., Marois, C., et al. 2007a, ApJ, 670, 1367

* Lafrenière, D., Marois, C., Doyon, R., Nadeau, D., & Artigau, É 2007b, ApJ, 660, 770

* Lanotte, A. A., Gillon, M., Demory, B.-O., et al. 2014, A&A, 572, A73

* Lissauer, J. J., Ragozzine, D., Fabrycky, D. C., et al. 2011, ApJS, 197, 8

* Lopez, E. D., & Fortney, J. J. 2014, ApJ, 792, 1

* Lovis, C., Mayor, M., Bouchy, F., et al. 2005, A&A, 437, 1121

* Lovis, C., Snellen, I., Mouillet, D., et al. 2017, A&A, 599, A16

* Luger, R., Sestovic, M., Kruse, E., et al. 2017, NatAs, 1, 0129

* Luhman, K. L., Stauffer, J. R., Muench, A. A., et al. 2003, ApJ, 593, 1093
Macintosh, B., Troy, M., Doyon, R., et al. 2006, Proc. SPIE, 6272, 62720N

Mahadevan, S., Ramsey, L., Bender, C., et al. 2012, Proc. SPIE, 8446, 84461S

Mahmud, N. I., Crockett, C. J., Johns-Krull, C. M., et al. 2011, ApJ, 736, 123

Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556

Martín, E. L., Guenther, E., Zapatero Osorio, M. R., Bouy, H., & Wainscoat, R. 2006, ApJL, 644, L75

Mayor, M., Pepe, F., Queloz, D., et al. 2003, Msngr, 114, 20

Mayor, M., & Queloz, D. 1995, Natur, 378, 355

Meunier, N., Desort, M., & Lagrange, A.-M. 2010, A&A, 512, A39

Meunier, N., Lagrange, A.-M., Mbemba Kabuiku, L., et al. 2017, A&A, 597, A52

Montet, B. T., Morton, T. D., Foreman-Mackey, D., et al. 2015, ApJ, 809, 25

Morin, J., Donati, J.-F., Petit, P., et al. 2008, MNRAS, 390, 567

Moutou, C., Hébrard, E. M., Morin, J., et al. 2017, MNRAS, 472, 4563

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

Newton, E. R., Irwin, J., Charbonneau, D., et al. 2016, ApJ, 821, 93

Nutzman, P., & Charbonneau, D. 2008, PASP, 120, 317

O'Neal, D., Saar, S. H., Neff, J. E., & Cuntz, M. 2005, in ESA Special Publication 560, 13th Cambridge Workshop on Cool Stars, Stellar Systems and the Sun, ed. F. Favata, G. A. J. Hussain, & B. Battrick (Paris: ESA), 853

Pan, M., & Schlichting, H. E. 2018, ApJ, submitted, arXiv:1704.07836

Pont, F., Sing, D. K., Gibson, N. P., et al. 2013, MNRAS, 432, 2917

Prato, L., Huerta, M., Johns-Krull, C. M., et al. 2008, ApJL, 687, L103

Quillen, A. C., Morbidelli, A., & Moore, A. 2007, MNRAS, 380, 1642

Quirrenbach, A., Amado, P. J., Caballero, J. A., et al. 2014, Proc. SPIE, 9147, 91471F

Rajpaul, V., Aigrain, S., Osborne, M. A., Reece, S., & Roberts, S. 2015, MNRAS, 452, 2269

Rajpaul, V., Aigrain, S., & Roberts, S. 2016, MNRAS, 456, L6

Rasmussen, C. E., & Williams, C. K. I. 2005, Gaussian Processes for Machine Learning (Cambridge, MA: MIT Press)

Reiners, A., & Basri, G. 2007, ApJ, 656, 1121

Reiners, A., Bean, J. L., Huber, K. F., et al. 2010, ApJ, 710, 432

Reiners, A., Shulyak, D., Anglada-Escudé, G., et al. 2013, A&A, 552, A103

Robertson, P., Bender, C., Mahadevan, S., Roy, A., & Ramsey, L. W. 2016, ApJ, 832, 112

Rodríguez-López, C., Gizis, J. E., MacDonald, J., Amado, P. J., & Carosso, A. 2015, MNRAS, 446, 2613

Rogers, L. A. 2015, ApJ, 801, 41

Route, M. 2016, ApJL, 830, L27

Saar, S. H., & Donahue, R. A. 1997, ApJ, 485, 319

Scargle, J. D. 1982, ApJ, 263, 835

Schmidt, S. J., Kowalski, A. F., Hawley, S. L., et al. 2012, ApJ, 745, 14

Shulyak, D., Reiners, A., Engeln, A., et al. 2017, NatAs, 1, 0184

Shulyak, D., Reiners, A., Seemann, U., Kochukhov, O., & Piskunov, N. 2014, A&A, 563, A35

Snellen, I., de Kok, R., Birkby, J. L., et al. 2015, A&A, 576, A59

Snellen, I. A. G., de Kok, R. J., de Mooij, E. J. W., & Albrecht, S. 2010, Natur, 465, 1049

Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJL, 755, L28

Stewart, G. R., & Ida, S. 2000, Icar, 143, 28

Suárez Mascareño, A., Rebolo, R., & González Hernández, J. I. 2016, A&A, 595, A12

Tamura, M., Suto, H., Nishikawa, J., et al. 2012, Proc. SPIE, 8446, 84461T

Trauger, J., Gordon, B., Krist, J., & Moody, D. 2015, Proc. SPIE, 9605, 96050N

Valencia, D., Guillot, T., Parmentier, V., & Freedman, R. S. 2013, ApJ, 775, 10

Van Eylen, V., & Albrecht, S. 2015, ApJ, 808, 126

Vanderburg, A., Plavchan, P., Johnson, J. A., et al. 2016, MNRAS, 459, 3565

Vogt, S. S., Allen, S. L., Bigelow, B. C., et al. 1994, Proc. SPIE, 2198, 362

Wargelin, B. J., Saar, S. H., Pojmański, G., Drake, J. J., & Kashyap, V. L. 2017, MNRAS, 464, 3281

Weiss, L. M., & Marcy, G. W. 2014, ApJL, 783, L6

West, A. A., Weisenburger, K. L., Irwin, J., et al. 2015, ApJ, 812, 3

Yu, L., Donati, J.-F., Hébrard, E. M., et al. 2017, MNRAS, 467, 1342