The occurrence rate of giant planets orbiting low-mass stars with TESS

ABSTRACT

We present a systematic search for transiting giant planets (⁠|$0.6 \mbox{$R_{\rm J}$}\le \mbox{$R_{\rm P}$}\le 2.0 \mbox{$R_{\rm J}$}$|⁠) orbiting nearby low-mass stars (⁠|$\mbox{$M_{*}$}\le 0.71 \mbox{${\rm M}_{\odot }$}$|⁠). The formation of giant planets around low-mass stars is predicted to be rare by the core-accretion planet formation theory. We search 91 306 low-mass stars in the TESS 30 min cadence photometry detecting fifteen giant planet candidates, including seven new planet candidates which were not known planets or identified as TOIs prior to our search. Our candidates present an exciting opportunity to improve our knowledge of the giant planet population around the lowest mass stars. We perform planet injection-recovery simulations and find that our pipeline has a high detection efficiency across the majority of our targeted parameter space. We measure the occurrence rates of giant planets with host stars in different stellar mass ranges spanning our full sample. We find occurrence rates of 0.137 ± 0.097  per cent (0.088–0.26 M_⊙), 0.108 ± 0.083 per cent (0.26–0.42 M_⊙), and 0.29 ± 0.15 per cent (0.42–0.71 M_⊙). For our full sample (0.088–0.71 M_⊙), we find a giant planet occurrence rate of 0.194 ± 0.072 per cent. We have measured for the first time the occurrence rate for giant planets orbiting stars with |$\mbox{$M_{*}$}\le 0.4\, \mbox{${\rm M}_{\odot }$}$| and we demonstrate this occurrence rate to be non-zero. This result contradicts currently accepted planet formation models and we discuss some possibilities for how these planets could have formed.

planets and satellites: detection, planets and satellites: formation, stars: low-mass

1 INTRODUCTION

Since the detection of 51-Pegasi b (Mayor & Queloz 1995) revealed the existence of giant planets on short orbits (P ≲ 10 d), known as hot Jupiters, the prevalence of these planets has been a question of great interest to the community. As such, there have been numerous studies that have focused on measuring the occurrence rate of these planets. Different planet formation mechanisms should leave signatures in the occurrence rates of the different populations of planets (e.g. Emsenhuber et al. 2021a,b; Schlecker et al. 2021). For example, it has been demonstrated that gas giant planets are found predominantly around high metallicity host stars (Fischer & Valenti 2005; Osborn & Bayliss 2020). It is believed that the metallicity of a star is inherited from its primordial cloud and thus shared by its protoplanetary disc (Fischer & Valenti 2005), and it has been shown that gas giant planets from through core-accretion more efficiently in these metal-rich protoplanetary discs (e.g. Ida & Lin 2004; Thommes, Matsumura & Rasio 2008). As such, this giant planet metallicity correlation has been interpreted as evidence that giant planets predominantly form through the core-accretion mechanism (e.g. Johnson et al. 2010).

Over the past two decades, the number of known exoplanets has exploded, particularly as a result of the Kepler mission (Borucki et al. 2010), resulting in huge advances in our knowledge of hot Jupiter occurrence rates. Using early Kepler results for a sample of G- and K-type stars (T_eff = 4100–6100  K and log g = 4.0–4.9), Howard et al. (2012) derived a frequency of 0.4 ± 0.1 per cent for planets with |$\mbox{$R_{\rm P}$}\ge 8$| R_⊕ and P < 10 d. Improving upon the estimates of the rate of false positives in the Kepler data, Fressin et al. (2013) refined the Kepler occurrence rate estimate to 0.43 ± 0.05 per cent for planets with |$6\, \mbox{$R_{\hbox{$\oplus $}}$}\le \mbox{$R_{\rm P}$}\le 22\, \mbox{$R_{\hbox{$\oplus $}}$}$| and 0.8 ≤ P ≤ 10 d.

Despite providing a deep insight into the occurrence rates of planets around FGK stars, Kepler has not been able to provide similar results for the occurrence rate of giant planets orbiting low-mass stars. Kepler only monitored roughly 4000 stars with T_eff|$\le 4000\,$| K (Dressing & Charbonneau 2013). With the detection of 91 transiting planet candidates with |$\mbox{$R_{\rm P}$}\le 4.0\, \mbox{$R_{\hbox{$\oplus $}}$}$| around these stars. Dressing & Charbonneau (2013) could strongly constrain the occurrence rates of small radii planets with low-mass host stars and periods under 50 d to be |$90^{+4}_{-3}$| per cent. Unfortunately, the low number of M-dwarfs monitored by Kepler lead to the discovery of just three planet candidates with |$\mbox{$R_{\rm P}$}\gt 0.6\mbox{$R_{\rm J}$}$| (Dressing & Charbonneau 2013). Based on their three candidates Dressing & Charbonneau (2013) derived an occurrence rate for giant planets on short orbital periods (⁠|$P \le 10\, d$|⁠) of |$1^{+0.8}_{-0.2}$|  per cent. However, of these three candidates just one – Kepler-45 b (Johnson et al. 2012) – has been confirmed as a giant planet. We note that the other two have since been shown to be a |$1.97\, \mbox{$R_{\hbox{$\oplus $}}$}$| planet (Kepler-1124 b; Morton et al. 2016; Cloutier & Menou 2020) and an M-dwarf–white dwarf binary system (KOI-256; Muirhead et al. 2013). If we instead consider now the fact that, from these three candidates, just Kepler-45 b is real, then from the Dressing & Charbonneau (2013) results we would obtain an occurrence rate of |$0.3^{+0.44}_{0.1}$|  per cent, which is statistically consistent with the hot Jupiter occurrence rate measured by Kepler (e.g. Fressin et al. 2013; Petigura et al. 2018). Therefore, the Kepler results alone are unable to conclude that the giant planet occurrence rate for low-mass stars is significantly lower than for solar-like stars.

Determining the frequency of giant planets around low-mass stars will act as a key test for planet formation theories. It has been shown that giant planets form less readily through core accretion around low-mass stars than Sun-like stars (Laughlin, Bodenheimer & Adams 2004; Ida & Lin 2005). Theoretical population synthesis models from Burn et al. (2021) predict that the occurrence rate of giant planets (⁠|$\mbox{$M_{\rm P}$}\ge 30$| M_⊕) with 0.7 M_⊙ host stars would be approximately half that for a 1 M_⊙ host star. For 0.5 M_⊙ host stars they predict an occurrence rate roughly 13 per cent that of a 1 M_⊙ host star, and for host stars with M_* < 0.5 M_⊙ their results find giant planet occurrence rates of zero, predicting that the formation of giant planets through core-accretion is impossible for these low-mass host stars.

Recent discoveries of hot Jupiter planets orbiting M-dwarfs (e.g. Hartman et al. 2015; Bayliss et al. 2018; Bakos et al. 2020) have shown that giant planets can form around low-mass stars. However, despite these recent discoveries, there are still very few known hot Jupiters transiting low-mass stars and the overall population is still not well understood. From the 536 transiting giant exoplanets known to date, only sixteen orbit a star with a stellar mass less than 0.65 M_⊙, and only one – HATS-71 b (Bakos et al. 2020) – orbits a star with a stellar mass less than 0.5 M_⊙¹ (see Fig. 1). The true impact of these discoveries on our knowledge of planet formation cannot be determined without a robust measurement of the occurrence rate of these systems.

Figure 1.

Known transiting exoplanets as a function of the mass of the host star (data accessed from the NASA Exoplanet Archive on 2022 November 16). The sample shown is made up of planets with masses between 0.1 and 13 M_J and radii larger than 0.6 R_J.

Attempts have been made to measure the occurrence rates for giant planets with low-mass host stars, but the majority have not yielded conclusive results. The Pan Planets survey (Obermeier et al. 2016) placed a 95 per cent confidence level upper limit of 0.34 per cent on the occurrence rates for hot Jupiters with host stars with an effective temperature T_eff < 3900 K and surface gravity log g > 4, concluding that they could not confirm a significant difference in the occurrence rate of hot Jupiters compared to those with solar-type host stars. More recently, Sabotta et al. (2021) attempted to derive planet occurrence rates around stars with masses M_* ≤ 0.7 M_⊙using data from CARMENES. This study suffers from a small stellar sample size, including data for only 71 low-mass stars, and did not detect any giant planets with P ≤ 10 d. Therefore, they could only place an upper limit of 3 per cent on giant planet frequency orbiting low-mass stars. Gan et al. (2022) have used the TESS photometry to constrain the occurrence rate of hot Jupiters (⁠|$7 \mbox{$R_{\hbox{$\oplus $}}$}\le \mbox{$R_{\rm P}$}\le 2 \mbox{$R_{\rm J}$}$|⁠; 0.8 ≤ P ≤ 10 d) around early M-dwarfs (⁠|$0.45-0.65 \mbox{${\rm M}_{\odot }$}$|⁠) to 0.27 ± 0.09 per cent.

In this paper, we present a systematic search for transiting gas giant planets (⁠|$0.6 \mbox{$R_{\rm J}$}\le \mbox{$R_{\rm P}$}\le 2.0 \mbox{$R_{\rm J}$}$|⁠) orbiting nearby (d ≤ 100 pc), low-mass (⁠|$0.088\, \mbox{${\rm M}_{\odot }$}\le \mbox{$M_{*}$}\le 0.71\, \mbox{${\rm M}_{\odot }$}$|⁠) stars in the TESS 30-min cadence light curves (Caldwell et al. 2020). This search has been motivated and designed to measure the occurrence rates of giant planets on close-in orbits (⁠|$1\, {\rm d} \le P \le 10\, {\rm d}$|⁠) around low-mass stars.

Through this study, we have measured giant planet occurrence rates for lower mass host stars than previous studies. In addition to measuring these occurrence rates, our search has also yielded the discovery of seven new giant planet candidates which were not previously known, either as confirmed planets or TOI planet candidates, for which follow-up work to confirm their true natures is underway. We present our new candidates in this work and we interpret our results in the context of giant planet formation.

We detail our sample selection in Section 2. We present our transit search in Section 3, the automated light-curve vetting in Section 4, and the transit fitting analysis in Section 5. We discuss some final vetting of our candidates, including identifying nearby blend scenarios, in Section 6. Our giant planet candidates are presented in Section 7 and discuss some preliminary spectroscopic follow-up in Section 7.1. We discuss the injection-recovery simulations used to determine the detection efficiency of our pipeline in Section 8 and present our occurrence rate measurements in Section 9, comparing them to previous occurrence rate studies.

2 TESS LOW-MASS STAR SAMPLE

Launched in 2018 the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) has been surveying the whole sky for transiting planets, providing us with a new way to study the planetary occurrence rates for different types of host star. TESS observes the sky in sectors. During each sector a 24 degree × 96 degree area of the sky is observed for 27.4 d. Across the primary mission the Northern and Southern hemispheres were observed for thirteen sectors each. These sectors overlapped, especially around the ecliptic poles, and so while most objects receive one or two sectors of observations, objects close to the ecliptic poles can receive up to 13 sectors of coverage.

The primary science goal of the TESS mission is to discover small radius planets – |$\mbox{$R_{\rm P}$}\le 4$| R_⊕–amenable to further characterization including with JWST (Gardner et al. 2006). The discovery of these small exoplanets was primarily achieved through 2-min cadence light curves produced for a selected sample of stars from each sector. In addition to these high cadence light curves, the TESS full-frame-images (FFIs) are supplied at a longer cadence. During the primary mission this cadence was 30 min. These FFIs allow for the production of long cadence light curves for large samples of stars across the whole sky. These large datasets allow TESS photometry to be used for occurrence rate studies.

