Searching For Transiting Planets Around Halo Stars. ii. Constraining the Occurrence Rate of Hot Jupiters

The evolution of the Milky Way and planet formation are connected (Forgan et al. 2017). Metallicity links these seemingly unrelated topics by playing an important role in the formation and evolution of stars and therefore affects the formation of giant planets (Johnson & Apps 2009; Choi & Nagamine 2009). Host-star metallicity is thought to reflect the metallicity of the protoplanetary disk from which planets form (Wyatt et al. 2007; Haworth et al. 2016). Given that the first generation of stars was metal poor, the first planets must also be lacking in metals. By probing planet occurrence in the metal-poor regime, we can gain insight into the first generation of giant planet formation and determine the metallicity at which no planet can form.

A metal-poor environment poses many challenges for planet formation. At low metallicities, the protoplanetary disk lifetime is significantly shorter which decreases the likelihood of planet formation (Kornet et al. 2005; Ercolano & Clarke 2010; Yasui et al. 2010). The vast majority of the gas within the disk is hydrogen and helium, and heavy elements compose only a fraction ∼10⁻⁵ of the total mass in the metal-poor regime (Johnson & Li 2012). A shortened disk lifetime can be detrimental to planet formation, especially at short orbital distances where planet formation is thought to be dominated by core accretion (Miller & Fortney 2011; Bailey & Batygin 2018). Although we know that planet formation is increasingly improbable with decreasing metallicity, the metallicity for which no planet can form is still uncertain.

Since the discovery of four Jovian planets orbiting metal-rich stars in 1997 by Gonzalez et al. (2001), the correlation between planet occurrence and host-star metallicity has been a point of focus within the exoplanet community. Following this discovery, various studies also noted that Jovian planets preferentially orbit metal-rich stars (Santos et al. 2004; Fischer & Valenti 2005; Mortier et al. 2013). The occurrence of Jovian planets has been shown to increase exponentially with the increase of metallicity for solar and supersolar metallicities (Udry & Santos 2007; Johnson et al. 2010; Mortier et al. 2013); however, the hot-Jupiter (HJ) occurrence trend within the metal-poor regime is still uncertain. Using radial velocities, Sozzetti et al. (2009) and Mortier et al. (2012) probed this regime to determine Jovian planet occurrence and provide an upper limit for HJ occurrence in the metal-poor regime. Since the study conducted by Mortier et al. (2012), there has been little investigation to further constrain HJ planet occurrence at low metallicities.

With the advent of the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015), it is the opportune time to revisit HJ occurrence in the metal-poor regime as TESS provides a crucial advantage compared to Kepler to complete this work. The major advantage of TESS is its larger field of view, which allows for 30-minute cadence full-frame images (FFI). These FFIs provide observations of a wide range of stellar types. TESS also covers ∼85% of the sky, which is a significant improvement to the ∼0.25% monitored by Kepler. The capabilities of TESS have enabled a more robust determination of the HJ occurrence within the metal-poor regime that was not possible until now.

In this paper, we focus on further constraining HJ planet occurrence in the metal-poor regime (−2.0 ≤ [Fe/H] ≤ −0.6) using a population of halo stars, as they provide important insights into planet-formation processes in the first generation of stars within the Milky Way. Our approach to deriving HJ occurrence rates is largely inspired by the method first outlined by Dressing & Charbonneau (2015). We perform a full transit search of our halo star sample, using both the Box Least Squares algorithm (Kovács et al. 2002) and Transit Least Squares algorithm (Hippke & Heller 2019) to search the light curves. Using an injection-recovery scheme to calculate our detection efficiency, we use our results to constrain the occurrence rate of HJs orbiting halo stars.

