Occurrence rate of hot Jupiters orbiting red giant stars

Hot Jupiters form an enigmatic class of object with yet unclear formation pathways. Determination of their occurrence rates as function orbit, planet and stellar mass, and system age, can be an important ingredient for understanding how they form. To date, various Hot Jupiters have been discovered orbiting red giant stars and deriving their incidence would be highly interesting. In this study we aim to determine the number of Hot Jupiters in a well-defined sample of red giants, estimate their occurrence rate and compare it with that for A, F and G-type stars. A sample of 14474 red giant stars, with estimated radii between 2 and 5 $R_\odot$, was selected using Gaia to coincide with observations by the NASA TESS mission. Subsequently, the TESS light curves were searched for transits from Hot Jupiters. The detection efficiency was determined using injected signals, and the results further corrected for the geometric transit probability to estimate the occurrence rate. Three previously confirmed Hot Jupiters were found in the TESS data, in addition to one other TESS Object of Interest, and two M-dwarf companions. This results in an occurrence rate of $0.37^{+0.29}_{-0.09}$%. Due to the yet large uncertainties, this cannot be distinguished from that of A-, F- and G-type stars. We argue that it is unlikely that planet engulfment in expanding red giants plays yet an important role in this sample.


Introduction
Hot Jupiters (HJs) are very interesting objects to study.The close proximity to their host stars makes them perfect candidates to observe with the transit method and causes their atmospheres to reach temperatures that are similar to some M-dwarf stars.However, it is still unknown how these objects form.To date, various hypotheses have been developed that could possibly explain their formation.The most promising are in situ formation, disk migration, and high-eccentricity tidal migration (Dawson & Johnson 2018).
One possible way to gain more information on the formation and the evolution of HJs is through the study of the occurrence rates as a function of planet and stellar masses, orbits, and system ages.Zhou et al. (2019) determined the HJ occurrence rate to be 0.43 ± 0.15 and 0.71 ± 0.31% for F-and G-type stars, respectively, while that around Kepler stars was estimated to be 0.43 ± 0.05 (Fressin et al. 2013), 0.43 +0.07  −0.06 (Masuda & Winn 2017), and 0.57 +0.14  −0.12 (Petigura et al. 2018).Zhou et al. (2019) found the occurrence rate around A stars to be 0.26 ± 0.11.Furthermore, Grunblatt et al. (2019) estimated the occurrence rate of giant planets around low-luminosity red giants to be 0.51 ± 0.29%.
It is particularly interesting to determine the occurrence rate of HJs around red giant stars.As evolved stars have gone through their initial hydrogen core burning phase, an estimate for the occurrence rate of HJs may tell us more about the evolution of these planetary systems and the survivability of close-in planets.To date, various HJs have been found to orbit red giants, Kepler-91b (Lillo-Box et al. 2014) being one of the most prominent.
In addition to HJs other kinds of exoplanetary populations have been found orbiting red giant stars as well.For example, both hot Neptune (e.g.K2-39b;Van Eylen et al. 2016) and warm eccentric Jupiter (e.g.Kepler-432b; Ciceri et al. 2015;Quinn et al. 2015) subpopulations have also been found.In addition, the detections of Kepler-56 b and c have shown that red giants can contain multiplanetary systems (Borucki et al. 2011;Huber et al. 2013).
In this paper we estimate the occurrence rate of HJs around red giants.In Sect. 2 we describe how we created a well-defined sample of red giants in the solar neighbourhood observed by the NASA TESS mission.In Sect. 3 we explain how we searched for HJs transiting these stars.Section 4 describes how we estimated the occurrence rate of HJs around red giants by combining the geometric probability of observing transits and the fraction of detected HJs.Our discussion and conclusions can be found in Sect. 5.