We use the 30-min cadence light curves from Years 1 and 2 (sectors 1–26) of the TESS mission produced from the TESS FFIs by the TESS-SPOC team (Caldwell et al. 2020), which we access from the Mikulski Archive for Space Telescopes (MAST) as a High Level Science Product² These light curves are selected as they provide a high quality reduction and a light-curve data product that is largely free of instrumental systematics. In addition, the target selection criteria used by the TESS-SPOC team prioritises low-mass stars such as those we are targeting in this study. The red-sensitive design of the instrument enables TESS to achieve high precision photometry for M and K spectral type stars. Therefore, in this work we seek to use these TESS data to place the best constraints to date on the occurrence of giant planets orbiting low-mass stars.

To select a statistical sample of low-mass stars, we used the TESS Input Catalog version 8 (TICv8; Stassun et al. 2019). The Gaia DR2 parameters are incorporated into TICv8 allowing for low-mass stars to be more reliably identified and our sample to be free of giant stars. We select stars for our sample which meet all the following criteria:

• T_eff ≤ 4500  K,

• R_* ≤ 0.75 R_⊙,

• d ≤ 100  pc, and

• T ≤ 16  mag,

where T_eff is the effective temperature of the star, R_* is the stellar radius, d is the distance of the star from the Sun, and T is the stellar apparent magnitude in the TESS filter. The distance cut is imposed as the TESS-SPOC target selection process uses the same selection cut. We find a total of 168 837 stars in the TIC that satisfy all these criteria, and using tess-point (Burke et al. 2020) we identify that 121 402 of those stars were observed during the TESS Prime Mission. Cross-matching our full TIC sample with the target lists for the TESS-SPOC 30-min cadence FFI light curves yields 91 306 stars with a light curve. We present our low-mass star sample in Fig. 2. From this plot, as well as the histogram in Fig. 3, we see that our low-mass star sample spans a wide range of spectral types and stellar masses (≈0.1–0.7 M_⊙), with the majority of a sample being lower mass than 0.4 M_⊙ and so expected to be incapable of forming giant planets (e.g. Burn et al. 2021). Therefore our stellar sample is ideal for studying the extremes of giant planet formation and investigating the lowest mass host stars around which it is possible to form giant planets. The stars which were observed by TESS but do not have a TESS-SPOC FFI light curve are preferentially faint (T ≲ 14 mag) and small (⁠|$\mbox{$R_{*}$}\ \lesssim 0.25\, \mbox{$R_{\odot }$}$|⁠) stars. While it is important to be aware of this selection we do not expect this to have a substantial impact on our occurrence rate measurements.

Figure 2.

Top: Gaia colour–magnitude diagram showing the stellar type distribution of the 91 306 low-mass stars in our sample (density heat map). Bottom: Histogram showing the Gaia colour distribution of the sample.

Figure 3.

Histograms displaying the apparent magnitude (in the TESS bandpass), stellar mass, and distances for the 91 306 stars in our low-mass star sample. We also show the number of TESS sectors for which as TESS-SPOC 30-min cadence light curve exists for each object, plotting this histogram on a log scale.

From the histogram in Fig. 2, we see that the majority of the stars in our sample have the spectral type M2–M4. This distribution arises from the apparent magnitude and distance criteria, as well as the distribution of stars across the sky. The left-hand side of the spectral type distribution follows the mass function for stars and is a result of us using a volume-limited sample and later spectral type stars being more common in the Galaxy (Chabrier & Baraffe 2000). The turnover and decrease in numbers to very late spectral types is a result of the magnitude limit of our sample (T ≤ 16 mag) and the very low luminosity of M5 and later spectral types stars.

The apparent magnitude distribution of our sample peaks around T = 13 mag (see top-left histogram in Fig. 3). With a magnitude-limited sample, we would expect this distribution to continue to rise until the magnitude limit. The turnover arises as a result of the distance limit imposed on our sample. Early M-dwarf stars with apparent magnitudes of T ≳ 13.5 mag have distances >100 pc. The distribution of the number of observing sectors for each star, shown in the bottom right histogram in Fig. 3, arises from the TESS observing strategy. As can be seen from the TESS observations footprint³, the majority of stars receive one or two sectors of coverage. Moreover, due to the overlapping of sectors in the continuous viewing zones surrounding the Ecliptic poles, a larger fraction of stars receive ≥11 sectors of coverage than receive between 6 and 10 sectors of observations.

3 TRANSIT SEARCH

Using the TESS-SPOC 30-min cadence light curves, we firstly exclude any data points with a quality flag q > 0. These flags highlight any observations which are adversely affected by cosmic rays and scattered light, among others. For this analysis we use the PDC_SAP flux time series, which has been largely corrected for instrumental systematics but for which real stellar variability has been preserved (Smith et al. 2012; Stumpe et al. 2012, 2014). To remove stellar variability prior to the planet search we smooth the light curves using a Savitzky–Golay filter, using the implementation from the lightkurve package (Lightkurve Collaboration 2018). The filter uses a window length of 45 data points, approximately 1 d for 30-min cadence data. This window length is long enough to preserve the transits produced by our target systems, which have durations ≪1 d.

We then search through the TESS-SPOC light curves for periodic planet candidate transit signals using the Box-fitting Least-Squares algorithm (BLS; Kovács, Zucker & Mazeh 2002). We use the astropy implementation of the algorithm.⁴ We searched for signals with periods between 1 and 10 d and durations between 0.03 and 0.3 d. We identify objects with a Signal Detection Efficiency (SDE; see Kovács et al. 2002) ≥8 and a transit signal-to-noise ratio (SNR) of greater than 8 as good transiting planet candidates. In total, our search yields the detection of 3930 periodic transit-like signals.

4 FALSE POSITIVE IDENTIFICATION

Along with planet candidates, our BLS detections also included a large number of false positives. These can be either real astrophysical false positive scenarios, such as eclipsing binaries or variable stars, or spurious detections. We performed a number of vetting checks to identify these false positive detections and exclude them from our candidate list. Previous transit searches have used similar criteria to identify false positives (e.g. Montalto et al. 2020). For each of the checks used to identify false positives we employ a strict quantitative cut-off to distinguish between false positives and planetary candidates. Applying these checks to simulated planet signals (see Section 8), we find that each check falsely identifies <1 per cent of the simulated planets as false positives. Therefore we are confident that these checks are not excluding a significant number of real planets from our candidate list. To summarize the results of these checks we identify 1828 clear false positives, leaving us with 2102 remaining candidates. Further details on the checks performed and the number of false positives identified by each check are provided in the following section. Note that in the pipeline these checks are performed sequentially in the order they are presented here, and any objects identified as false positives by one check are not passed to the following checks. So the numbers of identified false positives do not include any false positives identified by previous checks.

4.1 Secondary eclipse events

One of the main astrophysical false positive cases for transiting planet searches are eclipsing binary stars (e.g. Almenara et al. 2009; Santerne et al. 2012). The main signature of an eclipsing binary is the presence of a secondary eclipse. For a circular orbit, this eclipse is present at phase 0.5, however, for an eccentric orbit this is not always the case. We search for secondary eclipse events present at orbital phases between 0.2 and 0.8. To identify the secondary eclipse events we use a window of width equal to the duration of the BLS transit model and estimate the depth of any eclipse occurring within the window. The median value of the flux values within the window is taken as the depth of the eclipse. We compare the difference between the median flux level inside the window and the median overall out-of-transit flux baseline to the overall RMS scatter of the light curve, σ_LC, to assess the significance of any detected signals. While designed to look for eclipsing binaries, this check also identifies a large number of variable stars. This check has the possibility of identifying real planets as false positives in two scenarios. The first is a planet which is detected by BLS at the wrong orbital period, usually at a multiple of two or one and a half times the orbital period. However, while some real planets may be rejected by this check, when considering our injected transiting planets (see Section 8) less than 1 per cent of the simulated planets recovered by BLS are detected at incorrect orbital periods which would trigger this check. As such, we do not expect to have missed many, if any, real planets due to this. Moreover, all these effects are captured by the injection-recovery tests (Section 8). The second scenario is where the planet itself has a detectable secondary eclipse when it is occulted by the star. However, these planetary eclipses are typically very shallow. As an example, for a 1 R_J planet in a 1.5 d orbit around a 0.4 M_⊙ star, using the equations from Alonso (2018) we estimate a planetary secondary eclipse depth of just 320 ppm. This is less than the RMS scatter level for the vast majority of our stellar sample. Therefore, while it is important to be aware that these two scenarios could trigger this check, these would be very rare occurrences and we do not expect them to significantly affect our occurrence rate results.

Any object which displays a flux difference greater than 2.0 times σ_LC, occurring between phase values of ϕ_Sec = 0.47 and ϕ_Sec = 0.53 is taken to be an eclipsing binary. For the other phases searched for secondary eclipse signals, a minimum flux difference of 3.5 times σ_LC is required to identify the object as a false positive. A total of 246 objects are identified as false positives by this check.

4.2 Odd-even depth differences

Eclipsing binaries can also result in BLS detections with a difference in-depth between the odd and even transit events. This occurs when the BLS algorithm folds the light curve on half the true period of the binary, resulting in the primary and secondary eclipses both falling at phase 0. To check for such cases, we compare the depth of the odd transits with the depth of the even. The transit depths were estimated by comparing the median of the lowest 20 per cent in-transit flux values to the median overall out-of-transit flux baseline. Any objects with an odd-even depth difference greater than 3.5 σ_LC are identified as likely eclipsing binaries. Poorly sampled transit events can lead to inaccurate depth estimates, causing real transiting planets to be mis-identified as eclipsing binaries. To avoid this, we only apply this check to objects with at least ten in-transit points for both the odd and even transit events. A total of 118 objects are identified as likely eclipsing binaries by this check.

4.3 Sector depth differences

In some cases, the depth of the transit events is observed to be different in different sectors, but constant in a given sector. We identified this to be as a result of differing levels of contamination from a nearby eclipsing binary. For different sectors, the positions of neighbouring stars relative to a target star will change due to the rotation of the TESS spacecraft (see Fig. 4). Moreover, in different sectors a given star will fall in a different location in the camera field of view and so will have a differently shaped point-spread-function. The pixel mask used by the TESS-SPOC pipeline to extract photometry also varies from sector to sector. The result of these effects is that the level of contamination from neighbouring stars into the target pixel aperture is different for different sectors. The TESS-SPOC pipeline computes this dilution value and corrects the target light curve for it, so the depths of transits on the target stars should be unaffected by this. However, if one of these neighbouring stars is itself an eclipsing binary then in sectors in which it contributes more light to the target aperture the eclipse event seen in the light curve will be deeper. We present an example of such a scenario in Fig. 4.

Figure 4.

Left column: Individual sector light curves for TIC-233684011, which was identified as a nearby eclipsing binary through the difference in eclipse depth observed for different sectors. We plot the spline-smoothed PDC_SAP light curves phase folded on the best BLS period of P = 3.64 d. Right column: Target pixel cutouts for TIC-233684011 for each sector plotted. The red boxes denote the pixels included in the target aperture used in each sector.

To identify any clear nearby eclipsing binary cases, we compared the depths of the transit events between sectors, for those objects which were observed in more than one sector. The transit depths were calculated using the same method as in Section 4.2. Any object in which a depth difference between sectors of greater than 3.5 σ_LC was identified as a nearby eclipsing binary. A total of 44 objects are identified as likely nearby eclipsing binaries by this method and an example of such an object is shown in Fig. 4.