The outline of our paper is as follows. Section 2 provides a description of our stellar sample. In Section 3, we discuss the method from which we acquire TESS data. We explain our planet detection pipeline in Section 4. Section 5 contains a description of our planet injection pipeline. In Section 6, we assess the completeness of our pipeline. Section 7 is dedicated to a statistical analysis of the data and estimation of the planet occurrence rate. We compare our results to previous studies in Section 8 and their implications for the frequency of HJs in the metal-poor regime before concluding in Section 9.

We selected ∼16,940 halo stars, which are typically subdwarfs, to search for HJs using TESS data. Details on the sample selection and validation can be found in a companion paper (Kolecki et al. 2021); however, we briefly describe the sample selection here. The sample is selected using the transverse kinematic information, which was constructed following Koppelman et al. (2018). Using radial velocities and proper motions from Gaia DR2 (Gaia Collaboration et al. 2018), we select stars within 1 kpc and with 2D velocities greater than 210 km s⁻¹ with respect to the local standard of rest. The sample was chosen using heuristic tangential motion cuts based on simulations of the Milky Way from Galaxia (Sharma et al. 2011), which is outlined in Kolecki et al. (2021). By conducting a population study of ∼120 stars in our sample with APOGEE spectra (Kolecki et al. 2021), we confirmed that 70% of the selected targets are genuine halo stars. Table 1 list the ranges of various parameters of the sample, which were obtained via the TESS Input Catalog (Stassun et al. 2018). An H–R diagram of the sample is also provided in Figure 1. The discontinuity in the distribution is due to the in-homogeneous treatment to low-mass stars and other stars in the TESS Input Catalog (see Appendix E.1.1 in Stassun et al. 2018).

Zoom In Zoom Out Reset image size

From the initial sample, we obtained 30-minute cadence data for ∼65% (11,125 out of 16,940) from TESS through sector 23 (Ricker et al. 2015). Given that only a small fraction (∼1%) of the sample had TESS 2 minute cadence data, we chose to take advantage of a majority of the sample using the 30-minute cadence data. These light curves were obtained from TESS FFI with Eleanor (Feinstein et al. 2019). The FFI light curves are then used for planet detection and injection-recovery tests.

Eleanor (Feinstein et al. 2019), software that was designed for downloading and analyzing TESS FFI data, is employed to acquire the 30-minute cadence TESS FFIs for our sample. Using the TESS magnitude of each target, we determine the optimal aperture size, which is then input into Eleanor. We download the four different data files: PCA, Correlated, Raw, and PSF. From there, we determine the noise level of each light curve. We normalize the flux data by dividing the light curve by the median flux. The noise of each file is determined by calculating the RMS noise of the light curve. The data file with the least amount of noise is saved.

To determine the occurrence rates and completeness of our search, we must search within our sample for planet candidates. We developed an automated pipeline that encompasses all the data processing to include detrending, noise removal, and transit search. Given that our target stars are relatively faint, we must minimize the noise where possible. From acquiring the FFI data, we choose the light curve with the least amount of noise. We also apply the most robust detrending method available within the detrending program used.

We begin by detrending the light curves. The data is initially detrended using Wotan (Hippke et al. 2019) with the Tukey's biweight method, which is indicated to be the most robust detrending method within their work (Hippke et al. 2019). Since the TESS data has sharp peaks at the beginning and ends of an observation which we consider to be "data gaps," we perform sigma clipping to remove any systematics that Wotan was unable to remove. We determine the outliers at the edges of the light curves and at the edges of the data gaps. We only consider 5% of the data points on the edges when conducting the sigma clipping as they most negatively impact the transit search algorithms. We then calculate the median of the complete light curve and any points that are over 3σ within the range considered to be the edges of the data are removed. From here our pipeline splits into two separate transit search algorithms: transit least squares (TLS; Hippke & Heller 2019) and box least squares algorithm (BLS; Kovács et al. 2002).

TLS, an open source transit search python package, is noted to be more efficient at conducting transit searches for small planets (Hippke & Heller 2019). To expand upon their work, we chose to test whether it also performs better than BLS for giant planets around halo stars. To test how well TLS performs, we used the default transit template. For each star, we also input the limb-darkening coefficients, mass, and radius to optimize for our sample.

