Occurrence Rate of Hot Jupiters Around Early-type M Dwarfs Based on Transiting Exoplanet Survey Satellite Data

In order to obtain a high S/N transit detection and better understand the architecture of each system, we make use of all available QLP light curves of our stellar sample from both TESS primary and extended missions. We retrieve the light curves of each target from the Mikulski Archive for Space Telescopes (MAST ²⁷ ) via astroquery (Ginsburg et al. 2019). To improve the precision of the light curves, we ignore entries where the quality flag is assigned nonzero, which indicates anomalies in the data or images (Huang et al. 2020a). Despite this, the raw light curves of most stars still have a few data points with abnormally high flux values. We thus calculate the 99.5th percentile of each light curve and exclude 0.5% of the points with the highest flux for all stars in our sample, which might be related to instrumental or systematic noise. We show the TESS baseline length distribution of our final stellar sample in Figure 3.

Zoom In Zoom Out Reset image size

Second, we perform a uniform detrending by fitting basis spline models for each light curve. During the construction of the spline models, we use a running 3σ-clipped median filter. We divide the full light curve into several bins with a binning size of 0.3 day. Within each bin, we calculate the median flux after removing the 3σ outliers. Next, we interpolate the 0.3 day binned full spline that we obtained onto the full observation time stamps with a cubic interpolator. We finally produce the detrended light curve by dividing the original light curve by this interpolated spline function. We use these detrended light curves for candidate search.

Our planet detection pipeline is mainly based on the box least square (BLS; Kovács et al. 2002) algorithm. ²⁸ Following the methodology described in Dressing & Charbonneau (2015), we first perform a low-resolution transit search. We explore 1000 uniformly spaced period grids between 0.8 and 10 days. To determine the best duration searching grid, we randomly select 10,000 stars from our sample, generating arbitrary physical parameters with period P between 0.8 and 10 days, impact parameter b below 0.9, and planet size between 7 R_⊕ and 2 R_J, and compute the transit duration assuming a circular orbit (Seager & Mallén-Ornelas 2003). We find that more than 99% of the durations are located between 0.05 and 0.17 day. Consequently, to be conservative, we conduct our search for 10 uniformly spaced transit durations between 0.02 and 0.2 day, where all of our simulated duration values are located.

We first compute a BLS periodogram for each detrended light curve. To ensure a relatively clean sample without too many false positives, we require the selected candidates to have a BLS reported maximum S/N (S/N_transit) greater than 10. As we use a relatively sparse spline model to detrend the light curves, short-timescale sinusoidal-like stellar variations may be left in the data and cause false alarms. However, we do not expect a strong periodic brightening effect with a similar amplitude comparable with the dimming signal for real transit events. Therefore, we conduct an antitransit search by constructing another BLS periodogram for the flipped light curve (Wang et al. 2014) to identify stellar variability. Similarly, we record the BLS maximum S/N of the antitransit (S/N_antitransit). We define the S/N ratio (S/N_transit/S/N_antitransit) as an indicator that reflects the robustness of a real transit detection, and we only keep candidates with an S/N ratio greater than 1.5.

If a candidate passes the above low-resolution transit search, we next refine the period P and mid-transit timing T₀ by performing high-resolution BLS runs to examine alias signals. Starting with the period P_raw found in the previous step, we calculate all aliasing periods P_alias as P_raw/N and P_raw × N, where N is a positive integer, and save period values between 0.4 and 12 days. During the high-resolution runs, we focus on the period range of [P_alias − 0.1, P_alias + 0.1] with 1000 intervals and the same transit duration grid as in the previous search, and loop the BLS search for all aliasing periods. We regard the P_BLS–T_0,BLS–duration set with the highest BLS S/N as the final transit ephemeris.

Eclipsing binary systems generally have a significant secondary dip, manifesting as a depth difference between odd and even transits. On the contrary, we expect identical odd/even depths for real planetary signatures. To reject such astrophysical false positives, we compare the odd and even depths (δ_odd and δ_even) reported by the BLS algorithm. We calculate the odd/even depth difference Δ = ∣δ_odd − δ_even∣ and require all candidates to have δ_odd/Δ and δ_even/Δ greater than 3. This conservative threshold is set based on a test on several selected binary systems. A more careful investigation of the odd/even difference is performed in Section 4.3.

Our detection pipeline alerts a total of 437 events in the end. For each candidate, a diagnostic plot is generated as in Figure 4, which includes the raw QLP photometry, spline model detrended data, BLS periodogram of the low-resolution search as well as the phase-folded light curve to the transit ephemeris found in the high-resolution search. We note that our detection pipeline only examines the transit-like signal with the highest S/N, and it may miss giant planets around young M dwarfs with strong short-timescale variations. Such incompleteness will be characterized by the injection-recovery test. Additionally, in this work we do not deal with multiplanet cases, which have not been detected so far, as M-dwarf systems with close-in hot Jupiters and additional planets are rare.

Zoom In Zoom Out Reset image size

Except for the 437 candidates studied in this work, there are also known TOIs that were missed by our detection pipeline. We match the catalog of 60,382 stars without detection with a list of known TOIs, and we found a total of 67 TOIs that were not alerted. We present the full catalog of missed TOIs and report our search results in Table 5. Figure 5 shows their period–radius distribution. Most of these candidates were found in the 2 minutes cadence light curves extracted by the TESS Science Processing Operations Center (SPOC; Jenkins et al. 2016). Given their small companion size around 1–4 R_⊕ and short duration time, the transit depth will be diluted in long-cadence FFI data. Thus, the BLS S/Ns of these small planet candidates are generally low. Among all of these missed candidates, we note that there is one target (TIC 168751223/TOI-2331) that is within the parameter space we have searched in this work. It was not alerted by our detection pipeline because the S/N of our BLS search does not match our minimum threshold (S/N = 10). Further ground-based follow-up observations ruled out this candidate as by confirming that the eclipse signal is from a nearby binary star system on TIC 168751224 (ΔT = 3.1 mag) at 7'' away.