4.4 Transit phased variability

In addition to secondary eclipses and eclipse depth differences, many eclipsing binaries show out-of-transit variability in phase with the eclipses. This variability can arise from a number of effects, including ellipsoidal modulation of reflection from the an orbiting body (Faigler & Mazeh 2011). While these effects do occur for close orbiting planets, the amplitude of the resulting variability is much larger for eclipsing binaries.

We search for these variability signals in the pre-flattened light curves as we expect them to be removed by our pre-transit search light-curve flattening. We also use the results from the BLS search to mask out the transit events and perform this analysis only on the out-of-transit data. We use the harmonic model presented by Montalto et al. (2020)

$$\begin{eqnarray} f_{\rm Harmonic} = A\cos \left(\frac{2\pi \, t}{P}\right) + B\sin \left(\frac{2\pi \, t}{P}\right) + C \end{eqnarray}$$

(1)

and fit for the coefficients A, B, and C for each light curve. We then compare the quality of this harmonic fit to a simple flat line fit by computing the value

$$\begin{eqnarray} R^2_{\rm Harmonic} = 1-\frac{\Sigma ^N_{i=1}\left(f_i-f_{\rm Harmonic, {\it i}}\right)^2}{\Sigma ^N_{i=1}\left(f_i-C\right)^2} . \end{eqnarray}$$

(2)

The better the harmonic fit is to the data, the closer to 1 this value becomes. As with the secondary eclipse check (Section 4.1) this check finds purely variable stars, as well as eclipsing binaries. We identify any object with |$R^2_{\rm Harmonic} \gt 0.5$| as an eclipsing binary or a variable star, and in total identify 231 false positives using this method. An example of an eclipsing binary identified using this check is shown in the top left panel of Fig. 5.

Figure 5.

Examples of false positives identified by the variability analysis checks detailed in Sections 4.4–4.5.3. The titles of each panel give the object plotted as well as the vetting check that identified the object as a false positive. Note that for TIC-4957914 the pre-flattened light curve is plotted, but for the rest the light-curve smoothing has been applied. All objects are phase folded using the BLS period.

4.5 Variable stars

Besides eclipsing binary stars, the other major astrophysical contaminant in our BLS detections are variable stars. This stellar variability arises from two main mechanisms. The first is the combination of active regions on the stellar surface and the rotation of the star. These active regions are either regions which are less luminous than the majority of the stellar surface, which are called spots, or more luminous, known as faculae. These spots and faculae rotate into and out of view as the star rotates, resulting in a fluctuation in the observed brightness of the star (Boisse, Bonfils & Santos 2012). Periodic variability in the brightness of stars can also arise from the radial oscillation of the outer layers of the star, known as stellar pulsations (Cox 1980). Multiple types of pulsating stars are found, including δ-Scuti (e.g. Rodríguez & Breger 2001) and RR-Lyrae (e.g. Simon & Teays 1982), and can exhibit variability with an amplitude of the order of a magnitude. While our low-mass target stars themselves will not exhibit these pulsations, if a nearby star that is contaminating the photometric aperture is such a pulsator then it can imprint such variability on to the target light curve.

4.5.1 Light-curve symmetry

The phase-folded light curve of an exoplanet transit will be symmetric around phase of zero. Variable star light curves which trigger BLS will often not share this symmetry. As such, by determining the symmetry of the folded light curve we can identify variable stars. We calculate the mean point-to-point scatter of the folded light curve in two configurations. The first is a standard fold running from phase −0.5 to 0.5. The second runs from phase 0.0–0.5, where we have taken the absolute value of negative phases. We calculate the ratio of these two values and identify any objects with a ratio greater than 2.5 as variable stars. An example of a variable star identified through the asymmetry of its light curve is shown in the top right panel of Fig. 5. A total of 8 objects are identified as variable stars using this analysis.

4.5.2 Lomb–Scargle power

We also check for continuous variability by performing a Generalized Lomb–Scargle (Lomb 1976; Scargle 1982) analysis on the light curve, after masking out the flux data points which are detected as ‘in-transit’ by BLS. Any star with a normalized Lomb–Scargle power of greater than 0.3 is identified as a variable star and excluded from our candidate list. An example of such a variable star is shown in the bottom left panel of Fig. 5. A total of 1039 objects are identified as variable stars using this test. This check alone removes 26.4 per cent of the transiting candidates detected by the BLS search. However when applied to the simulated planet light curves just 0.22 per cent of the simulated planets are excluded by it. Therefore, we can be confident that we are not excluding a large number of real planets with this analysis.

4.5.3 Excess standard deviation scatter metric

A further way to identify stellar variability is to compare the RMS scatter of the phase-folded light curve to the RMS of the point-to-point scatter. For a transiting planet light curve, the out-of-transit section of the light curve will be close to flat and dominated by Gaussian noise. Therefore, these two values will differ by a factor of |$\sqrt{2}$|⁠, with the raw RMS being the smaller of the two. For a variable star, the RMS scatter will be significantly larger than the point-to-point scatter, as the RMS will be dominated by the variability. Therefore, we calculate following metric

$$\begin{eqnarray} \sigma _{\rm excess} = \frac{1}{\sigma _{\rm LC}\sqrt{2}}\sqrt{\frac{\Sigma _{1}^{N} (f_{i}-f_{i-1})^2}{N}} , \end{eqnarray}$$

(3)

where σ_LC is the RMS scatter and the square-root term gives a measure of the mean point-to-point scatter of the flux. We identify any star with σ_excess ≤ 0.5 as a variable star. An example of a variable star identified using this check is shown in the bottom right panel of Fig. 5. A total of 12 objects are identified as variable stars through this analysis.

4.5.4 Depth metric

As well as astrophysical false positives, our BLS sample also contains a few spurious detections. The most common of these spurious detections are those which are driven by the presence of a small number of outlying points. These outliers can be produced by sharp changes in the scattered light on the camera, or from a decrease in the quality of the spacecraft pointing. To identify these signals, we calculate the median flux level of the ‘in-transit’ data points of the BLS event and compare this value to the reported BLS depth. The BLS depth is determined by a fit to the data, and so is biased to larger values by the outlying points. For real transit events, these two values will be comparable, however, for spurious detections being driven by a few outlying points the BLS depth will be significantly larger than the median flux level. We calculate the ratio of the median flux level to the BLS depth and exclude any object with this ratio less than 0.5 as a spurious event. A total of 130 objects are identified as likely spurious events using this method.

5 TRANSIT FITTING ANALYSIS

Having identified and rejected the false positives, we then fit the transit events for the remaining candidates. We use the emcee package (Foreman-Mackey et al. 2013) to perform a Markov chain Monte Carlo (MCMC) analysis.

The free parameters we use in our analysis are: a reference mid-transit time, T_C, the planetary orbital period, P, the planet-to-star radius ratio, R_P/R_*, the scaled semimajor axis, a/R_*, and the orbital inclination, i. We also fit for a free flux baseline offset, f₀, which is defined such that the out-of-transit flux is equal to 1 + f₀. We use a quadratic limb-darkening law and for the limb-darkening coefficients we fit for the q₁, q₂ parameters from the parametrization of Kipping (2013). This parametrization ensures a physically realistic limb-darkening model for the star. The priors used for each parameter are provided in Table 1. These priors are selected to ensure that we fit physically realistic transit models to each light curve but also that we do not bias the results in any other way.

Table 1.

Priors used for each free parameter in the transit fitting detailed in Section 5. Note |$\mathcal {U}\left(A; B\right)$| denotes a uniform prior with lower bound A and upper bound B.

Parameter	Prior
Transit centre time, TC	|$\mathcal {U}(T_{\rm C, \ BLS}-0.1\, P_{\rm BLS}; T_{\rm C, \ BLS}+0.1\, P_{\rm BLS})$|
Orbital period, P	|$\mathcal {U}\left(0.95\, P_{\rm BLS}; 1.05\, P_{\rm BLS}\right)$|
Planet-to-star radius ratio, RP/R*	|$\mathcal {U}\left(0; 1\right)$|
Scaled semimajor axis, a/R*	|$\mathcal {U}\left(1.1; \infty \right)$|
Orbital inclination, i	|$\mathcal {U}\left(0^{\circ }; 90^{\circ }\right)$|
Limb-darkening parameter, q1	|$\mathcal {U}\left(0; 1\right)$|
Limb-darkening parameter, q2	|$\mathcal {U}\left(0; 1\right)$|
Constant flux offset, f0	|$\mathcal {U}\left(-0.05; 0.05\right)$|

For each object, 24 independent chains were each sampled for 7500 steps, following a burn-in phase of 2500 steps. This resulted in posterior distributions of 180  000 samples for each candidate. We use the results from our transit MCMC analysis to assess the likely planetary nature of our candidates.

As a first step, we consider the set of parameters for each object which resulted in the highest log likelihood values during the fitting; referred to in this paper as the ‘best-fitting’ parameters. We include only objects with a best-fitting planetary radius in the range |$0.6\, \mbox{$R_{\rm J}$}\le \mbox{$R_{\rm P}$}\le 2.0\, \mbox{$R_{\rm J}$}$| in our candidate list. Using the transit models computed using the best-fitting parameters we also calculate the ratio q_trans = T_dur/P, where T_dur is the transit duration. We accept only objects with 0.0 < q_trans ≤ 0.1 as likely planetary candidates. We test the threshold chosen for this criterion by calculating the q_trans value for all the planets we simulate for the Injection-Recovery tests (see Section 8). The highest q_trans value is 0.079, with 90 per cent of the simulated planets having q_trans < 0.03. Therefore, our chosen threshold is such that real transiting giant planets are unlikely to be excluded by this criterion. Using these two parameter threshold criteria, we identify just 93 of the 2102 candidates as being strong giant planet candidates.

For a final automated step, we use the posterior distributions from the MCMC sampling to assess the likely planetary nature of our candidates. We determine a posterior distribution for R_P using the R_P/R_* posterior and the host star radius from the TIC. Similarly, we use the posteriors for P and a/R_* along with Kepler’s 3rd law to determine a posterior distribution for the stellar density, ρ_*. We note that these ρ_* values are calculated assuming a circular orbit. With these distributions in R_P and ρ_* we then determine what fraction of the posterior falls in the giant planet regime, which we define with the following parameter ranges: |$0.6\, \mbox{$R_{\rm J}$}\le \mbox{$R_{\rm P}$}\le 2.0\, \mbox{$R_{\rm J}$}$| and 1.5 g cm⁻³ ≤ ρ_* ≤ 200 g cm⁻³. We then exclude from our candidate list any object with this fraction f_GP < 15 per cent. From these automated steps outlined above, we identified a list of 44 candidate giant exoplanets.

6 PLANET CANDIDATE MANUAL VETTING

We now perform a number of manual checks for each of the candidates found by the automated search in order to ensure we have a final list of strong candidates.

6.1 Visual inspection

For the first manual check, we visually inspected the TESS-SPOC 30-min cadence light curves for each of our 44 candidates. From this visual inspection, we identified 22 of our candidates which displayed clear signs of being astrophysical false positives or due to systematic effects. Examples of this included secondary eclipses and odd-even transit depth differences which had amplitudes below the automated thresholds and sharp ramps in flux at the end of TESS orbits which were folded into data gaps to mimic transit signals. These ramps are the result of increased scattered light at the end of the TESS orbit being imperfectly corrected for. These 22 systems are then removed from our candidate list.