For our BLS algorithm, we used the reference implementation of BLS provided by astropy 4.0.1 (Astropy Collaboration et al. 2013). With this implementation the false alarm probability (FAP) and signal-to-noise ratio (S/N) are not provided; however, they are required within the vetting process. In Section 4.1, we expand upon how we calculate the FAP and S/N given that these values are not calculated by the BLS program that we employed.

4.1. BLS

Given our choice to use two transit search algorithms, we make every effort to conduct a fair comparison. Following Hippke & Heller (2019), we empirically determined the log likelihood for a FAP value of 0.01%. Similar to TLS, we generated 10,000 synthetic light curves with only white noise with a 1 hr Combined Differential Photometric Precision of 290 ppm. We choose this noise level to replicate the typical noise level for TESS for stars with a TESS magnitude of 12 (Ricker et al. 2015), which is the mean TESS magnitude for our sample. It is important to note that TESS light curves have correlated noise resulting from instrumental trends; however, we use white noise to ensure an equal comparison between TLS and BLS. Once the white-noise-only light curves have been evaluated by BLS, we impose a FAP value of 0.01% determining the corresponding log-likelihood threshold for BLS.

Within BLS, the S/N is not automatically determined; therefore, we calculate a S/N for each fit as

$$\begin{eqnarray}&&\mathrm{SNR}=\displaystyle \frac{d}{\sigma }\sqrt{{n}_{\mathrm{transit}}}\end{eqnarray} \tag{ 1 }$$

where d is the difference between the median of the data points in the transit and the flux level outside the transit, n_transit is the number of points within transit, and σ is the standard deviation of the detrended light curve.

4.2. Edi-vetter Unplugged

4.3. Vetting

4.4. Follow-up Observations & Results

To accurately measure the planet occurrence rate based on the results of our planet search, we need to determine the probability of detecting a planet using our pipeline. We measure the detection efficiency of our planet detection pipeline by injecting transiting planet signals into the TESS light curves and running the light curves through our detection pipeline where they will be detrended and vetted. We generate synthetic transit signals for each light curve with orbital parameters drawn from a uniform distribution using pylightcurve (Tsiaras et al. 2016). Each synthetic planet signal that is injected into a light curve has a randomly determined planetary radius, orbital period, midtime, and inclination (Rp, P, t₀, i) taken from the following uniform distributions:

$$\begin{eqnarray*}&&\displaystyle \frac{P}{\mathrm{day}}\sim U(0.5,20)\end{eqnarray*}$$

$$\begin{eqnarray*}&&\displaystyle \frac{{R}_{p}}{{R}_{J}}\sim U(0.08,2)\end{eqnarray*}$$

$$\begin{eqnarray*}&&{t}_{0}\sim U(0,P)\end{eqnarray*}$$

$$\begin{eqnarray*}&&i\sim U({i}_{\min },{i}_{\max })\end{eqnarray*}$$

where i_max and i_min are calculated for each star as the maximum and minimum inclination for a transit using the semimajor axis and the radius of each star. The eccentricity is set to 0, given that HJs are subject to strong orbital circularization and are expected to typically have low eccentricities (e.g., Alvarado-Montes & García-Carmona 2019). Using stellar parameters, the limb-darkening coefficients are determined by interpolating the Claret (2017) limb-darkening tables. Specifically, they are estimated with the effective temperature, surface gravity, and the median metallicity of our sample (i.e., [Fe/H] = −1.27), which obtained from the TESS Input Catalog (Stassun et al. 2018). After each light curve is modified to contain a synthetic planet signal, it is then sent through our planet detection pipeline.