Gaia sample of red giant stars
All stars with available photometric measurements in the G, G BP , and G RP passbands, with estimated absolute magnitudes of M G ≤ 5.5 and with parallaxes of ϖ ≥ 2 milliarcseconds were selected using the Gaia DR2 catalogue (Gaia Collaboration 2016).To ensure that the sample contained only stars with sufficiently precise photometric and astrometric properties, the sample was subject to the following selection criteria, similar to those imposed in Gaia Collaboration ( 2018 First we remove all the stars with high astrometric excess noise by imposing the constraint where χ2 is the ASTROMETRIC_CHI2_AL parameter and ν ′ is the ASTROMETRIC_N_GOOD_OBS_AL parameter.The values for χ 2 and ν ′ were obtained from the Gaia DR2 catalogue.Next, to be able to correctly estimate the absolute magnitudes, we only retained the stars which have a precision of at leat 10% on their measured parallax, To only have stars in the sample with accurate photometric quantities, we imposed the following constrains on the signal-to-noise in the different Gaia passbands: As the Gaia G BP and G RP passbands were not corrected for blending of background sources.We remove the stars in our sample whose fluxes might suffer from contamination from background sources, all stars with PHOT_BP_RP_EXCESS_FACTOR falling outside the following two limits were removed from the collection: (4) Finally, we only retained the stars with well-measured astrometric quantities and remove those with a renormalised unit weight error of above 1.4.We reduced the sample further by ensuring that all stars had been observed during the first two years of TESS observations and for which it is expected that transits can be detected.First, we used the Python package TESS-POINT was used to check if the stars, based on their positional coordinates (right ascension and declination), are located in one of the first 26 TESS sectors.Next, we imposed an upper limit of R ≤ 10 R ⊙ on the stellar radii, as the transit depth scales inversely with the stellar radius squared.Furthermore, a lower limit of R ≥ 2 R ⊙ was imposed as we assumed this to be the typical stellar radius at the base of the red giant branch (RGB).Finally, we only retained the stars with apparent magnitudes of T ≤ 10, as the transit signals might be lost in the noise of the lightcurves of faint stars.The apparent magnitudes were estimated using Eq. 1 of Stassun et al. (2019): Finally, stars were selected based on their position in the colour-magnitude diagram (CMD) to remove stars from the sample that do not reside on the RGB.This criterion ensures that massive subgiant stars of similar radii do not contaminate the sample of red giant stars.The selection region has the following boundaries: The resulting CMD of the stars surviving the imposed selection criteria, together with the selection boundaries, is displayed in Fig. 1.The final sample consisted of 35 250 red giants.

TESS observations
If available, the two-minute cadence TESS Pre-search Data Condition SAP (PDCSAP; Jenkins et al. 2016) light curves were downloaded, otherwise the light curves were extracted from the 30-minute full frame images (FFIs) using the Python package TESSERACT 1 .About 83% of the light curves were obtained using TESSERACT and all these light curves were obtained using the default aperture option, which is based on the K2P2 method (see Lund et al. 2015).
The baselines of the light curves were detrended using an iterative spline fitting method (Vanderburg & Johnson 2014).The splines consisted of polynomials of degree 3, and in each iteration we excluded the 3σ outliers until no outliers remained.We further cleaned the light curves by removing all data points lying 3σ above the detrended baselines.Following Vanderburg et al. (2016), the two lowest data points were removed as well.

BLS search for hot Jupiters
To search the cleaned light curves for transits arising from HJ candidates on periods between 4 and 25 days, we used the BoxLeastSquares (BLS) implementation of the Python package ASTROPY (Astropy Collaboration 2013, 2018).The periodograms, yielded by the BLS method, were normalised following the steps given in Vanderburg et al. (2016).First, we subtracted the changing baseline level of the periodograms, to avoid spurious detections at longer periods.Next, we converted the BLS power measurements to signal-to-noise ratios by using Eq. 4 of Vanderburg et al. (2016), Here MAD stands for the median absolute deviation.
To validate any signals arising in the normalised periodograms, the periodograms were subject to a self-imposed algorithm, consisting of the following steps: first, we set a limit of S /N ≥ 9 on the highest peak found in the periodogram.If the highest peak did not meet this limit, it was assumed that the light curve did not contain any transits.Subsequently, the algorithm would consider the combination of transit parameters (orbital period, midpoint of first transit, transit depth, and transit duration) belonging to the signal yielding the highest peak in periodogram.
Secondly, if the corresponding transit duration was found to be longer than 25% of the orbital period, the found signal was rejected.If only a single data point was found in the fitted boxes of the BLS method, the signal was considered to be a false positive.In these cases, the algorithm would remove the data points yielding the signal from the light curves, generate a new periodogram, and return to step 1.
Finally, the radius of the candidate was estimated using where ∆F is the transit depth and R is the stellar radius.The candidate was assumed to be a potential The periodograms that passed all the steps were subsequently subjected to an individual visual inspection.The visual inspection allowed us to distinguish between real candidates and false positives.Among the false positives, we included binary stars, features in the light curves that clearly did not have a transit shape, and objects that were found to transit background stars.The objects orbiting background stars were investigated by thoroughly studying the TESS target pixel files.

Transit search results
The search for HJ candidates lead to 1392 potential detections, of which all except 6 were rejected upon visual inspection.The other 1386 detections were found to be false positives, binary stars, or objects orbiting stars in the background of the target star.The objects orbiting background stars were found by thoroughly studying the TESS target-pixel files and the light curves associated with each pixel.
As mentioned before, 83% of the light curves were obtained using TESSERACT.The light curves of TYC-5456-76-1, CD-25 2180, and both light curves of TYC-4004-1387-1 were obtained using TESSERACT.The light curves of KELT-11 and HD 221416 were both obtained from the Science Processing Operations Center (SPOC) pipeline, while one of the light curves of HD 1397 was obtained using the SPOC pipeline (sector 1) and the other using TESSERACT (sector 2).As the two pipelines resulted in similar numbers of observed candidates, both pipelines appear to reach similar sensitivities.Notes.Here, m dilution is a dilution factor and m flux is the mean out-oftransit flux; σ ω is a (unknown) jitter term; σ GP is the amplitude of the GP; and ρ denotes the Matern timescale.Notes.The subscript 'BLS' indicates that the value obtained from the BLS method was used.

Fitting the transits
To extract as much information about the found objects and their orbits, we fitted a transit model to the light curves using a Markov chain Monte Carlo (MCMC) algorithm, using the Python package EMCEE (Foreman-Mackey et al. 2013).The light curves were first re-detrended using a Gaussian processes (GP) method, which was implemented using the JULIET package (Espinoza et al. 2019).To ensure that the transit would not influence the GP trend, the out-of-transit data was removed by phase-folding the light curves, based on the period and midpoint of first transit found with our BLS search, and excluding all points that fall within 0.02 of the absolute phases.Subsequently, the GP trend was determined using the default Matern kernel; the hyperparameters used can be found in Table 1.The detrended light curves can be found in Appendix A.
The model light curves were created using the Python package PYTRANSIT (Parviainen 2015), assuming for simplicity linear limb-darkening with a value of 0.5 and circular orbits.The remaining parameters (orbital period, midpoint of first transit, ratio of object radius to stellar radius, semi-major axis in stellar radii, and inclination) were fitted with the MCMC algorithm.The a priori distributions for the various parameters are in Table 2.The best fitting transit models and the corresponding values for the transit parameters can be found in Appendix B.

Occurrence rate
To estimate the occurrence rate of HJs around red giants, we derived the geometric probability and the probability of actually detecting transits of the planets as functions of orbital period and planetary and stellar sizes.

Geometric probability
The geometric probability is the probability of observing a transit in a randomly oriented system.It depends mostly on the ratio of the semi-major axis of the orbit (a) to the stellar radius.The A26, page 3 of 12 A&A 670, A26 (2023) Table 3. Stellar masses, stellar radii, planetary radii, planetary masses, and semi-major axes used in the calculation for the geometric probability.
Star  geometric probability can be calculated using Eq. 2 in Kipping (2014), where after the right-hand equal sign, for simplicity the orbits are assumed to be circular by setting the eccentricity to e = 0 (see Table 3).
For our occurrence rate estimate, we used the average geometric probability of 0.196 +0.014 −0.057 .

Observational probability
To estimate the fraction of HJs detected, we performed artificial injection tests.The injected transits were based on the orbit of Kepler-91b (Lillo-Box et al. 2014).We carried out the injection tests for planetary radii of R p = 1.0 R Jup , 1.5 R Jup , and 2.0 R Jup , assuming an inclination of 90%.Our injected planets have periods between 4 and 15 days, randomly picked from a log-normal distribution.This period range covers the same range as for the detected planets in the geometric probability.
For each planetary radius we conducted five different tests, where in each test 150 stars were randomly chosen from our sample of red giants.Figure 2 shows the number of correctly detected and non-detected injected transits (top row), with the percentages of the correctly detected transits (bottom row).
As the confirmed exoplanets all orbit red giants with stellar radii of R < 5 R ⊙ , we only focus on these stellar radii in our calculation for the observational probability to ensure that both geometric and observational probabilities account for a similar population of red giant stars.This means that in the subsequent calculation of the total occurrence rate, we will only consider the 14 474 stars in our sample with radius R < 5 R ⊙ .To account for the fact that a greater number of smaller HJs with radii 1.0 R Jup ≤ R ≤ 1.5 R Jup have been observed than the larger HJs (R > 1.5 R Jup ), we have added weights based on the observed HJs with radii 1.0 R Jup ≤ R ≤ 2.0 R Jup and orbital periods P ≤ 15 days 2 .
We acquired a total number of 501 observed HJs, of which 323 have radii between 1.0 R Jup ≤ R ≤ 1.33 R Jup , 136 have radii between 1.33 R Jup ≤ R ≤ 1.66 R Jup and 42 have radii between 1.66 R Jup ≤ R ≤ 2.0 R Jup .Based on these values, we establish that about two times more HJs have radii close to 1.0 R Jup compared to 1.5 R Jup and three times more HJs have radii close to 1.5 R Jup compared to ∼2.0 R Jup .Subsequently, we determined our fraction of observed HJs using a weighted average, where we used weights of 2, 1, and 1 3 , respectively.Our weighted averaged yields a percentage of 0.38 ± 0.13 detected HJs. 2 The observed hot Jupiters were acquired from Exoplanet.eu(http: //exoplanet.eu/,on July 18, 2022) using the query: RADIUS:RJUP >=1 AND RADIUS:RJUP <= 2 AND PERIOD:DAY <= 15.

Total probability and occurrence rate
Combining the geometric and detection probability, we find the combined probability of seeing a HJ with our method to be 0.07 ± 0.03.In our search for HJs orbiting red giants, we find four systems.Correcting for the geometric probability and detection fraction, we estimate that there are 40-96 planets orbiting our sample stars.For our sample of 14 474 red giants with R ≤ 5 R ⊙ , this yields an occurrence rate of 0.37 +0.29  −0.09 .

Occurrence rate
The occurrence rate derived here of 0.37 +0.29 −0.09 % agrees, within the uncertainties, with those for the A-, F-, and G-type stars of, respectively, 0.26 ± 0.11, 0.43 ± 0.15, and 0.71 ± 0.31% (Zhou et al. 2019).These occurrence rates were derived for HJs with orbital periods between 0.9 and 10 days.No clear distinction can yet be made between the occurrence rates of these different stellar types.The occurrence rate also agrees with that for HJs orbiting the Kepler stars of 0.43 ± 0.05% (0.8 ≤ P ≤ 10 days; Fressin et al. 2013), 0.43 +0.07 −0.06 (P ≤ 10 days Masuda & Winn 2017), and 0.57 +0.14 −0.12 (1 ≤ P ≤ 10 days Petigura et al. 2018).Finally, our found value is similar to that reported by Grunblatt et al. (2019) for giant planets orbiting low-luminosity red giants (P ≤ 10 days) in the K2 observations, 0.51 ± 0.29%.All in all, the found occurrence rate for HJs orbiting red giant stars cannot (yet) be distinguished from the occurrence rates reported for the above-mentioned stellar types.
The masses of the three host stars already discussed in the literature are estimated to be in the range 1.2−1.5 M ⊙ , meaning that they used to be F-type main-sequence stars.As mentioned above, Zhou et al. (2019) found an occurrence rate of 0.43 ± 0.15 for these types of stars.A larger sample and/or more TESS data are needed to investigate the potential differences between these populations.
In addition, when more HJs have been found orbiting red giants, the occurrence rate can also be improved by including the non-zero eccentricities of the detected planets in the calculation of the geometric probability and in the injection-retrieval tests.

Possible planet engulfment
A reason that might explain a lower occurrence rate of HJs orbiting red giants and the low number of planets that have been detected is planet engulfment (Lillo-Box et al. 2016).Planet engulfment follows from the tidal dissipation of the orbit of the planet, due to interactions between the planet and the convective envelope of the red giant.
A26, page 4 of 12 We argue, however, that planet engulfment is not yet of importance during these early stages of the red giant phase.Figure 3 shows the orbital evolution of Jupiter-mass planets, on various initial semi-major axes, during the early red giant phase of a star of mass 1.5 M ⊙ .
The evolution was modelled using the code Modules for Experiments in Stellar Astrophysics (MESA v15140;Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015Paxton et al. , 2018Paxton et al. , 2019)).The star was evolved from the pre-main-sequence phase up to the core helium burning phase, the end of the RGB phase.During the evolution, the mass loss was estimated using Reimers' empirical formula (Reimers 1975): Here, L, M, and R are the stellar luminosity, mass and radius, respectively, all given in solar units.η is a parameter of order unity, whose value usually varies between 0.2 and 1.0.A lower value of η is often used for metal-poor stars as the coupling of photons to the stellar gas is lower when there are fewer heavy elements present (Kippenhahn et al. 2012).In the simulations we used a fixed value of η = 0.5.The orbital decay of the planet was modelled using the rate of change of the semi-major axis (Eq.(1) in Kunitomo et al. 2011): Here M p denotes the mass of the orbiting planet, a is the semimajor axis of the orbit, and k is a dimensionless factor that can be calculated using Eq. ( 2) of Kunitomo et al. (2011): Here M env is the mass that is contained within the convective envelope of the RGB star, while P denotes the orbital period.The orbital period was calculated using Kepler's third law: Furthermore, τ d is the eddy turnover timescale, which was calculated using Eq. 4 of Kunitomo et al. (2011), where R env is the radius at the base of the convective envelope.
The simulation of the orbital decay was halted as soon as the semi-major axis became smaller than the stellar radius.
As can be seen, the orbital decay for different initial semimajor axes becomes important once the star starts filling up a large portion of the planetary orbit.None of the confirmed exoplanets in this sample are experiencing significant orbital decay and they are, currently, not in danger of getting engulfed.For Kepler-91b, on the other hand, orbital decay might become important within tens of Myr (Lillo-Box et al. 2014), which is also visible in Fig. 3.In addition, ultra-HJs (periods <2 day), which account for a small subsample of the HJ population, will likely not have survived this early red-giant phase.It should be noted that orbital decay becomes more important for heavier planets around less massive stars, as can be inferred from Eq. ( 11).In addition, the confirmed exoplanets detected in this work shown in Fig. 3 are all less massive than 1 M Jup (see Table 3).Subsequently, their orbital decay will be even less severe than what is suggested in Fig. 3.

Future searches
The occurrence rate estimate presented in this work is based on only a few objects and may be prone to biases due to the detection probability being a strong function of stellar apparent magnitude and stellar radius.Since TESS continues to observe, the available data to search for HJs around red giant stars will only increase and mitigate this issue.It may also help to conduct a more sensitive study using difference imaging analysis (Oelkers & Stassun 2018;Montalto et al. 2020) or a principle component analysis (see, e.g., the eleanor Feinstein et al. 2019;or the Giants Saunders et al. 2022 pipelines) to the light curves obtained from the FFIs.These pipelines could yield cleaner light curves, an increase in detection probability, and a more precise estimate for the occurrence rate.For example, our method has missed one possible candidate, TOI-4551.013, due to transit-like signals in the light curve, which are likely related to the momentum dumps occurring during the observations.Using one of these two pipelines could lead to the inclusion of such candidates.
In the future, the ESA PLATO mission (Rauer 2017) may result in high-precision light curves for a large sample of red giant stars from which HJs can be discovered and their occurrence rate further constrained.
): A26, page 1 of 12 Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0),which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.This article is published in open access under the Subscribe-to-Open model.Subscribe to A&A to support open access publication.