6.2 Independently confirmed objects

Three of our candidates have had their natures confirmed prior to this work, which we independently detected with our pipeline. Two of these are the giant planets TOI-1130 c (TIC-254113311; 0.974 M_J; Huang et al. 2020b) and WASP-107 b (TIC-429302040; 0.12 M_J; Anderson et al. 2017). The third object is the known brown dwarf LP 261-75 b (TIC-67646988; 68 M_J; Irwin et al. 2018). Therefore, we exclude TIC-67646988 from our candidate list.

6.3 Blend scenario checks

Due to the large pixel scale of TESS (21 arcsec pix⁻¹), we must ensure that the transit signals we observe are not due to a nearby star blending into the target aperture. Using the target pixel files (TPFs) from the TESS-SPOC pipeline, it is possible to inspect each candidate to search for clear signs of a blend scenario. We download the TPFs for each candidate using the lightkurve python software (Lightkurve Collaboration 2018). We then perform two tests to search for evidence of the source of the signal being a neighbouring star and not the target.

The first of these involves generating light curves using different aperture sizes. We compare three different aperture sizes. The first is the aperture used by the TESS-SPOC pipeline. The second is a small aperture of just the single pixel at the location of the target star, and the third is a large 5 × 5 pixel aperture also centred at the target star location. By comparing the depths of these three light curves, we can reveal blend scenarios. These blend scenarios will in general be cases where a larger aperture results in a deeper transit, signifying that the signal is arising from a nearby star that is only partially within the target aperture, but more fully within the larger aperture. If all three light curves have equal depths then it is unlikely the signal arises from a nearby star. Similarly, a shallower transit for the larger aperture also points towards the signal arising from the target star. This is because the result of increasing the aperture size here is simply the inclusion of an increased amount of dilution, and so the signal being on any of the nearby stars is unlikely. The second test we performed is to generate light curves for each individual pixel in a 7 × 7 pixel grid centred on the target star location. With these light curves we can investigate whether the transit signals are clustered around the target star, or are offset.

From these tests, one of our candidates – TIC-174440134 – displayed clear signs of being a nearby eclipsing binary. The depth of the transit events was clearly correlated with the size of the photometric aperture used, and the individual pixel light curves showed a clear offset of the centre of the event from the target star. This object was therefore removed from our candidate list.

6.4 TESS higher cadence light curves

For some of our remaining candidates, we have access to higher cadence TESS photometry, which can often provide more information on the nature of a candidate.

For five of our candidates this higher cadence photometry is in the form of a TESS 2-min cadence SPOC light curve, from which we determined the following about these candidates. The 2-min cadence light curve for TIC-156067195 revealed a subtle difference in both the depth and the phase position of the odd and even transit events, identifying this candidate as an eclipsing stellar binary system with a slightly eccentric orbit. For TIC-190885165 the 2-min cadence photometry reveals more clearly the shape of the transit events. By refitting this photometry, we find a best-fitting radius of 2.59 R_J for the companion revealing its nature as most likely a stellar companion. TIC-389900760 and TIC-38460940 were initially considered as candidates because from the 30-min cadence alone they could have been giant planets on grazing orbits. The 2-min cadence photometry shows that they are not grazing, and so are most likely real planets but with |$\mbox{$R_{\rm P}$}\lt 0.6\, \mbox{$R_{\rm J}$}$|⁠, and so smaller than the giant planets we are searching for in this work. This is consistent with them being identified as TOI-2120.01 (⁠|$\mbox{$R_{\rm P}$}= 2.58\pm 0.24$| R_⊕) and TOI-805.01 (⁠|$\mbox{$R_{\rm P}$}= 1.28\pm 0.62$| R_⊕). Based on this, we exclude them both from our giant planet candidate list. For TIC-178709444, the 2-min cadence light curve contains no evidence that the candidate is a false positive. In fact the higher cadence confirms the transit as non-grazing (b < 1) providing us with increased confidence that the transit signal is produced by a planet and is not the result of an eclipsing binary.

For eight more of our candidates – TIC-46432937, TIC-67512645, TIC-95112238, TIC-165227846, TIC-202468443, TIC-243641947, TIC-271321097, and TIC-335590096 – there are TESS-SPOC 10-min cadence FFI light curves from the extended mission available. Using these data, we can identify TIC-271321097 as an eclipsing binary system. The extended mission light curve clearly displays primary and secondary eclipses. The secondary eclipses are what is seen in the primary mission 30-min cadence light curve, with the primary eclipses falling into the data gaps. TIC-271321097 is therefore excluded from our candidate list. For TIC-67512645, the 10-min cadence photometry revealed the transit as being flat-bottomed, whereas in the 30-min cadence photometry there was the possibility of a grazing eclipse, thereby strengthening the likelihood of this candidate being a real planet. Similarly for TIC-243641947 and TIC-335590096, the 10-min cadence photometry confirmed the transits as flat-bottomed, further supporting the nature of these two objects as planetary candidates. The 10-min cadence photometry for the remaining four objects contains no evidence that they are false positives, and so they remain as good transiting giant planet system candidates.

7 GIANT PLANET CANDIDATES

From our full planet search we have identified a final selection of 15 giant planet candidates. The numbers of objects from our sample at each step of the planet search pipeline are summarized in Fig. 6. The planetary parameters derived for our candidates in Section 5 are provided in Table 2 and we display the transit events for each candidate in Fig. 7. We provide details on the host stars of our candidates in Tables A1 and A2. Two of our giant planet candidates are already confirmed as giant planets and four have been identified as TOIs, for which we provide details in Table 3. A further two candidates have been independently reported as CTOIs; these are TIC-95112238 (Montalto et al. 2020) and TIC-77490011 (Montalto 2023). The remaining seven of our giant planet candidates are new candidates that were not previously known as planets or identified as candidates prior to our search.

Figure 6.

Flow chart of the steps of the planet search pipeline. The numbers indicate the number of stars from our sample for each step.

Figure 7.

Light curves for the 15 transiting giant planet candidates phase folded using the ephemerides from the fitting detailed in Section 5 and zoomed in to show the transit features. The orange lines show the best-fitting models. The TIC IDs and orbital periods are given for each system.

Table 2.

Planetary parameters for our giant planet candidates derived from the transit fitting in Section 5. The quoted values are the 50th percentiles of the posterior distributions, and the uncertainties are defined by the 16th and 84th percentiles and represent the 1σ uncertainty.

	TIC-95112238	TIC-311555090	TIC-46432937
TC (BJD – 2457000)	1468.539538 ± 0.000541	1817.344136 ± 0.001122	1468.630246 ± 0.000325
Period (d)	1.362066 ± 0.000060	1.416762 ± 0.000112	1.440455 ± 0.000037
RP (RJ)	|$0.79^{0.03}_{-0.04}$|	|$1.33^{0.61}_{-0.34}$|	|$1.57^{1.23}_{-0.45}$|
RP/R*	|$0.14^{0.01}_{-0.01}$|	|$0.58^{0.27}_{-0.15}$|	|$0.30^{0.23}_{-0.09}$|
ρ* (g cm−3)	|$17.03^{4.61}_{-5.24}$|	|$30.54^{8.95}_{-6.54}$|	|$5.45^{1.38}_{-0.46}$|
a/R*	|$11.9^{1.0}_{-1.4}$|	|$14.8^{1.3}_{-1.1}$|	|$8.4^{0.7}_{-0.2}$|
i (deg)	|$88.1^{1.4}_{-1.5}$|	|$87.0^{1.8}_{-1.1}$|	|$83.5^{2.1}_{-1.8}$|
	TIC-335590096	TIC-406717909	TIC-202468443
TC (BJD – 2457000)	1571.796562 ± 0.000817	1684.513652 ± 0.001532	1712.155769 ± 0.002109
Period (d)	1.522755 ± 0.000095	1.877180 ± 0.000214	2.068584 ± 0.000027
RP (RJ)	|$0.76^{0.14}_{-0.04}$|	|$1.55^{1.38}_{-0.33}$|	|$1.38^{1.09}_{-0.80}$|
RP/R*	|$0.22^{0.04}_{-0.01}$|	|$0.37^{0.33}_{-0.08}$|	|$0.46^{0.37}_{-0.27}$|
ρ* (g cm−3)	|$8.37^{2.82}_{-4.36}$|	|$13.64^{6.59}_{-3.55}$|	|$6.11^{4.12}_{-2.03}$|
a/R*	|$10.1^{1.0}_{-2.2}$|	|$13.6^{1.9}_{-1.3}$|	|$11.1^{2.1}_{-1.4}$|
i (deg)	|$88.0^{1.4}_{-3.7}$|	|$86.6^{2.5}_{-1.8}$|	|$83.6^{2.8}_{-1.8}$|
	TIC-165227846	TIC-60910638	TIC-73692250
TC (BJD – 2457000)	1572.424716 ± 0.000288	1764.926115 ± 0.001158	1791.488380 ± 0.001541
Period (d)	2.096696 ± 0.000041	2.149809 ± 0.000176	2.572953 ± 0.000300
RP (RJ)	|$1.14^{0.60}_{-0.10}$|	|$0.94^{1.29}_{-0.17}$|	|$0.76^{1.40}_{-0.06}$|
RP/R*	|$0.36^{0.19}_{-0.03}$|	|$0.21^{0.28}_{-0.04}$|	|$0.23^{0.42}_{-0.02}$|
ρ* (g cm−3)	|$12.13^{1.86}_{-1.21}$|	|$5.67^{10.35}_{-1.49}$|	|$21.11^{19.55}_{-8.58}$|
a/R*	|$14.1^{0.7}_{-0.5}$|	|$11.1^{4.6}_{-1.1}$|	|$19.5^{4.8}_{-3.1}$|
i (deg)	|$87.4^{0.8}_{-1.4}$|	|$85.6^{3.0}_{-2.5}$|	|$88.2^{1.3}_{-2.6}$|
	TIC-243641947	TIC-77490011	TIC-67512645
TC (BJD – 2457000)	1603.695100 ± 0.000444	1791.395003 ± 0.004345	1901.057603 ± 0.001274
Period (d)	2.592465 ± 0.000109	2.687991 ± 0.000772	3.370806 ± 0.000296
RP (RJ)	|$1.07^{0.03}_{-0.03}$|	|$0.86^{0.98}_{-0.52}$|	|$1.33^{0.57}_{-0.43}$|
RP/R*	|$0.29^{0.01}_{-0.01}$|	|$0.38^{0.43}_{-0.23}$|	|$0.61^{0.26}_{-0.20}$|
ρ* (g cm−3)	|$16.55^{2.01}_{-2.90}$|	|$17.15^{78.12}_{-11.84}$|	|$20.47^{5.18}_{-4.02}$|
a/R*	|$18.0^{0.7}_{-1.1}$|	|$18.7^{14.4}_{-6.1}$|	|$23.1^{1.8}_{-1.6}$|
i (deg)	|$89.3^{0.5}_{-0.7}$|	|$86.7^{2.5}_{-2.8}$|	|$87.7^{1.1}_{-0.6}$|
	TIC-178709444	TIC-429302040	TIC-254113311