For a synthetic planet signal to be considered a detection, we used a nearly identical vetting process to the detection pipeline. We require the S/N to be greater than 5σ. The modified light curves must be determined by Edivetter Unplugged not to be a false positive. We also require the FAP value to be ≤0.0001. After the modified light curve undergoes this vetting process, the remaining synthetic planet signals are considered a detection if the detected period of the planet matches within 1% of the generated period. Period aliases of the injected signal are not treated as detections. It is important to note that, unlike real data, we do not conduct visual inspections as it would be highly impractical given the large number of injections.

5.1. Assessing our Pipeline

To assess the detection efficiency of our pipeline, we conducted extensive injection-recovery tests. For each star, we generated 200 synthetic transit signals per light curve. In total, we injected 2,225,000 transiting planets into the light curves of the 11,125 stars for both TLS and BLS. Our pipeline successfully recovered ∼80% of injected planet signals with an S/N ≥ 5 and FAP ≤ 0.0001 within the range of 0.8–2 R_J and periods 0.5–10 days with both TLS and BLS. We used both TLS and BLS and compare their results; however, we will focus on TLS within this section and elaborate on the difference between TLS and BLS within Section 5.2.

We calculate the detection efficiency by dividing the number of the synthetic planet signals detected by the total planet signals injected. We take this fraction as the detection efficiency, which is binned as a function of planet period and radius. As shown in Figure 2, we found that our pipeline is sensitive to injected planets with radii larger than 0.8 R_J . Most of those planet signals were detected with nearly 75% efficiency out to the orbital period of ∼10 days. Our pipeline displayed a significant decrease in the detection efficiency of planets with radii <0.5 R_J with the recovery fraction ∼19% and 0.5–2.0 R_J planets with periods longer than 10 days with the recovery fraction ∼12%; however, the decrease in detection efficiency for smaller planets and longer periods does not affect our ability to constrain the occurrence rate for HJs. While we focus HJs in this work, we note that it is possible to detect close-in Neptune-sized planets given our detection efficiency.

Zoom In Zoom Out Reset image size

5.2. Comparison of TLS versus BLS

We define the search completeness to be the probability of finding a planet with a random orbital alignment within our sample using our detection pipeline. Therefore, the overall search completeness depends both on the detectability, P_det(R_p , P), of a particular planet and the likelihood that it will be observed to transit. We determine the detection efficiency of our pipeline by conducting extensive injection-recovery tests using our planet injection pipeline (see Section 5).

Once the detection efficiency was determined, we use the geometric transit probability to account for the likelihood of a planet being oriented so that a transit can be observed. For the geometric transit probability we assume an eccentricity of zero and use:

$$\begin{eqnarray}&&{P}_{t}({R}_{p},P)=\displaystyle \frac{{R}_{p}+{R}_{\star }}{a}\end{eqnarray} \tag{ 2 }$$

where R_⋆ is the star radius and a is the orbital semimajor axis. a is determined from Kepler's third law assuming a stellar mass of 1.01 M_⊙ and radius of 0.711 R_⊙, which are the median values for the halo stars within our sample.

For a given orbital period and planet radius, we computed the corresponding semimajor axis for a planet orbiting each of the stars in our sample. With the transit probability for each radius and orbital period, we calculated the search completeness by weighting the detection efficiency by the geometric transit probability for each bin. In Figure 3, we show the search completeness.

Zoom In Zoom Out Reset image size

To determine the upper limit of planet occurrence given a null detection, we follow the procedure originally described in the appendix of Burgasser et al. (2003) and in McCarthy & Zuckerman (2004) using the binomial distribution. The probability P(f_p ) of finding d detections in a sample of size N can be calculated as a function of the true planet frequency f_p :

$$\begin{eqnarray}&&P(d,N,{f}_{p})={f}_{p}^{d}{\left(1-{f}_{p}\right)}^{N-d}\displaystyle \frac{N!}{(N-d)!d!}.\end{eqnarray} \tag{ 3 }$$