Zoom In Zoom Out Reset image size

The QLP produces light curves of each source using several apertures and identifies an optimal light curve based on the target brightness. Generally, the size of this optimal aperture is larger than two TESS pixels for M dwarfs with 10.5 ≤ T_mag ≤ 13.5 so the light from nearby eclipsing binaries within $1^{\prime}$ may pollute the aperture and cause transit-like signals on the target light curve. We reject such scenarios using the difference image technique (Bryson et al. 2013). We perform a pixel-level centroid offset analysis in the difference image of each candidate with TESS-plots ²⁹ (Kunimoto et al. 2022). TESS-plots downloads 20 × 20 pixel cutout of TESS FFIs obtained in a certain sector, generates a difference image based on the flux of in- and out-of-transit images and calculates the S/N of each pixel. The light centroid of the difference image should be located around the source position in the direct image (i.e., FFI) if the signal is on target. Otherwise, the nearby eclipsing binary scenario is favored when a large centroid shift happens. We compute a S/N-weighted light centroid (x_c, y_c) in a 7 × 7 pixel difference image centered on target for every candidate using:

$$\begin{eqnarray}&&{x}_{{\rm{c}}}=\displaystyle \frac{{\sum }_{i=1}^{7}{\sum }_{j=1}^{7}{\rm{S}}/{{\rm{N}}}_{{x}_{i},{y}_{j}}^{2}\times {x}_{i}}{{\sum }_{i=1}^{7}{\sum }_{j=1}^{7}{\rm{S}}/{{\rm{N}}}_{{x}_{i},{y}_{j}}^{2}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{y}_{{\rm{c}}}=\displaystyle \frac{{\sum }_{i=1}^{7}{\sum }_{j=1}^{7}{\rm{S}}/{{\rm{N}}}_{{x}_{i},{y}_{j}}^{2}\times {y}_{j}}{{\sum }_{i=1}^{7}{\sum }_{j=1}^{7}{\rm{S}}/{{\rm{N}}}_{{x}_{i},{y}_{j}}^{2}},\end{eqnarray} \tag{ 2 }$$

where x_i and y_j are the pixel indices, ${\rm{S}}/{{\rm{N}}}_{{x}_{i},{y}_{j}}$ represents the signal-to-noise of pixel (x_i , y_j ) in the zoomed difference image.

We compare the target position on the difference image with the measured light centroid in the FFI and calculate the centroid shift (d_c). We generate difference images from different sectors for each star and accept the cases where the light centroid shift of the highest S/N difference image is smaller than 1 pixel. An example, TIC 14081980 with d_c = 1.2 pixels, which we excluded in this step, is shown in the left panel of Figure 7. In total, 238 candidates survive after this step. For confirmed planets detected by TESS (Table 1), all of their centroid shifts are <0.5 pixels. The 1 pixel (∼21'') centroid offset is a conservative threshold to rule out signals from nearby binary or planetary systems, which is also previously used in alerting TOIs by TESS teams (Guerrero et al. 2021; Kunimoto et al. 2022). In the SPOC validation reports (Twicken et al. 2018), a 25 error term is added in quadrature to the propagated uncertainty in the difference image centroid offsets. The 3σ centroid offset level for single-sector observations is roughly 75 (0.35 TESS pixel) for a majority of target stars. This is much smaller than the 1 pixel choice here. The SPOC centroid shifts of confirmed or known planets alerted with transit depths between 5000 and 25,000 ppm are smaller than 1σ (<3''). Excluding candidates with centroid shifts larger than 1 pixel (∼21'') thus give a completeness much higher than 99.7% so the false-negative rate of this step is negligible for our study. Although Twicken et al. (2018) pointed out that the light centroid determination is sometimes unreliable for targets within crowded fields as nearby bright stars will have an effect, we note that we excluded M stars with dilution A_D > 0.3 (see Section 2). Therefore, the targets in our sample are relatively isolated and have, in principle, precise light centroid measurements.

Zoom In Zoom Out Reset image size

Next, we further exclude 44 events with centroid shifting to the same nearby star in an adjacent pixel but with shifts smaller than 1 TESS pixel, and the difference images from different TESS sectors give consistent results. We show the example difference images and FFIs of a target TIC 470988013 in our candidate list that we removed with d_c = 0.9 pixels in the right panel of Figure 7. All other 43 targets show a similar degree of centroid shifts like this object. We note that we keep negligible false-negative rates during this step and only exclude obvious nearby eclipsing binary systems. Figure 8 displays the centroid shift distribution of all candidates as well as the 194 remaining candidate events.

Zoom In Zoom Out Reset image size

We next remove 44 false positives through visual inspection. Among them, the alerts of 11 candidates are caused by systematics, and they show an apparent flux drop after detrending. Most of these false signals happen at the edges of TESS light-curve gaps where the flux changes sharply due to instrumental systematics. Such signals neither show transit-like shapes nor appear periodically in the light curves. A total of 31 candidates are alerted due to stars with residual stellar variations after detrending. Our uniform basis spline model failed to fully remove the stellar variability and caused the alarms. We also removed two binary systems (TIC 334790937 and TIC 446963308) with orbital period larger than 10 days and deep primary transit. We show the light curves of all these excluded candidates in Figure 15.

We note that the exact false-negative rate of this step has a negligible effect on our final result, so in our calculation of the occurrence rate, we assume that the false-negative rate (or detection completeness) for this sample of 44 stars is the same as the whole sample of ∼60,000 stars. We justify this assumption as follows: First, the false-negative rate of this sample of 44 stars is probably not significantly higher than the whole sample (although likely a bit higher due to the increased systematics or stellar variations). We have carefully examined the light curves of these candidates to minimize the possibility of missing bona fide planets. If there are additional periodic transit signals with depths of about 1% or larger in the variability or systematics, they would have been picked out easily by eye. Second, even if the false-negative rate of this sample is as high as, for example, a few times higher than the whole sample, it would still have a negligible effect on our final statistics. The sample size of 44 is very small compared to the size of the whole sample of ∼60,000, so the additional uncertainty in the false-negative rate of this sample is negligible in the equation calculating the occurrence rate, as it is a negligible fraction in the denominator of Equation (12). Therefore, we conclude that these 44 stars will not affect the final statistics.

Next, we examine the secondary eclipse signals more carefully to identify additional false positives in the candidates. Though we have already placed a constraint on the odd/even transit depth in the detection pipeline, BLS only reports the depths for two models where the period is twice the fiducial period (odd transit, 2P_BLS, T_0,BLS) and the same period but the phase is offset by one fiducial period (even transit, 2P_BLS, T_0,BLS + P_BLS). In this case, we note that: (1) Candidates with secondary eclipses happen at the fiducial period but have a half fiducial period phase shift (P_BLS, T_0,BLS + P_BLS/2) will be missed. ³⁰ (2) Our detection pipeline cannot handle cases when the odd/even depths are close to each other but different (i.e., do not satisfy the criteria we set in Section 3.2: either δ_odd/Δ or δ_even/Δ smaller than 3). Therefore, we investigate the phase-folded light curves (see Figure 16) at a suite of transit ephemerides, i.e., P_BLS, T_0,BLS; P_BLS, T_0,BLS + P_BLS/2; 2P_BLS, T_0,BLS and 2P_BLS, T_0,BLS + P_BLS. We exclude 83 targets with significant secondary eclipses (>1%) or the difference between odd and even depth (∣δ_odd − δ_even∣) is higher than 1% in this step.

We note that our odd/even vetting is unlikely to reject real planets with bona fide secondary eclipse signals. As we set a very conservative threshold on the secondary eclipse (<1%) and the odd/even difference (>1%), an imperfect detrending is unlikely to cause such a large difference between different transits, which is significant and comparable with the transit depth. We visually checked the light curves of these 83 targets and confirmed that the depth differences are astrophysical instead of systematics due to detrending issues. The photometric noise could not cause false odd/even or secondary eclipse signals as deep as 1% given the photometric precision of these 83 targets, which is all much better than 1%. Some ultrashort-period hot Jupiters have detected secondary signals (e.g., WASP-18b, Shporer et al. 2019; TOI-2109b, Wong et al. 2021). However, such secondary signals (Twicken et al. 2018) would be buried in the noise of the QLP light curves of our M-dwarf sample—planets around M dwarfs are much less irradiated by their host star compared to Sun-like stars, leading to a lower equilibrium temperature T_eq and a shallower secondary eclipse in the optical (<1 mmag) band. Therefore, our 1% secondary depth cut above is much larger than the expected <1 mmag signal, and it will keep all real planets in this step, meaning a negligible false-negative rate. Moreover, the typical standard deviation of the QLP light curves of our sample is around 3–4 mmag, which is insufficient to detect the secondary eclipse signal in our cases.

Another way of identifying false-positive signals is to compare the stellar variation periodicity with the transit signal's periodicity. Candidates with eclipse signals synchronized with out-of-transit phase variation are unlikely to be real planetary systems, because it is rare to have the stellar rotation period synchronized with the planet orbital period especially for M dwarfs. Based on the empirical relation derived by Engle & Guinan (2018), we estimate that early-type M dwarfs with rotation periods smaller than 10 days would have ages below 0.9 Gyr. However, the expected time for a planet to enforce its host star to spin at the same period with the planet's orbital period is much longer than one Hubble time (Zahn 1977). Assuming a 0.5 M_J hot Jupiter around a typical early-type M dwarf with a mass and radius of 0.5 M_⊙ and 0.5 R_⊙, the synchronization timescale can be approximated by

$$\begin{eqnarray}&&{\tau }_{\mathrm{sync}}\sim {q}^{-2}\,{\left(\displaystyle \frac{a}{{R}_{* }}\right)}^{6}\mathrm{yr}\sim {10}^{5}\,{\left(\displaystyle \frac{a}{0.05\ \mathrm{au}}\right)}^{6}\ \mathrm{Gyr},\end{eqnarray} \tag{ 3 }$$

where q = M_p /M_* is the mass ratio between the planet and star, and a/R_* is the planet orbital semimajor axis in units of the stellar radius. This is much longer than any astrophysical timescale; hence the correlation between the rotation modulation and the eclipse signal is likely due to ellipsoidal variations caused by tidal distortions and gravity brightening in stellar binaries in these cases.

In order to remove such false positives with ellipsoidal variation, we mask out all transit signals in the raw QLP light curves of each candidate and perform a Lomb–Scargle periodogram (Lomb 1976; Scargle 1982) analysis between 0.4 and 12 days to measure the stellar rotation period P_rot. We regard the highest peak of the periodogram as the rotational period P_rot. Following the methodology described in Coughlin et al. (2014), we examine the significance of the match between P_BLS and P_rot by calculating:

$$\begin{eqnarray}&&{\rm{\Delta }}P=\displaystyle \frac{{P}_{\mathrm{BLS}}-{P}_{\mathrm{rot}}}{{P}_{\mathrm{BLS}}}\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&{\rm{\Delta }}P^{\prime} =\mathrm{abs}({\rm{\Delta }}P-\mathrm{int}({\rm{\Delta }}P)),\end{eqnarray} \tag{ 5 }$$

where abs returns the absolute value, and int yields the nearest integer. This method examines and accounts for any possible period ratios between P_rot and P_BLS. We then transform ${\rm{\Delta }}P^{\prime}$ to a value quantifying the significance of the similarities between these two periods by computing the inverse of the complementary error function:

$$\begin{eqnarray}&&{\sigma }_{{P}_{\mathrm{match}}}=\sqrt{2}\times \mathrm{erfcinv}({\rm{\Delta }}P^{\prime} ).\end{eqnarray} \tag{ 6 }$$

A larger ${\sigma }_{{P}_{\mathrm{match}}}$ value means that P_BLS and P_rot are more likely to be from the same origin. Figure 9 displays the ${\sigma }_{{P}_{\mathrm{match}}}$ distribution of our sample. We remove candidates with ${\sigma }_{{P}_{\mathrm{match}}}$ greater than 2.5 (roughly corresponding to a 2.5σ significance) and with the peak in the periodogram having a false-alarm probability below 0.1%. We reduce the candidate number of our sample to 44 after this step. We note that the choice of 2.5σ is somewhat arbitrary, and the key point here is to exclude obvious eclipsing binary systems with ellipsoidal variations. We carry out an independent test with a 3σ threshold. With a stricter threshold, it requires a better match between the orbital and rotation periods, which will exclude fewer candidates. Using a 3σ cut, we find 10 new planet candidates left in the sample. However, all these additional candidates are excluded according to the planet radius, orbital period, and impact parameter cut in the final step (see Section 4.6). We thus consider that the choice of selection cut has little effect on our statistics. More importantly, the false-negative rate of setting a 2.5σ cut here will be calculated and considered in the injection and recovery tests (see Section 6).

Zoom In Zoom Out Reset image size

We next crossmatch these 44 candidates with the TESS objects-of-interest (TOI) catalog, and we find that 19 out of 44 candidates are known TOIs. Based on publicly available observational notes on ExoFOP, ³¹ we further exclude two targets from our planet candidate sample. TIC 305478010 (TOI-3580) is confirmed to be a nearby eclipsing binary through the Gaia time series (Panahi et al. 2022). Additionally, we also retire TIC 7439480 as the ground observations have confirmed that the signal is on the nearby star TIC 7439481 (TOI-4339). Although we utilize outside studies to reject false positives here, we emphasize that all candidates in our final catalog are vetted through ground-based photometric observations (see Section 5).

Finally, we derive the best physical parameters of each companion. First, we apply the celerite package (Foreman-Mackey et al. 2017) to re-detrend the raw QLP light curve by fitting a Gaussian process (GP) model with a simple Matern 3/2 kernel. The out-of-transit part of the total light curve is selected in the phase space using

$$\begin{eqnarray}&&{\rm{\Phi }}\geqslant \displaystyle \frac{\tau }{2\times {P}_{b}}+\varphi ,\end{eqnarray} \tag{ 7 }$$

where Φ represents the orbital phase of the planet candidate, and P_b and τ are the orbital period and duration. We account for the uncertainties on period P_b , mid-transit time T_0,b and duration τ by including an additional factor φ, which was set to 0.02.

After detrending, we conduct a uniform transit fit across our sample. We utilize the juliet package (Espinoza et al. 2019) to perform the fit, which employs the dynamic nested sampling approach to determine the posterior distributions of each parameter based on the dynesty package (Higson et al. 2019; Speagle 2020). The transit is modeled by batman (Kreidberg 2015). We set Gaussian priors that center around the orbital period P_b and mid-transit time T_0,b found by our detection pipeline with a width of 0.2 day. We adopt a quadratic limb-darkening law for the TESS photometry, as parameterized by Kipping (2013), as well as an informative Gaussian prior on the stellar density based on TIC v8 stellar parameters (Stassun et al. 2019). In addition, juliet makes use of the new parameterizations r₁ and r₂ to efficiently sample points in the planet-to-star radius ratio and impact parameter space (Espinoza 2018), and we place wide uniform priors on both of them. Regarding the light contamination, we set a tight truncated normal prior on the dilution factor D_TESS ³² with a standard deviation of 0.05, and allow it to vary between 0 and 1. We include a photometric jitter term to account for additional white noise in the TESS photometry and fit circular orbits with eccentricity fixed at 0. We summarize our prior settings in Table 6.

We accept candidates with orbital periods of 0.8 ≤ P_b ≤ 10 days, impact parameters of b ≤ 0.9, and companion sizes of 7 R_⊕ ≤ R_b ≤ 2 R_J, which is the parameter space of concern in this work. We end up with a final sample of nine candidates, all of which are previously alerted TOIs. Among them, four are confirmed planets, and the other five are planet candidates. Table 1 lists their basic information. The other 33 targets that do not satisfy the selection limits of R_b , P_b , and b are listed in the Appendix. Table 7 along with their exclusion reasons. We show the light curves and best-fit transit models of these 33 objects in Figure 18.

The injection is carried out as follows:

• 1.  
  We uniformly divide the period–radius space (0.8 ≤ P_b ≤ 10 days, 7 R_⊕ ≤ R_p ≤ 2 R_J) into a 5 × 5 grid;
• 2.  
  Within each cell, we draw 20 sets of physical parameters P_b , R_p as well as the impact parameter b ≤ 0.9 from uniform distributions, and randomly generate mid-transit times T_0,b between the start time of a light curve t_begin and t_begin + P_b ;
• 3.  
  We build artificial transit models using batman assuming circular orbits, during which we correct the dilution effect for each star. ³⁴ We fix the limb-darkening coefficients [μ₁, μ₂] to [0.3, 0.3] in this step for simplicity;
• 4.  
  We initialize the model at a high cadence level (100,000 points) and resample it to the real observation time stamps, and superimpose the transit model with the detrended light curve.

We randomly choose 3000 stars from the 60,382 stars in our parent sample without transit alerts and apply the above injection process. We put all synthetic light curves through our detection pipeline and record signals that are recovered. As there are inevitably variable stars in the randomly selected 3000 stars, we also require all recovered planets to pass the "synchronization" test as we did in Section 4. We did not perform the odd/even and secondary eclipse analysis because (1) the false-negative rate of this analysis in the vetting step is negligible (see Section 4.3); (2) here we only inject periodic signals with consistent transit depth.

Consequently, we insert and test 1,500,000 signals in total. We show the distribution of a total of 10,000 recovered or missed planets randomly drawn from this simulation as a function of the planet period and radius in Figure 10. Based on the fraction of recovered planets in each cell, we generate the individual sensitivity maps (${p}_{\det ,i}$) for these 3000 random stars, and combine all of them to provide an average transit detection sensitivity ($\left\langle {p}_{\det }\right\rangle$) map as in the left panel of Figure 11. We also conduct a study on injecting the planet signals into the raw QLP light curves instead of the detrended data sets, followed by spline model detrending and planet searching. We find a minor average BLS S/N decrease of ≲1. The difference in the final mean sensitivity map is around 0.003, which is within the errors we consider below (see Sections 6.2 and 7).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Before deriving the planet occurrence rate, we have to correct for the geometric probability of transit for the detectability map to find out our sample completeness. Based on Kepler's third law, the transit probability is defined as

$$\begin{eqnarray}&&{p}_{\mathrm{tr}}=0.9\displaystyle \frac{{R}_{* }}{a}=0.9{R}_{* }{\left(\displaystyle \frac{{{GM}}_{* }{P}^{2}}{4{\pi }^{2}}\right)}^{-1/3},\end{eqnarray} \tag{ 8 }$$

where R_* and M_* are the radius and mass of the star, and P is the orbital period of the companion in a circular orbit. We include a factor of 0.9 as we only take planet candidates with b ≤ 0.9 into consideration in this study. For each injected planet of every randomly selected star in Section 6.1, we compute the transit probability and multiply this factor in the detectability map to account for the geometric effect. We generate an individual completeness (${p}_{\det ,i}{p}_{\mathrm{tr},i}$) map for each star and show the resulting average ($\left\langle {p}_{\det }{p}_{\mathrm{tr}}\right\rangle$) map in the right panel of Figure 11. We rerun the injection and recovery process using another two different sets of 3000 stars and find that the differences in the completeness map are all within 0.004. We account for this uncertainty on $\left\langle {p}_{\det }{p}_{\mathrm{tr}}\right\rangle$ in the occurrence rate computation.

Compared with the occurrence rates of hot Jupiters around AFGK stars, we find that the value 0.27% ± 0.09% for early-type M dwarfs deviates from the majority of measurements. We note that our result is within the occurrence rate upper limits of Jovian-size planets around M dwarfs previously reported by Endl et al. (2006), Kovács et al. (2013), and Sabotta et al. (2021). Figure 13 shows the hot-Jupiter occurrence rates from different works as a function of the stellar type. We caution the readers that these studies use different methods (transit or RV) and have slightly different definitions for hot Jupiters. A summary of these results is presented in Table 3. After adding a measurement at the low-stellar-mass end from this work, the occurrence rate of hot Jupiters appears to have a maximum peak around G stars and decrease toward M and A dwarfs, but actually most measurements still agree with each other within 1σ–2σ so we cannot draw definitive conclusions regarding the trend in the hot-Jupiter occurrence across stellar types.

Zoom In Zoom Out Reset image size

Table 3. Summary of Occurrence Rates of Hot Jupiters from Other Works

Work	focc(%)	Stellar Type	Method	Definition of Hot Jupiters
Endl et al. (2006)	<1.27	M	RV	${M}_{p}\sin i\sim 1\ {M}_{{\rm{J}}}$, a < 1 au
Gould et al. (2006)	${0.31}_{-0.18}^{+0.43}$	FGKM	Transit	1 ≤ Rp ≤ 1.25 RJ, 3 ≤ Pb ≤ 5 days
Cumming et al. (2008)	1.5 ± 0.6	FGK	RV	${M}_{p}\sin i\gt 0.3\ {M}_{{\rm{J}}}$, Pb < 11.5 days
Mayor et al. (2011)	0.89 ± 0.36	FGKM	RV	${M}_{p}\sin i\gt 50\ {M}_{\oplus }$, Pb < 11 days
Wright et al. (2012)	1.20 ± 0.38	FGK	RV	${M}_{p}\sin i\gt 0.1\ {M}_{{\rm{J}}}$, Pb < 10 days
Howard et al. (2012)	0.4 ± 0.1	GK	Transit	8 ≤ Rp ≤ 32 R⊕, Pb < 10 days
Fressin et al. (2013)	0.43 ± 0.05	FGKM	Transit	6 ≤ Rp ≤ 22 R⊕, 0.8 ≤ Pb ≤ 10 days
Kovács et al. (2013)	<1.7%–2.0	M	Transit	Rp ∼ 1.0 RJ, 0.8 ≤ Pb ≤ 10 days
Masuda & Winn (2017)	${0.43}_{-0.06}^{+0.07}$	FGK	Transit	0.8 ≤ Rp ≤ 2 RJ, Pb < 10 days
Petigura et al. (2018)	${0.57}_{-0.12}^{+0.14}$	FGK	Transit	8 ≤ Rp ≤ 24 R⊕, 1 ≤ Pb ≤ 10 days
Zhou et al. (2019) a	0.26 ± 0.11	A	Transit	0.8 ≤ Rp ≤ 2.5 RJ, 0.9 ≤ Pb ≤ 10 days
Zhou et al. (2019)	0.43 ± 0.15	F	Transit	0.8 ≤ Rp ≤ 2.5 RJ, 0.9 ≤ Pb ≤ 10 days
Zhou et al. (2019)	0.71 ± 0.31	G	Transit	0.8 ≤ Rp ≤ 2.5 RJ, 0.9 ≤ Pb ≤ 10 days
Wittenmyer et al. (2020)	${0.84}_{-0.20}^{+0.70}$	FGK	RV	${M}_{p}\sin i\gt 0.3\ {M}_{{\rm{J}}}$, 1 ≤ Pb ≤ 10 days
Sabotta et al. (2021)	<3	M	RV	100 < Mp < 1000 M⊕, Pb < 10 days
Zhu (2022)	2.8 ± 0.8	FGK	RV	${M}_{p}\sin i\gt 0.3\ {M}_{{\rm{J}}}$, a ≤ 0.1 au
Beleznay & Kunimoto (2022) b	0.29 ± 0.05	A	Transit	0.8 ≤ Rp ≤ 2.5 RJ, 0.9 ≤ Pb ≤ 10 days
Beleznay & Kunimoto (2022)	0.36 ± 0.06	F	Transit	0.8 ≤ Rp ≤ 2.5 RJ, 0.9 ≤ Pb ≤ 10 days
Beleznay & Kunimoto (2022)	0.55 ± 0.14	G	Transit	0.8 ≤ Rp ≤ 2.5 RJ, 0.9 ≤ Pb ≤ 10 days
This work	0.27 ± 0.09	M	Transit	7 R⊕ ≤ Rp ≤ 2 RJ, 0.8 ≤ Pb ≤ 10 days

Notes.

^a Zhou et al. (2019) also reported an average occurrence rate of 0.41% ± 0.10% within their full AFG sample. ^b Beleznay & Kunimoto (2022) also reported an average occurrence rate of 0.33% ± 0.04% within their full AFG sample.

Download table as:  ASCIITypeset image

However, if this occurrence rate trend is real, it might reflect the different formation histories of hot Jupiters around different types of stars. As the mass of the protoplanetary disk scales linearly with the stellar mass M_* (Andrews et al. 2013), theoretical works predict that Jupiters are more rare around M dwarfs (Laughlin et al. 2004; Ida & Lin 2005; Kennedy & Kenyon 2008; Liu et al. 2019) due to the shortage of solid materials in protoplanetary disks to support giant planet formation. Indeed, a simulation carried out by Burn et al. (2021) shows that gas giants (M_p > 100 M_⊕) cannot form around M dwarfs with M_* < 0.5 M_⊙ through core accretion. Such drawback could, in principle, be compensated by metal-rich stars (Maldonado et al. 2020). For A-type stars, if there is indeed a drop in the occurrence rate of hot Jupiters, it could be attributed to several potential reasons. First, the rapid rotation of A-type stars and their high surface temperatures would impede giant planet detection and confirmation through spectroscopic observations. In addition, hot Jupiters may be engulfed by their host A stars (Stephan et al. 2018). Finally, the disk lifetime of A stars tends to be shorter than that of FGK stars (Ribas et al. 2015), and there could be fewer successfully formed giant planets before disk dissipation. More detections and studies on hot Jupiters around A stars are required to draw conclusions.

As can be seen from Figure 13, the occurrence rates reported by RV surveys, although consistent within 1σ, are systematically higher than the values from transit studies (see Zhu & Dong 2021 and references therein). In particular, recent work by Zhu (2022) used the Sun-like sample from the California Legacy Survey (CLS; Rosenthal et al. 2021) and measured a hot-Jupiter frequency of 2.8% ± 0.8%, which is substantially higher than the rate obtained in our work around early-type M dwarfs. Wright et al. (2012) pointed out that such a difference between the RV and transit results might be partly owing to the different stellar metallicities between these two samples. However, a further study of the Kepler stellar sample from Dong et al. (2014) shows that they have a subsolar metallicity (∼–0.04 dex) similar to the RV sample (∼0.0 dex), which implies that metallicity may have a minor impact on this discrepancy. A similar conclusion was also drawn by Guo et al. (2017). Moreover, according to the statistics from Moe & Kratter (2021), RV surveys probably increase the detection rates of hot Jupiters by a factor of 1.8 ± 0.2 by removing spectroscopic binaries among their parent samples, which could result in this feature. A promising way to test this hypothesis is to search for close stellar companions of transiting hot Jupiters with high-contrast imaging (e.g., Ngo et al. 2016), excluding circumbinary systems, and compare the remaining sample with the RV sample. Finally, stellar age may also play a role although such an effect has not been thoroughly discussed (Donati et al. 2016).

Previous research found that cold Jupiters around M dwarfs with semimajor axis a ≳ 1 au have an occurrence rate of ∼4%. Long-term RV observations from the California Planet Survey showed that the frequency is around ${3.4}_{-0.9}^{+2.2} \%$ for an M dwarf (M_* < 0.6 M_⊙) harboring planets with M_p > 0.3 M_J within 2.5 au (Johnson et al. 2010). Although Johnson et al. (2010) did not claim the inner bound of their detection limit, planet GJ 876 c in their sample has the smallest semimajor axis, about 0.13 au (Marcy et al. 2001; Rivera et al. 2005). Furthermore, various microlensing studies reported a value around 5% at 1–10 au (Gould et al. 2010; Cassan et al. 2012; Suzuki et al. 2016; Shvartzvald et al. 2016). Although RV surveys mainly focus on giant planets within 2.5 au while the microlensing method is sensitive to planets beyond the snow line, the occurrence rates of cold Jupiters measured using these two methods are consistent with each other at 1σ. For hot Jupiters located at a distance of a ≲ 0.1 au from their early-type M-dwarf hosts, we measure an occurrence rate of 0.27% ± 0.09%. Our result is significantly smaller than the frequency of outer cold gas giants, indicating that cold Jupiters are more common than hot Jupiters around M dwarfs, which is consistent with solar-like stars (e.g., Wittenmyer et al. 2020). We show the occurrence rates per semimajor axis bin obtained using different methods as a function of the semimajor axis in Figure 14.

Zoom In Zoom Out Reset image size

Combining archival data from the Anglo-Australian Planet Search (Tinney et al. 2001), Wittenmyer et al. (2020) investigated the occurrence rate of giant planets (M_p > 0.3 M_J) around solar-like stars across a wide range of semimajor axes (0.02 ≲ a ≲ 9 au). More recently, Fulton et al. (2021) also looked into the same problem using an independent sample from the California Legacy Survey, of which the orbital separation spans 0.03–30 au. The results from both works infer that the occurrence rate decreases by about 6 times from cold (1 ≲ a ≲ 10 au, $\tfrac{{dN}}{d{\mathrm{log}}_{10}(a)}\sim 19 \%$) to hot Jupiters (0.01 ≲ a ≲ 0.1 au, $\tfrac{{dN}}{d{\mathrm{log}}_{10}(a)}\sim 3 \%$). In contrast, for equivalent systems around early-type M-dwarf hosts, we find a steeper decrease, of about 14 times, for Jupiters at 1 ≲ a ≲ 10 au and 0.01 ≲ a ≲ 0.1 au, from $\tfrac{{dN}}{d{\mathrm{log}}_{10}(a)}\sim 5 \%$ to $\tfrac{{dN}}{d{\mathrm{log}}_{10}(a)}\sim 0.34 \%$ (see Figure 14). The decrease we find hints that hot Jupiters around M dwarfs may be even more difficult to form than cold Jupiters when compared with G dwarfs. However, due to large uncertainties on the occurrence rates of cold Jupiters and especially because the measurements of M dwarfs come from different methods that might have sample biases, we cannot draw firm conclusions yet. Future homogeneous near-infrared spectroscopic surveys (e.g., Mahadevan et al. 2014; Fouqué et al. 2018; Reiners et al. 2018), which perform long-term RV observations, will shed some light on this puzzle.

Both Wittenmyer et al. (2020) and Fulton et al. (2021) reported that there exists an occurrence rate transition point around 1 au for giant planets around FGK stars (see Figure 14). This jump is suggested to be relevant to the location of the snow line (Ida & Lin 2008), as an enhanced solid density beyond the snow line will facilitate the formation of the solid core under the core accretion paradigm. Due to lower irradiation, the snow line of M dwarfs is closer to the star compared with FGK dwarfs. Combining the radial velocity and microlensing findings, we can see that the occurrence rate trend of giant planets around M dwarfs is a monotonic increase as a function of the semimajor axis, similar to that of the FGK dwarfs. If a sudden increase in the occurrence rate of giant planets around M dwarfs indeed exists, the exact position of this transition is still unclear due to the limited amount of data at the moment. If future occurrence rate studies on the warm Jupiters around M dwarfs confirm the existence of a transition, it would indicate that the formation of giant planets around M dwarfs is similar to FGK stars.

During the candidate search, we found seven new planet candidates that were not previously alerted as TOIs. They pass all vetting steps including centroid analysis, visual inspection for odd/even, and secondary signals as well as synchronization test. All of these candidates have a radius below 7 R_⊕ so they are not included in the statistical sample. We summarize their properties in Table 4. As we only report the results from our uniform fits (see Section 4), we emphasize that more detailed analyses are required to evaluate the robustness of these candidates but such works are beyond the scope of this study. Note that during the writing of this manuscript, TIC 291109653 was alerted by the QLP faint-star search program (Kunimoto et al. 2022) through a sector-combined analysis (Sectors 23 and 46), designated as TOI-5486. In this work, we independently find the signal using sector 23 data.

Table 4. New Planet Candidates Detected by Our Pipeline that were not Announced as TOIs Before

TIC	Tmag	Period (days)	T0 (BJD-2457000)	Rp (RJ)	Centroid Shift (pixels)	${\sigma }_{{P}_{\mathrm{match}}}$
32296259	12.28	2.77757	1492.9043	0.42	0.19	0.99
101736867	13.09	2.64795	1655.0761	0.59	0.83	0.72
115524526	12.94	4.65712	1956.4309	0.45	0.61	1.99
246974219	12.29	1.90943	1793.5922	0.42	0.48	0.96
291109653 a	12.29	2.02479	1929.7720	0.33	0.17	0.82
367411575	13.25	1.19342	1792.7369	0.57	0.48	1.91
371315491	13.26	0.40622	1571.8947	0.55	0.78	1.26

Note. The physical parameters come from our uniform transit fit.

^a TIC 291109653 was recently alerted as TOI-5486 by the QLP faint-star search program (Kunimoto et al. 2022).

Download table as:  ASCIITypeset image

We obtained two focused ground-based follow-up observations for TIC 20182780 (TOI-3984) using the 1 m Las Cumbres Observatory Global Telescopes (LCOGT; Brown et al. 2013) on 2022 April 14 and June 6 in the $i^{\prime}$ and $g^{\prime}$ bands with exposure times of 150 and 300 s. The Sinistro cameras have a $26^{\prime} \times 26^{\prime}$ field of view (FOV) as well as a plate scale of 0389 per pixel. The images of both observations have stellar point-spread-functions (PSF) with FWHMs of 22 and 23, respectively. We also acquired focused alternating V&I-band observations using the 1.5 m telescope at the Observatorio de Sierra Nevada on 2022 May 10, which has an FOV of $18^{\prime} \times 13^{\prime}$ and a pixel scale of 0455. The exposure times we set are 100 s and 50 s for V- and I-band observations. The estimated PSF is about 27.

In addition, we took an AO observation for TIC 20182780 using the Palomar High Angular Resolution Observer (PHARO) on the Palomar-5.1 m telescope on 2022 February 13 in the Br-γ band (λ_o = 2.2 μm) to search for stellar companions. The result reveals that it is an isolated star, with no companions 6.1 mag fainter than the target out to 05.

To constrain the companion mass, we also obtain seven RVs with the WIYN/NEID spectrograph (Evans et al. 2016) between 2022 March 13 and 2022 April 9 with a baseline of 27 days. The observations are queue scheduled. Each exposure took 1200 s, and the median RV precision is about 25 m s⁻¹. The NEID RVs put a 3σ mass upper limit of 0.32 M_J on the companion mass and rule out the stellar binary scenario. We refer the readers to Wang et al. (2022) for more information about NEID observation and data reduction. The final RV data are listed in Table 8.

We collected two ground focused photometric observations for TIC 71268730 (TOI-5375) using the GdP-0.4 m and CMO-0.6 m at Grand-Pra and Caucasian Mountain Observatory. The GdP-0.4 m has an FOV of $12\buildrel{\,\prime}\over{.} 9\times 12\buildrel{\,\prime}\over{.} 5$ and a pixel scale of 073. The observation was taken on 2022 March 5 in a clear filter with an exposure time of 180 s. The seeing is good with a light-curve rms of 0.0074. The CMO-0.6 m has an FOV of $22^{\prime} \times 22^{\prime}$ and a pixel scale of 067 (Berdnikov et al. 2020). The observation was taken on 2022 March 31 in R_c band with an exposure time of 120 s. The PSF of two observations are 45 and 24. We used 10 and 7 comparison stars with an aperture of 66 and 47 to do the photometric analysis for these two observations, respectively.

We collected three LCOGT light curves for TIC 79920467 (TOI-3288), one of them was done with the 0.4 m while two were done with the 1 m. LCOGT-0.4 m has an FOV of $29^{\prime} \times 19^{\prime}$ with a pixel scale of 0571. The LCOGT-0.4 m observation was carried out in the $i^{\prime}$ band on 2021 June 7 with an exposure time of 200 s. The LCOGT-1 m observations were carried out in the $i^{\prime}$ and $g^{\prime}$ bands on 2021 June 19 and 2022 May 16 with an exposure time of 100 s and 300 s. The images are all focused with a PSF of 26, 23, and 20, respectively. We used an aperture of 51 to reduce the LCOGT-0.4 m data while 16 for the LCOGT-1 m data. We also obtained two luminous band observations for this target using a CDK20-0.5 m at El Sauce Observatory, Chile on 2021 September 2 and 2021 October 28, under a good seeing condition. For these observations, the exposure time was set at 120 s. The CDK20-0.5 m has an FOV of $35\buildrel{\,\prime}\over{.} 87\times 35\buildrel{\,\prime}\over{.} 87$ and a pixel scale of 052.

We acquired a total of five focused LCOGT-1 m observations for TIC 95057860 (TOI-4201). The first observation was done in the $i^{\prime}$ band on 2021 September 1. Two $g^{\prime}$- and $i^{\prime}$-band alternating observations (four light curves) were carried out on 2021 September 26 and 2021 October 13. All $i^{\prime}$-band observations were taken with an exposure time of 180 s while $g^{\prime}$-band observations have a 300 s exposure time. The photometric apertures we used for the three observations are 85, 85, and 62. The PSFs of the three observations are 23, 51, and 33, respectively.

We collected two ground-based follow-up light curves for TIC 382602147 (TOI-2384). The first observation was obtained with the Evans telescope at the El Sauce Observatory, Chile, a 0.36 m Corrected Dall Kirkham, in the R_c band on 2020 November 9. The telescope was fitted with an SBIG 1603-3 CCD with 1536 × 1024 pixels binned 2 × 2 in camera for an image scale of 147 pixel⁻¹, giving a field of view of $18\buildrel{\,\prime}\over{.} 8\times 12\buildrel{\,\prime}\over{.} 5$. The calibrated data consisting of 105 exposures of 180 s were analyzed with a circular aperture of 59 radius in AstroImageJ. Another observation was done with LCOGT-1 m in the $g^{\prime}$ band on 2021 August 5. We reduced the data with a 47 circular aperture. The PSFs of the two observations are 33 and 24, respectively.