TC (BJD – 2457000)	1572.678689 ± 0.001707	1574.147012 ± 0.000308	1657.903824 ± 0.000449
Period (d)	3.471712 ± 0.000474	5.721615 ± 0.000157	8.351069 ± 0.000350
RP (RJ)	|$1.51^{1.83}_{-0.79}$|	|$1.04^{0.02}_{-0.02}$|	|$1.95^{2.59}_{-0.89}$|
RP/R*	|$0.36^{0.44}_{-0.19}$|	|$0.15^{0.00}_{-0.00}$|	|$0.27^{0.36}_{-0.12}$|
ρ* (g cm−3)	|$7.57^{10.32}_{-2.14}$|	|$3.31^{0.37}_{-0.45}$|	|$2.96^{0.75}_{-0.39}$|
a/R*	|$16.9^{5.6}_{-1.8}$|	|$17.9^{0.6}_{-0.9}$|	|$22.2^{1.7}_{-1.0}$|
i (deg)	|$86.1^{2.9}_{-1.4}$|	|$89.0^{0.6}_{-0.5}$|	|$87.1^{0.8}_{-0.9}$|

Table 3.

We provide a summary of the spectroscopic follow-up data available to date for our giant planet candidates, as well as providing any comments on their natures for the solved systems or the triceratops FFPs for the currently un-dispositioned candidates.

TIC	Radial velocites	Comments
95 112 238	NO	CTOI (Montalto et al. 2020)
		triceratops FPP = 0.24
311 555 090	NO	triceratops FPP = 0.62
46 432 937	NO	triceratops FPP = 0.025
335 590 096	NO	TOI-4860
		triceratops FPP = 1.00
406 717 909	NO	triceratops FPP = 0.67
202 468 443	NO	TOI-5268
		triceratops FPP = 0.94
165 227 846	NO	triceratops FPP = 1.00
60 910 638	ESPRESSO	ESPRESSO observations show this
		to be an eclipsing binary
73 692 250	NO	triceratops FPP = 0.35
243 641 947	ESPRESSO	TOI-3235
		Radial velocity observations confirm as
		a giant planet (Hobson et al. 2023)
77 490 011	NO	CTOI (Montalto 2023)
		triceratops FPP = 0.22
67 512 645	ESPRESSO	Radial velocity observations to date
		consistent with a giant planet
178709444	ESPRESSO	TOI-762
		Radial velocity observations to date
		consistent with a giant planet
429 302 040	Euler/CORALIE	WASP-107 (Anderson et al. 2017)
254 113 311	CHIRON	TOI-1130 (Huang et al. 2020b)

We also compare our sample of candidates to the full population of known transiting giant planets in Fig. 8 and to the subset of known transiting giant planets with low-mass host stars in Fig. 9. From these two figures, it is evident that our search has extended the population of known transiting gas planets to lower stellar mass hosts than ever before. From Fig. 8 we can see that our giant planet candidates have host stars spanning almost the full stellar mass range of our sample. We have detected eight candidates whose host stars have masses |$\mbox{$M_{*}$}\lt 0.4\, \mbox{${\rm M}_{\odot }$}$|⁠, with some as low mass as |$\mbox{$M_{*}$}\approx 0.2\, \mbox{${\rm M}_{\odot }$}$|⁠. If confirmed as real giant planets, these candidates will provide strong tension with our current understanding of planet formation (e.g. Pascucci et al. 2016; Burn et al. 2021).

Figure 8.

Radii of our giant planet candidates compared to the known transiting exoplanets plotted in Fig. 1 as a function of the mass of the host star. The magenta stars show the 15 gas giant planet candidates detected in this study, with the filled stars denoting the five confirmed giant planets. For these candidates, the R_P values plotted are the best-fitting values from the fitting detailed in Section 5.

Figure 9.

Comparison of our 15 giant planet candidates (magenta stars) to the Gaia colour–magnitude distribution of the 91 306 low-mass stars in our sample (density heat map). We also plot the 12 known transiting gas giant planets with hosts stars with masses |$\le 0.65\, \mbox{${\rm M}_{\odot }$}$| as the blue squares for comparison to our candidates.

From the colour–magnitude diagram in Fig. 9, we can see that some of our targets appear at slightly brighter absolute magnitudes than the bulk population, for a given Gaia colour. Such an effect could be indicative of these objects being binary stars, with the additional star resulting in the objects appearing much brighter than would be expected for a single star of the same spectral type. However, we also note that one of the known giant planets, HATS-74A b (Jordán et al. 2022), also sits raised above the line. This system consists of a low-mass star that hosts a transiting giant planet, as well as a bound stellar companion on a wide orbit. Such wide stellar companions have been found to be common for hot Jupiter host stars (Knutson et al. 2014; Ngo et al. 2016). As such, we cannot exclude these objects as giant planet candidates simply based on their position on the colour–magnitude diagram.

The Gaia astrometric re-normalized unit weight error (RUWE; Lindegren 2018) can also be used to indicate the potential presence of a stellar binary, with a values of RUWE > 1.4 often used to signify the presence of a stellar companion in the system (Lindegren et al. 2021a,b). For our candidates we identify two with RUWE > 1.4. The first of these is TIC-429302040 with RUWE = 1.54, which is the known exoplanet system WASP-107 (Anderson et al. 2017). There is a known outer planetary companion in the system (WASP-107 c; 1088 d; |$\mbox{$M_{\rm P}$}\, \sin i = 0.36\, \mbox{$M_{\rm J}$}$|⁠; Piaulet et al. 2021) although there is no known stellar companion in the system. The second is TIC-311555090 with RUWE = 2.29. Similarly with above, this object cannot be excluded as a giant planet candidate based solely on this RUWE value, but we report it here for reference.

7.1 ESPRESSO radial velocity monitoring

We note that not all of candidates may be true giant planets. In order to uncover their true natures, we have begun a program to obtain mass measurements for our candidates using the ESPRESSO spectrograph (Pepe et al. 2020) on the VLT. The current status of the spectroscopic follow-up of our candidates is summarized in Table 3. ESPRESSO has been successfully used over recent years to confirm similar systems (eg. HATS-71 b; Bakos et al. 2020). We have been reducing the ESPRESSO data throughout the program using the publicly available pipeline, which runs in the ESOReflex environment (Freudling et al. 2013).

The full spectroscopic follow-up for all of our candidates is beyond the scope of this work, and is likely a couple of years at least away, however, here we report some results to date that are consequential for this study. To date we have obtained spectroscopic observations for four of our candidates: TIC-243641947 (TOI-3235), TIC-178709444 (TOI-762), TIC-67512645, and TIC-60910638. For TIC-243641947, we have a completed spectroscopic orbit confirming the transiting companion as a giant planet. The full analysis of this system is outside the scope of this paper and is the subject of a recent paper (Hobson et al. 2023). For TIC-178709444 and TIC-67512645, the spectroscopic data we have obtained to date also shows the transiting companions to be giant planets, although a small amount of further spectroscopic observations are required to fully confirm both. For the occurrence rate analysis in this paper, we therefore take these three candidates to be genuine giant planets. For TIC-60910638, the spectroscopic observations reveal its true nature to be an eclipsing binary system and therefore exclude it from our occurrence rate analysis, although we include it in our reported candidate list for completion.

8 INJECTION AND RECOVERY TESTS

In order to derive occurrence rates for giant planets we must determine the detection efficiency of our search pipeline. To do this we performed injection-recovery simulations, which have been used successfully for studies into Kepler results (e.g. Christiansen et al. 2020).

In order to accurately replicate the noise properties and observing windows for the light curves we used real TESS-SPOC 30 min light curves, choosing as our base light curve set the selection of our input sample which did not yield a significant BLS detection. We then injected the transiting planet signals, simulating 10 different planets for each star. The parameters for each simulated planet were randomly selected from the uniform distributions. The parameters and corresponding distributions are as follows: orbital period, P – selected from a range of 1–10 d, transit centre time, T_C – selected such that the first transit lies with one orbital period of the start of the light curve, planet radius, R_P – selected from a range of 0.6–2.0 R_J, and impact parameter, b – selected from a range of 0 − 1 + R_P/R_*.

The stellar parameters, as obtained from the TIC, were then used along with the randomly drawn parameters to derive the required parameters to generate the batman light-curve model. The stellar radius, R_*, and R_P were used to derive the radius ratio, R_PR_*. The stellar mass, M_*, along with R_* and P were used to derive the scaled semimajor axis of the orbit, a/R_*using Kepler’s Third Law, assuming circular orbits. Finally, a/R_* and b were used to calculate the orbital inclination, i. Note that we assume circular orbits for all simulated planets.

Using this method, we simulated 873 760 transiting planets, which we then passed through our planet search pipeline, as if they were real light curves. Every simulated light curve was searched using BLS, and the simulated systems which yielded a significant detection at a period within 5 per cent of the injected period were then passed to the vetting steps we outline in Section 4. Our pipeline yielded 780 379 significant BLS detections from our simulated planet population, of which a total of 759 632 simulated planets then successfully passed the vetting checks described in Section 4. The transit fitting (Section 5) is a key step in the planet search pipeline, however we lack the computing resources to perform a transit fit on all 759 632 simulated planets. Instead, we selected at random a single simulated planet for each star and performed a transit fit for these. By doing this, we were able to obtain a representative estimate of the effect had by the transit fitting on the overall detection efficiency. We do not perform any visual inspection or further manual vetting (see Section 6) on our simulated planets. Therefore the false identification of any real planets as false positives due to this manual vetting is not captured in these detection efficiencies. Given the characteristically strong signals of transiting hot Jupiters compared with the photometric noise and any systematics we do not expect this to significantly affect our results.

The detection efficiencies we calculate for our planet search pipeline are presented in Fig. 10. We can see that across the majority of the {P; R_P} parameter space we study in this search we have high detection efficiencies, and that our sensitivity to transiting planets decreases for longer orbital periods. This is to be expected, as the signal-to-noise ratio of the phase-folded transit event decreases for these planets due to the smaller number of total transit events in the light curve.

Figure 10.

2D map showing the detection efficiency of our planet search pipeline and how this varies for different planet radii and orbital periods. These sensitivities take into account the transit search, light-curve vetting, and transit fitting analysis steps of the pipeline. The colour of each cell gives the detection efficiency withing that cell. The numbers printed onto each cell also give these detection efficiencies in per cent. The values here highlight the sensitivity of our pipeline to transiting giant planets.

In the large majority of cases where a simulated planet would be excluded from our giant planet candidate list by the transit fitting analysis this is due to the fitting algorithms finding an inaccurate impact parameter, b, for these planets, leading to non-grazing transits being fit with grazing transit models and vice-versa. The 30-min cadence of the TESS FFI observations reduces our knowledge of the true transit shape causing these inaccurate fits, resulting in either an over- or underestimation of the planet radius. This has more of an impact on simulated planets with longer orbital periods, as they exhibit fewer transit events, which combined with the 30-min cadence of the TESS FFI observations results in a sparse sampling of the transit event in the phase-folded light curve. This sparse sampling then amplifies the lack of knowledge available for the true transit shape, and therefore the true planet radius. Similarly, simulated planets with radii close to our giant planet radius limits of 0.6 R_J and 2.0 R_J are more susceptible to be misidentified as false positives due to this miscalculation of the impact parameter, as a smaller inaccuracy in the recovered impact parameter, and thus planet radius, is needed to result in best-fitting parameters for these simulated planets outside the giant planet regime.

From the overall detection efficiency of our pipeline, which is plotted in Fig. 11, it is clear that there is a strong dependence with orbital period. This arises as the geometric transit probability itself depends strongly on the orbital separation between the star and planet. The detection efficiency of the BLS search and transit fitting analysis also depends on the orbital period of the planet. It is therefore not surprising that the majority of our giant planet candidates are found with orbital periods P ≤ 3.5 d.

Figure 11.