In the case of a null detection d = 0, the planet frequency is constrained to the interval ($0,{f}_{p,\max }$) for a confidence interval, C by

$$\begin{eqnarray}&&{\int }_{0}^{{f}_{p,\max }}(N+1)\times P(0,N,{f}_{p}){{df}}_{p}=C\end{eqnarray} \tag{ 4 }$$

where (N + 1) is a normalization factor. Using this equation, the maximum planet frequency can be solved for

$$\begin{eqnarray}&&{f}_{p,max}=1-{\left(1-C\right)}^{{\textstyle \tfrac{1}{N+1}}}.\end{eqnarray} \tag{ 5 }$$

Therefore, we can calculate an upper limit of f_p at a given confidence level by setting C to the corresponding probability for the confidence level of interest.

However, not all planets are detectable using the transit method and N must be replaced by the effective sample size N_eff. To calculate the true upper limit of planet occurrence, we must account for the geometric transit probability and the detection efficiency of our pipeline. Therefore, we define our effective sample size as:

$$\begin{eqnarray}&&{N}_{\mathrm{eff}}({R}_{p},P)=N\times {P}_{t}({R}_{p},P)\times {P}_{\det }({R}_{p},P)\end{eqnarray} \tag{ 6 }$$

where N is the total sample size, P_t (R_p , P) is the geometric transit probability, ${P}_{\det }({R}_{p},P)$ is the detection efficiency (see Section 6), P the orbital period, and R_p the planet radius.

Using Equations (5) and (6), we calculate the upper limit of HJ occurrence by measuring the range in f_p that covers 68% of the integrated probability function. This is equivalent to the 1σ confidence level. We conservatively choose this confidence level as to not overestimate the upper limit of HJ occurrence.

In Figure 4, we show our results using the detection efficiency determined from BLS. For the upper limit of HJ occurrence, we consider a radius of 0.8–2 R_J and periods 0.5–10 days. We choose a radius of 2 R_J given that the majority of the confirmed planets to date are less than 2 R_J (Akeson et al. 2013). Within this region, we find the upper limit of occurrence range between 0.04% ≤ f_p ≤ 0.36% within the parameter space considered with the mean upper limit of occurrence being f_p < 0.18%.

Zoom In Zoom Out Reset image size

The occurrence of HJs has been studied by multiple surveys. For HJs orbiting solar-type stars, radial velocity (RV) surveys have found the occurrence to be between 0.9% and 1.2% (Marcy et al. 2005; Mayor et al. 2011; Wright et al. 2012; Petigura et al. 2018), which is in contrast to the significantly lower occurrence rates determined with transit surveys. The HJ occurrence observed using TESS and Kepler range between 0.4% and 0.6% (Howard et al. 2012; Fressin et al. 2013; Mulders et al. 2015; Santerne et al. 2016; Petigura et al. 2018; Zhou et al. 2019). Although TESS and Kepler observe different targets, the occurrence rates are consistent.

The discrepancy of occurrence rates between RV and transit surveys has been addressed by several works. Many studies have speculated that the difference may be a result of selection bias. Guo et al. (2017) investigated the differences in metallicity between Kepler and some radial velocity surveys; however, they found the difference in metallicity to be insufficient to explain the HJ occurrence discrepancy. Another possible explanation for the discrepancy may be attributed to the contamination of multiple stellar systems and evolved stars in transit samples (Wang et al. 2015). Most recently, Moe & Kratter (2019) concluded that selection effects are the most probable cause of the HJ occurrence discrepancy. They quantitatively demonstrate that RV surveys for giant planets can increase the occurrence rate by a factor of two by systematically removing spectroscopic binaries from their samples.