Fig. 1 .
Fig. 1.Colour-magnitude diagram of stars that have survived the imposed selection criteria.The black box represents the boundaries of the region enclosing the RGB.The colour (G BP − G RP ) of the stars can be found on the x-axis, while the y-axis denotes the absolute magnitudes.The logarithmic colour map indicates the number of stars in each two-dimensional bin.Stars removed by the selection criteria are shown in grey.

Fig. 2 .
Fig. 2. Results from the injecting artificial transit test.Top row: result for injected planets with radius R p = 1.0 R Jup (left), for planets with radius R p = 1.5 R Jup (middle), and for planets with radius R p = 2.0 R Jup (right).The green shaded areas shows the correctly detected artificial transits for different stellar radii, while the red shaded areas show the non-detections.The black line shows the total number of injected planets for the different stellar radii.Bottom row: percentages of correctly detected artificial transits for stellar radii in the range 2. R ⊙ ≤ R ≤ 5 R ⊙ .The red line indicates the average percentage and the red shaded area indicates the 1σ errors.

Fig. 3 .
Fig.3.Expected orbital evolution of Jupiter-mass planets orbiting a red giant of mass 1.5 M ⊙ , on orbits with initial semi-major axis within the range 0.02-0.10AU, followingKunitomo et al. (2011).The HJs in our sample are indicated; Kepler-91 is included for comparison.

Table 1 .
Distributions and hyperparameters used during the GP method.

Table 2 .
A priori parameter limits used in the MCMC algorithm.

Table B .
3. Same as TableB1, but for the star HD 221416.

Table B .
5. Same as TableB1, but for the star CD-25 2180.