2D colour map displaying the overall detection efficiency of the planet search pipeline, split across multiple bins for the injected orbital period and planet radius of the simulated planet. These detection efficiencies take into account all of the pipeline steps as well as the geometric transit probability. The numbers in each bin give the detection efficiency for the corresponding bin in per cent. The stars highlight the locations in {P; R_P} parameter space of our fifteen giant planet candidates.

9 GIANT PLANET OCCURRENCE RATES

Using the candidates from our planet search and results from our planet injection-recovery simulations, we can now place constraints on the occurrence rates of giant planets orbiting low-mass stars. We can calculate an occurrence rate, f_occ, using the following equation

$$\begin{eqnarray} f_{\rm occ} = \frac{n_{\rm pl}}{N_{\rm pr}}, \end{eqnarray}$$

(4)

where n_pl is the number of planets detected and N_pr is the number of stars amenable to the discovery of an exoplanet (see e.g. Bayliss & Sackett 2011; Gan et al. 2022). This factor is calculated as follows:

$$\begin{eqnarray} N_{\rm pr} = N_{\ast } \frac{\Sigma _{i=1}^{N_{\rm sim}} \delta _{{\rm vet}, i} \mathcal {P}_{{\rm tr}, i}}{N_{\rm sim}} \frac{\Sigma _{j=1}^{N_{\rm fit}} \delta _{{\rm fit}, j}}{N_{\rm fit}}, \end{eqnarray}$$

(5)

where N_* is the number of stars in the stellar sample, N_sim is the total number of planets simulated, N_fit is the number of simulated planets which were passed to the transit fitting, δ_vet,i is a detection delta function that equals one if the simulated planet was detected by BLS and passed the vetting checks, else equals zero, δ_fit,j is a similar delta function dependent on whether a simulated planet was identified as a giant planet candidate by the transit fitting, and |$\mathcal {P}_{{\rm tr}, i}$| is the geometric transit probability of the simulated planet. This transit probability is defined as

$$\begin{eqnarray} \mathcal {P}_{{\rm tr}, i} = \frac{R_{\rm P} + R_{\ast }}{a} . \end{eqnarray}$$

(6)

It is clear from Fig. 11 that both the detection efficiency and planet occurrence rates will be functions of the orbital period, P, and planet radius R_P. Therefore, to calculate the overall giant planet occurrence rate we first calculate individual f_occ values for each of the cells in {P; R_P} parameter space in Fig. 11, and then calculate the sum of these individual values. Similarly to estimate the uncertainty on the occurrence rate, |$\sigma _{f_{\rm occ}}$|⁠, for each cell in {P; R_P} parameter space in which at least one planet candidate is detected we calculate the uncertainties on n_pl and N_pr as |$\sqrt{n_{\rm pl}}$| and |$\sqrt{N_{\rm pr}}$|⁠, respectively, as estimated from Poisson counting statistics. The individual cell occurrence rate uncertainty is then computed for each cell as

$$\begin{eqnarray} \left(\frac{\sigma _{f_{\rm occ}}\left(P; R_{\rm P}\right)}{f_{\rm occ}\left(P; R_{\rm P}\right)}\right)^2 = \left(\frac{\sigma _{n_{\rm rpl}}\left(P; R_{\rm P}\right)}{n_{\rm rpl}\left(P; R_{\rm P}\right)}\right)^2 + \left(\frac{\sigma _{N_{\rm pr}}\left(P; R_{\rm P}\right)}{N_{\rm pr}\left(P; R_{\rm P}\right)}\right)^2 . \end{eqnarray}$$

(7)

Due to the large number of simulated planets injected in this work the first term dominates the occurrence rate uncertainty. The individual values for each cell are then summed in quadrature to calculate the overall occurrence rate uncertainty.

For each {P; R_P} bin, the value used for n_pl is calculated as

$$\begin{eqnarray} n_{\rm pl}\left(P; R_{\rm P}\right) = \Sigma _{i=1}^{N_{\rm cands}} 1-{\rm FPP}_i , \end{eqnarray}$$

(8)

where N_cands is the number of candidates in the parameter space range and FPP is a false positive probability that we assign to each candidate. The two previously known planets – TIC-429302040 and TIC-254113311 – as well as the three planets we have spectroscopically confirmed – TIC-243641947, TIC–178709444, and TIC-67512645 – are all assigned FPP = 0. TIC-60910638, which our spectroscopic follow-up reveals is an eclipsing binary, is assigned FPP = 1. For each of the remaining candidates, we calculate an FPP value using triceratops (Giacalone & Dressing 2020; Giacalone et al. 2021). triceratops estimates the relative probabilities that a given transit signal was produced by a transiting planet or a number of false positive scenarios using a Bayesian analysis. These false positive scenarios include nearby, background, and bound but unresolved eclipsing binaries. For the nine remaining candidates, we obtain a mean FPP of 0.56, and we report the FPP values for our candidates in Table 3. For our full low-mass star sample (⁠|$0.088\mbox{${\rm M}_{\odot }$}\le \mbox{$M_{*}$}\le 0.71\mbox{${\rm M}_{\odot }$}$|⁠), we calculate an occurrence rate of close-in giant planets of 0.194 ± 0.072 per cent.

It is important to note that one must be cautious when using such validation techniques as we use here for giant planets. This is due to the fact that the radii of such planets are very similar to brown dwarfs and very low mass stars (e.g. Mayo et al. 2018). This makes the models used by these validation techniques for giant planets indistinguishable from those for brown dwarfs and low-mass stars. triceratops also includes a prior on planet radius which penalizes giant planets with low-mass host stars. These two effects both impact the reliability of the numerical FPP estimates for each candidate, and as such they cannot be used to validate or reject a given individual candidate. However without spectroscopic observations for all candidates they provide the best estimates available to us for the overall false positive probability of our candidates, and therefore they also provide the current best estimate for the occurrence rate of these giant planets. We note that these spectroscopic observations are underway for our candidates. To investigate the sensitivity of our occurrence rate results on the FPP estimates we consider the cases in which all or none of these candidates are real to set the upper and lower limits of the occurrence rates we calculate. For our full low-mass star sample, we calculate such upper and lower limits of 0.267 ± 0.079 and 0.139 ± 0.066 per cent. We emphasize here that even in the case in which all the currently un-dispositioned candidates are false positives, we have a significantly non-zero occurrence rate for close-in giant planet with low-mass host stars.

The large size of our sample and the stellar mass range it covers allows us to investigate the dependence of giant planet occurrence on stellar mass for our low-mass stars. We separate our stellar sample into three sub-samples with mass ranges defined in order to result in a similar number of stars for each sub-sample, identifying these three similar-size sub-samples as having mass limits of 0.088–0.26, 0.26–0.42, and 0.42–0.71 M_⊙. We calculate occurrence rates for each sub-sample finding values of 0.137 ± 0.097 per cent (0.088–0.26 M_⊙), 0.108 ± 0.083 per cent (0.26–0.42 M_⊙), and 0.29 ± 0.15  per cent (0.42–0.71 M_⊙). We plot these occurrence rates, along with results from previous studies, in Fig. 12, and occurrence rates we calculate in this work are summarized in Table 4. These occurrence rate measurements are related to Poisson distributions governing the expected number of detected planets. From these distributions, we can calculate 95 per cent lower limits for the occurrence rates in each mass range. Doing so we find lower limits of 0.048 per cent (0.088–0.26 M_⊙), 0.042 per cent (0.26–0.42 M_⊙), and 0.10 per cent (0.42–0.71 M_⊙). Each of these three mass ranges also includes at least one of the five confirmed giant planets in our sample, and so in the case where none of the remaining nine candidates are real planets, we calculate occurrence rates of 0.081, 0.059, and 0.24 per cent, with 95 per cent lower limits of 0.029, 0.021, and 0.085 per cent. Therefore, our results show that a population of close-in giant planets exists even for the lowest range of stellar masses we study.

Figure 12.

Dependence of the occurrence rate of giant planets with short orbital periods (⁠|$P \le 10\, d$|⁠) on the stellar mass of the host star. We compare our results derived in this work (filled magenta squares) to previous TESS results from Beleznay & Kunimoto (2022) (filled black circles) and Gan et al. (2022) (unfilled black circle). We also compare to the Kepler results from Petigura et al. (2018) (blue triangle). The x-axis errorbars give the span of the stellar mass ranges used to determine each occurrence rate value. The y-axis errorbars for our results give the standard deviation of the derived Poisson distributions underlying the number of planets we detect. We also plot the 95 per cent lower occurrence rate limits in the case in which none of the remaining nine candidates are real (see text in Section 9), shown by the magenta arrow heads. We note that we have at least one confirmed giant planet in each of our mass ranges, and so we are confident that the occurrence rates in all of the mass ranges are non-zero.

Table 4.

Summary of the short period giant planet occurrence rates we measure in this work, for different stellar mass ranges. The ‘Max’ labels denote the occurrence rates that have been calculated assuming all nine of the currently unconfirmed candidates are real, and the ‘Min’ labels denote the occurrence rates that have been calculated assuming none of these nine candidates are real.

Stellar mass range (M⊙)	Occurrence rate (per cent)
0.088–0.71	0.194 ± 0.072
0.088–0.71 (Max)	0.267 ± 0.079
0.088–0.71 (Min)	0.139 ± 0.066
0.088–0.4	0.134 ± 0.069
0.088–0.4 (Max)	0.249 ± 0.090
0.088–0.4 (Min)	0.074 ± 0.053
0.088–0.26	0.137 ± 0.097
0.26–0.42	0.108 ± 0.083
0.42–0.71	0.29 ± 0.15

The distribution of our candidates in Fig. 11 suggests the possibility of a dependence of the occurrence rate on the planet radius, R_P. To test whether this distribution is indicative of the underlying population, or simply an effect of the detection efficiency, we calculate the occurrence rate as a function of R_P. These results are plotted in Fig. 13 and we see that the occurrence rates show clear variation with R_P. We see a peak in occurrence rates for |$\mbox{$R_{\rm P}$}\sim 1 \mbox{$R_{\rm J}$}$| and occurrence rates consistent with zero for |$\mbox{$R_{\rm P}$}\ge 1.4 \mbox{$R_{\rm J}$}$|⁠. Whereas for more massive host stars – |$\mbox{$M_{*}$}\ge 1.2\, \mbox{${\rm M}_{\odot }$}$| – there exist transiting planets with radii |$\mbox{$R_{\rm P}$}\gt 1.5\, \mbox{$R_{\rm J}$}$| (see Fig. 8). It has been previously suggested that extreme levels of irradiation of the planet by the star, of the order of 10⁵–10⁶ W m⁻², are required to cause the inflation of the planetary atmosphere necessary to reach planetary radii |$\gtrsim 1.1 \mbox{$R_{\rm J}$}$| (Sestovic, Demory & Queloz 2018). A planet in a 3 d orbit around a 0.45 M_⊙ star will receive an irradiation of around 5.7 × 10⁴ W m⁻². Therefore, it is possible that the derived lack of high radius giant planets for low-mass stars is because these stars are unable to inflate the atmospheres of the planets orbiting them. Previous occurrence rate studies have shown an increase in giant planet occurrence rate for longer orbital periods (e.g. Petigura et al. 2018). The distribution of our candidates in Fig. 11 shows that the majority of our candidates have short periods |$\lt 4\, d$|⁠. However, the detection efficiencies are also much higher for short orbital periods. Plotting the occurrence rate as a function of orbital period in Fig. 14 we see no clear trend in the occurrence rates. At this stage, the numbers of planets are too low in order to discern a significant trend with orbital period.