The first attempt to determine HJ occurrence within the metal-poor regime was Sozzetti et al. (2009). They used RV data with most stars having 4–10 measurements from HIRES on the Keck i telescope with a sample of 160 metal-poor stars (see Table 2). Consequently, strong constraints were derived for short-period gas giants where they had 95% completeness for periods ≤3 yr with a minimum companion mass of ∼0.75 M_J , which corresponds to a HJ radius of approximately 1.25 R_J (Chen & Kipping 2017). Given that the majority of the range is for cool Jupiters, we consider radius as a step function using the Chen & Kipping (2017) for periods <10 days and Thorngren et al. (2019) for periods >10 days. Within the HJ range, we expect the planets to be inflated thereby increasing the radius. Their work resulted in a 1σ upper limit of HJ occurrence of 0.67%.

Table 2. Hot-Jupiter Occurrence Comparison

Author	N	fp	[Fe/H]	P	Rp	Method
Petigura et al. (2018)	1883 a	${0.57}_{-0.12}^{+0.14} \%$	(−0.3 ≤ [Fe/H] ≤ 0.3)	1 ≤ P ≤ 10 days	0.71-2.14 RJ	Transit
Wright et al. (2012)	836	1.2 ± 0.38%	(−0.4 ≤ [Fe/H] ≤ 0.3)	≤10 days	>0.51 RJ	RV
Zhou et al. (2019)	2401 a	0.4 ± 0.10%	(−0.5 ≤ [Fe/H] ≤ 0.4)	0.9 ≤ P ≤ 10 days	0.8–2.5 RJ	Transit
Mayor et al. (2011)	822	0.89 ± 0.36%	(−0.5 ≤ [Fe/H] ≤ 0.3)	≤11 days	>0.75 RJ	RV
Masuda & Winn (2017)	1737 a	${0.43}_{-0.06}^{+0.07} \%$	(−0.6 ≤ [Fe/H] ≤ 0.2)	0.5 ≤ P ≤ 10 days	0.8-2 RJ	Transit
Mortier et al. (2012)	114	1.00% b	(−1.4 ≤ [Fe/H] ≤ −0.4)	≤10 days	>1.04 RJ	RV
Sozzetti et al. (2009)	160	0.67% b	(−1.6 ≤ [Fe/H] ≤ −0.6)	<3 yrs	>1.25 RJ	RV
This Work	576 a	0.18% b	(−2.0 ≤ [Fe/H] ≤ −0.6)	0.5 ≤ P ≤ 10 days	0.8-2 RJ	Transit

Notes. This table is arranged in terms of the minimum metallicity. We use the 5 to 95% quantiles for each study. For all radial velocity surveys, we use the Chen & Kipping (2017) relation to convert from mass to radius.

^a We estimate the median effective sample size using our detection efficiency. ^b This is the maximum 1σ upper limit of occurrence.

Download table as:  ASCIITypeset image

The most recent work to directly consider this regime for HJs is Mortier et al. (2012). Similar to Sozzetti et al. (2009), RV measurements were taken using HIRES and HARPS with a sample of 114 stars. To avoid biases, the sample was restricted to targets that had six or more measurements. By imposing this requirement, only 50 of the target from Sozzetti et al. (2009) were used in their work. For a completeness of 80% and period < 10 days, they found a minimum companion mass of ∼0.314 M_J , which corresponds to a radius of 1.04 R_J (Chen & Kipping 2017). Once again, we use the mass–radius relationship from Chen & Kipping (2017) because we expect the HJs to be inflated. Given the smaller sample size, they found the 1σ upper limit of HJ occurrence to be 1.00%.

Compared to the previous works using radial velocities, we place the most stringent upper limit on HJ occurrence. Using TESS FFI data, we allow for a larger effective sample size that ranges from N_eff = 313–2582 per bin with a median of N_eff = 576 (see Equation (6)). Our minimum effect sample size is ∼2–3 times larger than the sample sizes of Mortier et al. (2012) and Sozzetti et al. (2009; see Table 2). Another difference to note is that each study has varying ranges of radius and period (see Table 2); however, this does not significantly affect the comparison. For example, our upper limit of occurrence within the range Mortier et al. (2012) considered is between 0.04% ≤ f_p ≤ 0.34% with the mean upper limit of 0.17%.

