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

We present a systematic search for transiting giant planets ($0.6 R_{\rm J} \leq R_{\rm P} \leq 2.0 R_{\rm J}$) orbiting nearby low-mass stars ($M_{*} \leq 0.71 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 minute 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 \pm 0.097$% (0.088 - 0.26 $M_{\odot}$), $0.108 \pm 0.083$% (0.26 - 0.42 $M_{\odot}$), and $0.29 \pm 0.15$% (0.42 - 0.71 $M_{\odot}$). For our full sample (0.088 - 0.71 $M_{\odot}$) we find a giant planet occurrence rate of $0.194 \pm 0.072$%. We have measured for the first time the occurrence rate for giant planets orbiting stars with $M_{*} \leq 0.4 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.


INTRODUCTION
Since the detection of 51-Pegasi b (Mayor & Queloz 1995) revealed the existence of giant planets on short orbits ( 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 (eg.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 2005a;Osborn & Bayliss 2020).It is believed that the metallicity of a star is inherited from its primordial cloud and thus shared by its protoplanetary disk (Fischer & Valenti 2005b), and it has been shown that gas giant planets from through core-accretion more efficiently in these metal-rich protoplanetary disks (e.g.Ida & Lin 2004;Thommes et al. 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 ( eff = 4100 − 6100 K and log  = 4.0 − 4.9), Howard ★ E-mail:edward.bryant@ucl.ac.uk 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 Stellar Mass [M ] ).The sample shown is made up of planets with masses between 0.1 and 13  J and radii larger than 0.6  J .
et al. ( 2012) derived a frequency of 0.4 ± 0.1% for planets with  P ≥ 8  ⊕ and  < 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% for planets with 6  ⊕ ≤  P ≤ 22  ⊕ and 0.8 ≤  ≤ 10 d. currence rates around stars with masses  * ≤ 0.7  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  ≤ 10 d.Therefore, they could only place an upper limit of 3% 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 ⊕ ≤  P ≤ 2 J ; 0.8 ≤  ≤ 10 d) around early M-dwarfs (0.45 − 0.65 ) to 0.27 ± 0.09%.
In this paper we present a systematic search for transiting gas giant planets (0.6 J ≤  P ≤ 2.0 J ) orbiting nearby ( ≤ 100 pc), lowmass (0.088  ≤  * ≤ 0.71  ) stars in the TESS 30 minute 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 d ≤  ≤ 10 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.

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 We also show the number of TESS sectors for which as TESS-SPOC 30 minute cadence light curve exists for each object, plotting this histogram on a log scale.observed for 27.4 days.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 - P ≤ 4  ⊕ -amenable to further characterisation including with JWST (Gardner et al. 2006).The discovery of these small exoplanets was primarily achieved through 2 minute 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 minutes.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 minute 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 Product2 .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: (i)  eff ≤ 4500 K, (ii)  * ≤ 0.75  , (iii)  ≤ 100 pc, and (iv)  ≤ 16 mag, where  eff is the effective temperature of the star,  * is the stellar radius,  is the distance of the star from the Sun, and  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 - (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 minute cadence FFI light curves yields 91,306 stars with a light curve.We present our low-mass star sample in Figure 2. From this plot, as well as the histogram in Figure 3, we see that our lowmass star sample spans a wide range of spectral types and stellar masses (≈ 0.1 − 0.7  ), with the majority of a sample being lowermass than 0.4  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 ( 14 mag) and small ( * 0.25  ) 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.
From the histogram in Figure 2 we see that the majority of the stars in our sample have the spectral type M2 to 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 ( ≤ 16 mag) and the very low luminosity of M5 and later spectral types stars.
The apparent magnitude distribution of our sample peaks around  = 13 mag (see top-left histogram in Figure 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  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 Figure 3, arises from the TESS observing strategy.As can be seen from the TESS observations footprint 3 , 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.

TRANSIT SEARCH
Using the TESS-SPOC 30 minute cadence light curves, we firstly exclude any data points with a quality flag  > 0. These flags highlight any observations which are adversely affected by cosmic rays and scattered light, among others.For this anal ysis 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 (Stumpe et al. 2012(Stumpe et al. , 2014;;Smith et al. 2012).To remove stellar variability prior to the planet search we smooth the light curves using a Savitzky-Golay filter, using the implementation from the package (Lightkurve Collaboration et al.

2018
).The filter uses a window length of 45 data points, approximately 1 d for 30 minute 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 et al. 2002).We use the implementation of the algorithm4 .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 (SNR) of greater than 8 as good transiting planet candidates.In total our search yields the detection of 3930 periodic transit-like signals.

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 (eg.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 % 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 summarise 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.

Secondary Eclipse Events
One of the main astrophysical false positive cases for transiting planet searches are eclipsing binary stars (eg.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 % 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  J planet in a 1.5 d orbit around a 0.4  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.

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% 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 misidentified 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.

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 Figure 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 Figure 4.
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 Figure 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) and fit for the coefficients A, B, C for each light curve.We then compare the quality of this harmonic fit to a simple flat line fit by computing the value 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  2 Harmonic > 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 Figure 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 et al. 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 (eg.Rodríguez & Breger 2001) and RR-Lyrae (eg.Simon & Teays 1982), and can exhibit variability with an amplitude on 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 onto the target light curve.

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 to 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 Figure 5.A total of 8 objects are identified as variable stars using this analysis.

Lomb-Scargle Power
We also check for continuous variability by performing a Generalised Lomb-Scargle (Lomb 1976;Scargle 1982) analysis on the light curve, after masking out the flux data points which are detected as "intransit" by BLS.Any star with a normalised 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 Figure 5.A total of 1039 objects are identified as variable stars using this test.This check alone removes 26.4% of the transiting candidates detected by the BLS search.However when applied to the simulated planet light curves just 0.22% 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.

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 root-mean-square of the point-to-point scatter.For a transiting planet light curve, the out-oftransit 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 √ 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 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 Figure 5.A total of 12 objects are identified as variable stars through this analysis.

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.

TRANSIT FITTING ANALYSIS
Having identified and rejected the false positives, we then fit the transit events for the remaining candidates.We use the 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 midtransit time,  C , the planetary orbital period, , the planet-to-star radius ratio,  P / * , the scaled semi-major axis, a/ * , and the orbital inclination, .We also fit for a free flux baseline offset,  0 , which is defined such that the out-of-transit flux is equal to 1 +  0 .We use a quadratic limb-darkening law and for the limb-darkening coefficients we fit for the  1 ,  2 parameters from the parameterisation of Kipping (2013).This parameterisation ensures a physically realistic limbdarkening 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.
For each object, 24 independent chains were each sampled for 7,500 steps, following a burn-in phase of 2,500 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-fit" parameters.We include only objects with a best-fit planetary radius in the range 0.6  J ≤  P ≤ 2.0  J in our candidate list.Using the transit models computed using the best-fit parameters we also calculate the ratio  trans =   /, where   is the transit duration.We accept only objects with 0.0 <  trans ≤ 0.1 as likely planetary candidates.We test the threshold chosen for this criterion by calculating the  trans value for all the planets we simulate for the Injection-Recovery tests (see Section 8).The highest  trans value is 0.079, with 90% of the simulated planets having  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  P using the  P / * posterior and the host star radius from the TIC.Similarly, we use the posteriors for  and a/ * along with Kepler's 3 rd 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  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  J ≤  P ≤ 2.0  J and 1.5 g cm −3 ≤  * ≤ 200 g cm −3 .We then exclude from our candidate list any object with this fraction  GP < 15%.From these automated steps outlined above, we identified a list of 44 candidate giant exoplanets.

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.

Visual Inspection
For the first manual check we visually inspected the TESS-SPOC 30 minute 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.

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  J ; Huang et al. 2020b) and WASP-107 b (TIC-429302040; 0.12  J ; Anderson et al. 2017).The third object is the known brown dwarf LP 261-75 b (TIC-67646988; 68  J ; Irwin et al. 2018).Therefore we exclude TIC-67646988 from our candidate list.

Blend Scenario Checks
Due to the large pixel scale of TESS (21 /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 Python software (Lightkurve Collaboration et al. 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.

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 minute cadence SPOC light curve, from which we determined the following about these candidates.The 2 minute 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 minute cadence photometry reveals more clearly the shape of the transit events.By refitting this photometry we find a best fit radius of 2.59  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 minute cadence alone they could have been giant planets on grazing orbits.The 2 minute cadence photometry shows that they are not grazing, and so are most likely real planets but with  P < 0.6  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( P = 2.58 ± 0.24  ⊕ ) and TOI-805.01 ( P = 1.28 ± 0.62  ⊕ ).Based on this we exclude them both from our giant planet candidate list.For TIC-178709444 the 2 minute cadence light curve contains no evidence that the candidate is a false positive.In fact the higher cadence confirms the transit as non-grazing ( < 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 minute cadence full-frame-image 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 minute 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 minute cadence photometry revealed the transit as being flat-bottomed, whereas in the 30 minute 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 minute cadence photometry confirmed the transits as flat-bottomed, further supporting the nature of these two objects as planetary candidates.The 10 minute 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.

GIANT PLANET CANDIDATES
From our full planet search we have identified a final selection of fifteen giant planet candidates.The numbers of objects from our sample at each step of the planet search pipeline are summarised in Figure 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 Figure 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.
We also compare our sample of candidates to the full population of known transiting giant planets in Figures 8 and to the subset of known transiting giant planets with low-mass host stars in 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 Figure 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  * < 0.4  , with some as low-mass as  * ≈ 0.2  .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).
From the colour-magnitude diagram in Figure 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-normalised 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;  P sin  = 0.36  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.

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 summarised 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.parameters to generate the light curve model.The stellar radius,  * , and  P were used to derive the radius ratio,  P /  * .The stellar mass,  * , along with  * and  were used to derive the scaled semi-major axis of the orbit,  /  * using Kepler's Third Law, assuming circular orbits.Finally,  /  * and  were used to calculate the orbital inclination, .Note that we assume circular orbits for all simulated planets.
Using this method, we simulated 873760 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% of the injected period were then passed to the vetting steps we outline in Section 4. Our pipeline yielded 780379 significant BLS detections from our simulated planet population, of which a total of 759632 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 759632 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 Figure 10.We can see that across the majority of the {;  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  to the Gaia colour-magnitude distribution of the 91,306 low-mass stars in our sample (density heat map).We also plot the twelve known transiting gas giant planets with hosts stars with masses ≤ 0.65  as the blue squares for comparison to our candidates.
for these planets due to the smaller number of total transit events in the light curve.
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, , for these planets, leading to non-grazing transits being fit with grazing transit models and vice-versa.The 30 minute cadence of the TESS FFI observations reduces our knowledge of the true transit shape causing these inaccurate fits, resulting in either an over-or under-estimation 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 minute 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  J and 2.0  J are 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 FFPs for the currently un-dispositioned candidates.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 fit parameters for these simulated planets outside the giant planet regime.
From the overall detection efficiency of our pipeline, which is plotted in Figure 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  ≤ 3.5 d.

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,  occ , using t he following equation where  pl is the number of planets detected and  pr is the number of stars amenable to the discovery of an exoplanet (see e.g.   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. Sackett 2011; Gan et al. 2022).This factor is calculated as follows where  * is the number of stars in the stellar sample,  sim is the total number of planets simulated,  fit is the number of simulated planets which were passed to the transit fitting,  vet, 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,  is a similar delta function dependent on whether a simulated planet was identified as a giant planet candidate by the transit fitting, and P tr, is the geometric transit probability of the simulated planet.This transit probability is defined as It is clear from Figure 11 that both the detection efficiency and planet occurrence rates will be functions of the orbital period, , and planet radius  P .Therefore, to calculate the overall giant planet occurrence rate we first calculate individual  occ values for each of the cells in {;  P } parameter space in Figure 11, and then calculate the sum of these individual values.Similarly to estimate the uncertainty on the occurrence rate,   occ , for each cell in {;  P } parameter space in which at least one planet candidate is detected we calculate the uncertainties on  pl and  pr as √  pl and √︁  pr respectively, as estimated from Poisson counting statistics.The individual cell occurrence rate uncertainty is then computed for each cell as 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 } bin, the value used for  pl is calculated as where  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   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 ≤  * ≤ 0.71 ) we calculate an occurrence rate of close-in giant planets of 0.194 ± 0.072 %.
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.also includes a prior on planet radius which penalises 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 %.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  .We calculate occurrence rates for each sub-sample finding values of 0.137 ± 0.097 % (0.088 -0.26  ), 0.108 ± 0.083 % (0.26 -0.42  ), and 0.29 ± 0.15 % (0.42 -0.71  ).We plot these occurrence rates, along with results from previous studies, in Figure 12, and occurrence rates we calculate in this work are summarised 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 % lower limits for the occurrence rates in each mass range.Doing so we find lower limits of 0.048 % (0.088 -0.26  ), 0.042 % (0.26 -0.42  ), and 0.10 % (0.42 -0.71  ).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 %, with 95 % lower limits of 0.029 %, 0.021 %, and 0.085 %.Therefore, our results show that a population of close-in giant planets exists even for the lowest range of stellar masses we study.
The distribution of our candidates in Figure 11 suggests the possibility of a dependence of the occurrence rate on the planet radius,  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  P .These results are plotted in Figure 13 and we see that the occurrence rates show clear variation with  P .We see a peak in occurrence rates for  P ∼ 1 J and occurrence rates consistent with zero for  P ≥ 1.4 J .Whereas for more massive host stars - * ≥ 1.2  -there exist transiting planets with radii  P > 1.5  J (see Figure 8).It has been previously suggested that extreme levels of irradiation of the planet by the star, Figure 12.Dependence of the occurrence rate of giant planets with short orbital periods ( ≤ 10 ) 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 % 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.on the order of 10 5 −10 6 W m −2 , are required to cause the inflation of the planetary atmosphere necessary to reach planetary radii 1.1 J (Sestovic et al. 2018).A planet in a 3 d orbit around a 0.45 star will receive an irradiation of around 5.7 × 10 4 W m −2 .Therefore, it is possible that the derived lack of high radius giant planets for lowmass 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 Figure 11 shows that the majority of our candidates have short periods < 4 .However, the detection efficiencies are also much higher for short orbital periods.Plotting the occurrence rate as a function of orbital period in Figure 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.
We compare our results to some of these previous studies in Figure 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  ≤ 1  .Considering our result from our full sample -0.194 ± 0.072 % -we find it is less than the Petigura et al. (2018) result at a level of 2.69  and less than the Beleznay & Kunimoto (2022) 0.8 − 1.05  result at a level of 2.27 .Considering the case in which all nine of the remaining candidates are real planets, our upper limit occurrence rate of 0.267 ± 0.079 % is still less than both of these results, at a level of 2.10  for Petigura et al. (2018) and a level of 1.76  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  selection of our stellar sample -0.29 ± 0.15 % -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  * ≤ 0.4 we find an occurrence rate of 0.134 ± 0.069 %, which is less than the Gan et al. (2022) result at a level of 1.15 .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.

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  * < 1  and comparing the results for our full sample to those for solar-like stars from Kepler (Petigura et al. 2018) and TESS (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 coreaccretion formation theory, also predicts that the formation of giant planets becomes impossible for stars with  * ≤ 0.4  .Studies of protoplanetary disks 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 disk 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  * ≤ 0.4  of 0.134 ± 0.069 %, with a 95 % lower limit of 0.05 %.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  * < 0.4 .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 % with a 95 % lower limit of 0.026 %.
There must be some pathway through which these stars can form giant planets.Firstly, we consider that the current main inhibitor for the predicted formation of these planets is the limits on the mass of the protoplanetary disk.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 disk-hosting low-mass stars could allow the occurrence of such massive disks to be constrained.It is important here to note that the dust masses for protoplanetary disks 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 disks.We also note that to date protoplanetary disks 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 disks 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  (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 2  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.

CONCLUSIONS AND FUTURE OUTLOOK
We have presented a systematic search through TESS 30 minute 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 lowmass 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 Figure 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 % 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 % (0.088 -0.26  ), 0.108 ± 0.083 % (0.26 -0.42  ), and 0.29 ± 0.15 % (0.46 -0.71  ), 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  (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 coreaccretion 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 .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 - * ≤ 0.4  -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 disks 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  * < 0.4 .Therefore the new generation of stabilised 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 realised the limitations of the long cadence (30 minute) observations in characterising 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 minutes in the first extended mission and then 200 seconds in the second extended mission5 -will greatly improve our ability to classify giant planets transiting low-mass stars.

Figure 1 .
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  J and radii larger than 0.6  J .

Figure 2 .
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 .
Figure3.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 minute cadence light curve exists for each object, plotting this histogram on a log scale.

Figure 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  = 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

Figure 5 .
Figure5.Examples of false positives identified by the variability analysis checks detailed in Sections 4.4 to 4.5.3.The titles of each panel gives 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.

Figure 7 .
Figure 7.Light curves for the fifteen 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 fit models.The TIC IDs and orbital periods are given for each system.

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

Figure 9 .
Figure9.Comparison of our fifteen 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 twelve known transiting gas giant planets with hosts stars with masses ≤ 0.65  as the blue squares for comparison to our candidates.

Figure 10 .
Figure10.Two-dimensional 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 percent.The values here highlight the sensitivity of our pipeline to transiting giant planets.

Figure 11 .
Figure11.Two-dimensional 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 percent.The stars highlight the locations in {;  P } parameter space of our fifteen giant planet candidates.

Figure 13 .
Figure13.Dependence of our measured occurrence rate of giant planets with short orbital periods ( ≤ 10 ) and low-mass host stars on planet radius.Dashed horizontal line represents the 90% confidence level upper limit for the bin in which there are no planet candidates.

Figure 14 .
Figure 14.Dependence of our measured occurrence rate of giant planets with short orbital periods ( ≤ 10 ) and low-mass host stars on the orbital period of the planet.

Table 1 .
C U ( C, BLS − 0.1  BLS ;  C, BLS + 0.1  BLS ) Priors used for each free parameter in the transit fitting detailed in Section 5. Note U ( ; ) denotes a uniform prior with lower bound A and upper bound B.

GP ≥ 0.15 44 Giant Planet Candidates Manual Candidate Vetting and Additional Obs Sections 6 & 7 5 Confirmed Planets 1 Eclipsing Binary identified with spectroscopy 9 Giant Planet Candidates remaining 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.

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

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