Figure 13.

Dependence of our measured occurrence rate of giant planets with short orbital periods (⁠|$P \le 10\, d$|⁠) and low-mass host stars on planet radius. Dashed horizontal line represents the 90  per cent confidence level upper limit for the bin in which there are no planet candidates.

Figure 14.

Dependence of our measured occurrence rate of giant planets with short orbital periods (⁠|$P \le 10\, d$|⁠) and low-mass host stars on the orbital period of the planet.

9.1 Comparison with previous works

Many previous works have measured the occurrence rates of giant planets with varying stellar hosts. In particular, results from the Kepler survey – a sample of F, G, and K dwarf stars – yielded short period giant planet occurrence rates of 0.43 ± 0.05 per cent (Fressin et al. 2013, |$\mbox{$R_{\rm P}$}= 6 - 22 \mbox{$R_{\hbox{$\oplus $}}$}$|⁠; P = 0.8–10 d) and |$0.57^{+0.14}_{-0.12}$| per cent (Petigura et al. 2018, |$\mbox{$R_{\rm P}$}= 8 - 24 \mbox{$R_{\rm J}$}$|⁠; P = 1–10 d). The occurrence rates we derive here for low-mass stars are significantly lower (>1σ) than these Kepler results. This is the case even when we consider the upper limit scenario in which all of our candidates are real. More recently, Zhou et al. (2019) used results from the first seven TESS sectors to study the occurrence rate of hot Jupiters with host stars with varying masses, deriving values of 0.71 ± 0.31 per cent (⁠|$0.8-1.05\, \mbox{${\rm M}_{\odot }$}$|⁠), 0.43 ± 0.15 per cent (⁠|$1.05-1.4\, \mbox{${\rm M}_{\odot }$}$|⁠), and 0.26 ± 0.11 per cent (⁠|$1.4-2.3\, \mbox{${\rm M}_{\odot }$}$|⁠). Beleznay & Kunimoto (2022) built upon this early work using photometry from the first two years of the TESS mission to refine these giant planet occurrence rates to 0.55 ± 0.14 per cent (⁠|$0.8-1.05\, \mbox{${\rm M}_{\odot }$}$|⁠), 0.36 ± 0.06 per cent (⁠|$1.05-1.4\, \mbox{${\rm M}_{\odot }$}$|⁠), and 0.29 ± 0.05 per cent (⁠|$1.4-2.3\, \mbox{${\rm M}_{\odot }$}$|⁠). The first TESS giant planet occurrence rate measurement of 0.27 ± 0.09 per cent for M-dwarf host stars (0.45–0.65 M_⊙) was provided by Gan et al. (2022).

We compare our results to some of these previous studies in Fig. 12. The first clear result from this comparison is that the occurrence rates of short orbital period giant planets decreases with decreasing stellar mass for host stars with |$\mbox{${\rm M}_{\odot }$}\le 1\, \mbox{${\rm M}_{\odot }$}$|⁠. Considering our result from our full sample – 0.194 ± 0.072 per cent – we find it is less than the Petigura et al. (2018) result at a level of |$2.69\, \sigma$| and less than the Beleznay & Kunimoto (2022) |$0.8-1.05\, \mbox{${\rm M}_{\odot }$}$| result at a level of |$2.27\, \sigma$|⁠. Considering the case in which all nine of the remaining candidates are real planets, our upper limit occurrence rate of 0.267 ± 0.079 per cent is still less than both of these results, at a level of |$2.10\, \sigma$| for Petigura et al. (2018) and a level of |$1.76\, \sigma$| for Beleznay & Kunimoto (2022). Comparing to the Gan et al. (2022) results for early M-dwarfs we find that our derived result for the |$0.42-0.71\, \mbox{${\rm M}_{\odot }$}$| selection of our stellar sample – 0.29 ± 0.15 per cent – is fully consistent with their measured occurrence rate for a similar mass range. We note that we have used a different set of TESS FFI light curves than Gan et al. (2022); we have used those produced by the SPOC pipeline while they made use of data from the QLP pipeline (Huang et al. 2020a). This different set of light curves, along with the use of an independent planet search methodology, allows us to use our results for early M-dwarfs to corroborate the Gan et al. (2022) findings. The wider stellar mass range of our sample also allows us to measure the occurrence rates for giant planets with much lower mass host stars. If we consider the stars in our sample with |$\mbox{$M_{*}$}\le 0.4 \mbox{${\rm M}_{\odot }$}$| we find an occurrence rate of 0.134 ± 0.069 per cent, which is less than the Gan et al. (2022) result at a level of |$1.15\, \sigma$|⁠. This is some evidence, if still statistically marginal, that the occurrence rates of giant planets continues to decrease for mid- and late- M-dwarfs compared to early M-dwarfs. Comparing our results from the low and middle mass ranges of our sample to the upper mass range we also find evidence of this decrease for M-dwarfs with later spectral types at a similar statistical level. The opposite trend has been observed for small planets, which have been shown to be more common around M-dwarfs than more massive host stars (e.g. Dressing & Charbonneau 2015; Kunimoto & Matthews 2020; Pinamonti et al. 2022). This difference in occurrence rate variation with mass of the host star could imply a formation efficiency for planets around low-mass stars which favours the formation of multiple small planets over a single giant planet.

9.2 Implications for giant planet formation

The dependence of giant planet occurrence rates on the mass of the host star is a clear prediction of the core-accretion planet formation theory (e.g. Laughlin et al. 2004; Burn et al. 2021). This theory predicts a decrease in occurrence rates for |$\mbox{$M_{*}$}\lt 1\, \mbox{${\rm M}_{\odot }$}$| and comparing the results for our full sample to those for solar-like stars from Kepler (Petigura et al. 2018) andTESS (Beleznay & Kunimoto 2022) we indeed find such a decrease. This agreement in the observed trend with the theoretical prediction suggests that the core-accretion mechanism dominates the formation of these close-in giant planets. The Burn et al. (2021) planet population synthesis, which uses the core-accretion formation theory, also predicts that the formation of giant planets becomes impossible for stars with |$\mbox{$M_{*}$}\le 0.4\, \mbox{${\rm M}_{\odot }$}$|⁠. Studies of protoplanetary discs have shown their mass to decrease for lower masses of the central star (e.g. Andrews et al. 2013; Kurtovic et al. 2021) and therefore for these low-mass stars the predictions are that the disc mass is not sufficient to form a giant planet (Pascucci et al. 2016). However, our results show that close-in giant planets can and do exist for these low-mass host stars. From our study, we derived a giant planet occurrence rate for host stars with |$\mbox{$M_{*}$}\le 0.4\, \mbox{${\rm M}_{\odot }$}$| of 0.134 ± 0.069 per cent, with a 95 per cent lower limit of 0.05 per cent. In particular, two of our giant planet candidates for which we have obtained spectroscopic follow-up confirmation – TIC-243641947 and TIC-67512645 – both have host stars with |$\mbox{$M_{*}$}\lt 0.4 \mbox{${\rm M}_{\odot }$}$|⁠. So for the case in which we consider that none of our remaining nine giant planet candidates are real, we derive an occurrence rate of 0.074 per cent with a 95 per cent lower limit of 0.026 per cent.

There must be some pathway through which these stars can form giant planets. First, we consider that the current main inhibitor for the predicted formation of these planets is the limits on the mass of the protoplanetary disc. Therefore, if these low-mass stars were capable of supporting much more massive disks than currently expected, even rarely, this would allow for the formation of these giant planets. Observational studies of disc-hosting low-mass stars could allow the occurrence of such massive discs to be constrained. It is important here to note that the dust masses for protoplanetary discs are currently derived from the flux received at millimeter wavelengths, assuming the dust continuum emission is optically thin (e.g. Hildebrand 1983; Kurtovic et al. 2021). If this emission is in fact optically thick, or thicker than currently assumed, then it is possible that we are underestimating the masses of the observed protoplanetary discs. We also note that to date protoplanetary discs have not been observed at a very young age, as at these ages they are still embedded in gas clouds and hard to observe. As such, these discs could begin with a much higher mass than when we can observe them, and so it is possible that their initial masses are sufficient to support giant planet formation.

Another possibility is that these planets did not form through core-accretion, but instead formed by gravitational instability (Boss 1997). It has been shown that giant planets are capable of forming through this mechanism for host stars with masses as low as 0.1 M_⊙ (Boss 2006; Mercer & Stamatellos 2020). Mercer & Stamatellos (2020) show that the giant planets which form around low-mass stars through gravitational instability have high masses |$\gtrsim 2\, \mbox{$M_{\rm J}$}$|⁠. Spectroscopic follow-up of our candidates is required to measure their masses determine whether they are massive planets that could have formed through gravitational instability.

10 CONCLUSIONS AND FUTURE OUTLOOK

We have presented a systematic search through TESS 30-min cadence light curves for transiting gas giant planets orbiting low-mass stars. These systems are predicted to be rare by planet formation theory, and so the aim of this search is to derive a robust occurrence rate for these systems for the first time. We have presented our planet search and candidate vetting pipeline, using which we identified fifteen giant planet candidates from an initial sample of 91 306 low-mass stars. Of these candidates, two are previously known transiting planets (Anderson et al. 2017; Huang et al. 2020b) and through our own preliminary spectroscopic follow-up we have shown three more of our candidates to be giant planets, and one to be an eclipsing binary. Spectroscopic observations over the coming years will be required to uncover the nature of the remaining nine candidates. In addition, seven of our candidates had not been identified as planet candidates prior to this study, including TIC-67512645. Our candidate giant planets have stellar hosts which are in general lower mass than the hosts of currently known transiting giant planets (see Fig. 8). Therefore, this study is providing strong evidence that the population of transiting giant planets extends to lower stellar mass host stars than previously expected.

We have also performed planet injection-recovery simulations to estimate the detection efficiency of our pipeline. Using these results and our giant planet candidate list we have constrained the occurrence rate of giant planets around low-mass stars, deriving an occurrence rate of 0.194 ± 0.072 per cent for our low-mass star sample. We also derive occurrence rates for three separate host star mass ranges, calculating values of 0.137 ± 0.097 per cent (0.088–0.26 M_⊙), 0.108 ± 0.083 per cent (0.26–0.42 M_⊙), and 0.29 ± 0.15 per cent (0.46–0.71 M_⊙), extending our knowledge of giant planet populations to lower stellar mass hosts than previously studied. The occurrence rate we calculate for our highest mass bin is fully consistent with the results from an independent study focusing on host stars with masses in the range 0.45–0.65 M_⊙ (Gan et al. 2022). Comparing with occurrence rate studies for higher mass host stars (e.g. Fressin et al. 2013; Petigura et al. 2018; Beleznay & Kunimoto 2022), we demonstrate that giant planets are less common around low-mass stars than solar-type stars, as predicted by the core-accretion planet formation theory (e.g. Laughlin et al. 2004; Burn et al. 2021). Our results provide some strong early evidence that giant planets are even less common around late M-dwarfs than early M-dwarfs but that giant planets can exist with host stars as low mass as |$0.2-0.3 \mbox{${\rm M}_{\odot }$}$|⁠. It has previously been asserted that lower protoplanetary disk masses (e.g. Pascucci et al. 2016; Kurtovic et al. 2021) and longer Keplerian timescales (e.g. Laughlin et al. 2004) completely inhibits the formation of giant planets around stars with masses as low as this. Therefore, while our results for the higher mass stars in our sample are consistent with core-accretion, our results for the lower mass stars in our sample – |$\mbox{$M_{*}$}\le 0.4\, \mbox{${\rm M}_{\odot }$}$| – present some conflict with the current understanding of how giant planets form. Further observations, both to measure the masses of our candidates and study the masses of protoplanetary discs around low-mass stars, will help determine how these planets formed.