From the results of this work, a metallicity limit can be established below which no HJs can form. Within Mortier et al. (2012), the authors suggested that giant planets cannot form below [Fe/H] = −0.5. Similarly, we propose that the metallicity limit for HJs would be between −0.7 and −0.6, which is the upper limit of the metallicity range for this work. This value is consistent with observations and previous works (Mordasini et al. 2012; Mortier et al. 2012). To date, no HJ has been discovered below a metallicity of −0.6 (Akeson et al. 2013). Currently, the most metal-poor star to host a HJ has a stellar metallicity of −0.6 (Hellier et al. 2014).

Although there is a strong correlation between metallicity, it is important to note that we are referring to the iron abundance of the star. However, there is evidence that indicates the overall metallicity of iron-poor HJ hosts is higher than the iron abundance (Haywood 2008; Adibekyan et al. 2012a, 2012b). Stars of the thin and thick disks have an increase in α element as a function of metallicity (Haywood 2008; Adibekyan et al. 2012b). As a result, almost all giant planet hosts in the metal-poor regime have high [α/Fe] values, such as Mg and Si (Adibekyan et al. 2012b). This indicates that overall metallicity is, in fact, higher than the iron abundance and that the α element aids in the formation of HJs in the metal-poor regime.

In this paper, we presented an updated upper limit of HJ occurrence within the metal-poor regime using TESS data. We developed our own planet detection pipeline to search for transiting planets within our light curves. We then characterized the completeness of our pipeline by injecting simulated transiting planets into the light curves and attempted to recover them. Our search of the light curves of 11,125 halo stars revealed no planet candidates. We then calculated the upper limit of occurrence for HJs around halo stars by using the null detection hypotheses accounting for the geometric transit probability and the detection efficiency of our pipeline. The primary conclusions of this work are:

• 1.  
  In line with previous works, HJs are rare or nonexistent in the metal-poor regime (−2.0 ≤ [Fe/H] ≤ −0.6).
• 2.  
  We further constrain the upper limit HJ occurrence within the metal-poor regime using TESS data. We find the upper limit of HJ occurrence range from 0.04% ≤ f_p ≤ 0.36% with the mean occurrence being f_p < 0.18% for radii 0.8–2 R_J and periods 0.5–10 days at the 1σ confidence level, which places the most stringent constraint on HJ occurrence.
• 3.  
  TLS is comparable to the BLS when searching for HJs around metal-poor stars. It performs slightly better than BLS overall.

Further constraining the occurrence of HJs would require significantly larger samples; however, any planet discoveries would be of great interest within the metal-poor regime. They would expand our knowledge of the environment in which the first generation of planets formed.

Moving forward, our future work will build upon this work. Using the framework that was developed conducting this study, we plan to determine an accurate functional form of planet occurrence versus metallicity (Fischer & Valenti 2005; Udry & Santos 2007; Johnson et al. 2010). We will consider a larger sample size and utilize the pipelines created for this work. This is a starting point in constructing a galactic planet formation model, which relates the metal-enrichment of the Milky Way to planet formation.

We thank David Latham for the effort in organizing the TESS Follow-up Observing Program Working Group. We acknowledge the MEarth team for providing the follow-up observation. We thank Samuel Quinn for the follow-up observations as well. T.S.L. is supported by NASA through Hubble Fellowship grant HST-HF2-51439.001, awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. J.C.Z. is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2001869. We also thank the anonymous referee and Johanna Teske who both provided suggestions that improved the quality and credibility of the manuscript.

Software: eleanor (Feinstein et al. 2019), pylightcurve (Tsiaras et al. 2016), transit least squares algorithm (Hippke & Heller 2019), and astropy (Astropy Collaboration et al. 2013).