The constraints we have placed on giant planet occurrence are currently limited by the fact that two thirds of our candidates are unconfirmed, and they will be significantly strengthened by the confirmation of these candidates. The work to obtain these confirmations is underway and over the next few years we hope to be able to improve our constraints on the occurrence rates as the follow-up effort continues. We also provide full details on our candidates so that the wider community is able to assist in this effort. We note that some of the nine as yet unconfirmed candidates have host stars with |$\mbox{$M_{*}$}\lt 0.4 \mbox{${\rm M}_{\odot }$}$|⁠. Therefore the new generation of stabilized spectrographs operating at (near-)infrared wavelengths, such as NIRPS (Bouchy et al. 2017) or SPIRou (Thibault et al. 2012), can play a major role in refining the occurrence rates of these systems.

During the fitting analysis outlined in Section 5 we realized the limitations of the long cadence (30 mi) observations in characterizing the true shape of the transit signal. These transit shapes are incredibly important for giant planets transiting small radius stars because there is a much higher probability of the planet crossing the limb. Moving the FFIs to higher cadence – 10 min in the first extended mission and then 200 s in the second extended mission⁵ – will greatly improve our ability to classify giant planets transiting low-mass stars.

Acknowledgement

The contributions at the Mullard Space Science Laboratory by EMB and VVE have been supported by STFC through the consolidated grant ST/W001136/1. We thank Paola Pinilla for useful discussions regarding the planet formation interpretation of our results discussed in Section 9.2. This paper includes data collected with the TESS mission obtained from the MAST data archive at the Space Telescope Science Institute (STScI). Funding for the TESSmission is provided by the NASA Explorer Program. STScI is operated by the Association of Universities for Research in Astronomy, Ibc., under NASA contract NAS 526555. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

DATA AVAILABILITY

This work made use of the publicly available TESS-SPOC 30 min cadence FFI light curves (Caldwell et al. 2020). These are available as a High Level Science Product from the Mikulski Archive for Space Telescopes here: https://archive.stsci.edu/hlsp/tess-spoc.

Footnotes

1

Data accessed from NASA Exoplanet Archive on 2022 November 16.

2

TESS-SPOC light curves accessed as a bulk download from https://archive.stsci.edu/hlsp/tess-spoc

3

See https://tess.mit.edu/tess-year-2-observations/

4

Docs can be found at https://docs.astropy.org/en/stable/timeseries/bls.html

5

https://heasarc.gsfc.nasa.gov/docs/tess/second-extended.html

APPENDIX A: STELLAR PARAMETERS FOR CANDIDATE SYSTEMS

Here, we provide details on the stellar parameters of the host stars for our candidate giant planet systems.

Table A1.

Stellar magnitudes for the host stars of our giant planet candidates. The values are from TICv8 (Stassun et al. 2019), GAIA DR2 (Gaia Collaboration 2018), and 2MASS (Skrutskie et al. 2006).

	TIC-95112238	TIC-311555090	TIC-46432937
TESS (mag)	12.2745 ± 0.0073	14.8993 ± 0.0077	12.3718 ± 0.0074
Gaia G (mag)	13.5155 ± 0.0006	16.2222 ± 0.0009	13.433 ± 0.0004
Gaia GBp (mag)	14.9506 ± 0.0022	17.8488 ± 0.0076	14.5196 ± 0.0025
Gaia GRp (mag)	12.3406 ± 0.0012	14.9837 ± 0.0027	12.3996 ± 0.0015
J (mag)	10.657 ± 0.021	13.236 ± 0.023	11.011 ± 0.022
H (mag)	10.063 ± 0.024	12.717 ± 0.025	10.427 ± 0.023
K (mag)	9.807 ± 0.02	12.468 ± 0.023	10.195 ± 0.02
	TIC-335590096	TIC-406717909	TIC-202468443
TESS (mag)	13.7726 ± 0.0079	13.5451 ± 0.0080	14.5982 ± 0.0074
Gaia G (mag)	15.1091 ± 0.0013	14.8626 ± 0.0009	15.9658 ± 0.0009
Gaia GBp (mag)	16.7462 ± 0.0062	16.454 ± 0.0063	17.7402 ± 0.0072
Gaia GRp (mag)	13.8593 ± 0.0035	13.6273 ± 0.0037	14.6966 ± 0.0014
J (mag)	12.056 ± 0.022	11.897 ± 0.021	12.818 ± 0.021
H (mag)	11.431 ± 0.026	11.3 ± 0.027	12.218 ± 0.023
K (mag)	11.175 ± 0.026	11.079 ± 0.023	11.965 ± 0.022
	TIC-165227846	TIC-60910638	TIC-73692250
TESS (mag)	13.4581 ± 0.0074	13.0687 ± 0.0077	13.9856 ± 0.0077
Gaia G (mag)	14.8075 ± 0.0005	14.3144 ± 0.0011	15.3646 ± 0.0009
Gaia GBp (mag)	16.5098 ± 0.0053	15.7465 ± 0.0077	17.1267 ± 0.0057
Gaia GRp (mag)	13.5507 ± 0.0015	13.1344 ± 0.0028	14.0843 ± 0.0028
J (mag)	11.743 ± 0.021	11.487 ± 0.02	12.239 ± 0.025
H (mag)	11.169 ± 0.022	10.894 ± 0.023	11.681 ± 0.025
K (mag)	10.885 ± 0.019	10.654 ± 0.014	11.409 ± 0.02
	TIC-243641947	TIC-77490011	TIC-67512645
TESS (mag)	13.2486 ± 0.0078	14.1398 ± 0.0074	14.9046 ± 0.0078
Gaia G (mag)	14.4858 ± 0.0007	15.4618 ± 0.0006	16.3092 ± 0.0012
Gaia GBp (mag)	15.9209 ± 0.0023	17.0973 ± 0.0036	18.1135 ± 0.0234
Gaia GRp (mag)	13.3138 ± 0.0032	14.2251 ± 0.0015	15.0078 ± 0.0026
J (mag)	11.706 ± 0.025	12.479 ± 0.021	13.169 ± 0.023
H (mag)	11.099 ± 0.024	11.930 ± 0.023	12.486 ± 0.022
K (mag)	10.819 ± 0.021	11.622 ± 0.019	12.207 ± 0.021
	TIC-178709444	TIC-429302040	TIC-254113311
TESS (mag)	13.6615 ± 0.0073	10.4180 ± 0.0061	10.1429 ± 0.0061
Gaia G (mag)	14.9535 ± 0.0005	11.1740 ± 0.0009	10.9028 ± 0.0008
Gaia GBp (mag)	16.5039 ± 0.0031	11.9254 ± 0.0012	11.6530 ± 0.0028
Gaia GRp (mag)	13.7388 ± 0.0014	10.3741 ± 0.0007	10.0917 ± 0.0014
J (mag)	11.999 ± 0.026	9.378 ± 0.021	9.055 ± 0.023
H (mag)	11.348 ± 0.028	8.777 ± 0.026	8.493 ± 0.059
K (mag)	11.084 ± 0.027	8.637 ± 0.023	8.351 ± 0.033

Table A2.

Stellar parameters for the host stars of our giant planet candidates. The parameters have been taken from the TICv8 (Stassun et al. 2019).

	TIC-95112238	TIC-311555090	TIC-46432937
RA (deg)	90.4837848985143	87.0953906689957	83.8690390629687
Dec. (deg)	−16.8166146303267	53.491840564709	−14.59734983288
Parallax (mas)	12.543 ± 0.045	11.224 ± 0.105	11.016 ± 0.025
Distance (pc)	79.72 ± 0.29	89.09 ± 0.84	90.77 ± 0.21
R* (R⊙)	0.557 ± 0.017	0.235 ± 0.007	0.539 ± 0.016
M* (M⊙)	0.552 ± 0.020	0.205 ± 0.020	0.535 ± 0.020
Teff (K)	3374 ± 157	3247 ± 157	3673 ± 157
log g	4.6878 ± 0.0098	5.0074 ± 0.0158	4.7037 ± 0.0093
	TIC-335590096	TIC-406717909	TIC-202468443
RA (deg)	183.5647447	309.47724827762	229.006270911694
Dec. (deg)	−13.17481731	19.5764691725682	62.7183462414147
Parallax (mas)	12.450 ± 0.099	10.122 ± 0.046	10.273 ± 0.063
Distance (pc)	80.32 ± 0.64	98.79 ± 0.45	97.35 ± 0.59
R* (R⊙)	0.353 ± 0.011	0.428 ± 0.013	0.307 ± 0.009
M* (M⊙)	0.338 ± 0.021	0.422 ± 0.020	0.285 ± 0.020
Teff (K)	3237 ± 157	3266 ± 157	3162 ± 157
log g	4.87111 ± 0.00054	4.8014 ± 0.0052	4.9177 ± 0.0045
	TIC-165227846	TIC-60910638	TIC-73692250
RA (deg)	178.819448968594	5.68952194396963	40.0813749354698
Dec. (deg)	−40.1491860945041	22.9455977613927	35.4105323079258
Parallax (mas)	16.038 ± 0.074	10.909 ± 0.073	11.749 ± 0.087
Distance (pc)	62.35 ± 0.29	91.67 ± 0.62	85.11 ± 0.63
R* (R⊙)	0.321 ± 0.010	0.467 ± 0.014	0.339 ± 0.010
M* (M⊙)	0.300 ± 0.020	0.465 ± 0.020	0.322 ± 0.020
Teff (K)	3202 ± 157	3373 ± 157	3163 ± 157
log g	4.9034 ± 0.0032	4.7660 ± 0.0068	4.8847 ± 0.0011
	TIC-243641947	TIC-77490011	TIC-67512645
RA (deg)	207.4749061	42.8241708040705	173.4697891
Dec. (deg)	−46.06623829	30.2817681183302	12.45099531
Parallax (mas)	13.707 ± 0.099	16.749 ± 0.068	13.52 ± 0.0968
Distance (pc)	72.96 ± 0.53	59.70 ± 0.24	73.96 ± 0.53
R* (R⊙)	0.372 ± 0.011	0.233 ± 0.007	0.2221 ± 0.0068
M* (M⊙)	0.360 ± 0.020	0.203 ± 0.020	0.191 ± 0.020
Teff (K)	3375 ± 157	3244 ± 157	3133 ± 157
log g	4.8525 ± 0.0018	5.0103 ± 0.0171	5.027 ± 0.019
	TIC-178709444	TIC-429302040	TIC-254113311
RA (deg)	166.0757629	188.386850399595	286.375931663504
Dec. (deg)	−47.82139085	−10.1461733018634	−41.4375228598884
Parallax (mas)	10.129 ± 0.044	15.418 ± 0.062	17.136 ± 0.050
Distance (pc)	98.72 ± 0.43	64.86 ± 0.26	58.36 ± 0.17
R* (R⊙)	0.427 ± 0.013	0.720 ± 0.065	0.742 ± 0.068
M* (M⊙)	0.421 ± 0.021	0.660 ± 0.080	0.660 ± 0.076
Teff (K)	3296 ± 157	4251 ± 121	4236 ± 122
log g	4.8023 ± 0.0053	4.5435 ± 0.1034	4.5169 ± 0.1023