A sub-Neptune transiting the young field star HD 18599 at 40 pc

HD 18599 was observed by TESS in sectors 2 (2018 Aug 22-Sep 2 UT; see Figure 1) and 3 (2018 Sep 20-Oct 18 UT) during its first year with 2-minute cadence, and sectors 29 (2020 Aug 26-Sep 22 UT) and 30 (2020 Sep 22-Oct 21) during its third year with 20-sec cadence. The photometric data were processed by the Science Processing Operations Center (SPOC; Jenkins et al., 2016) data reduction pipeline which produced time series light curves for each campaign using Simple Aperture Photometry (SAP), and the Pre-search Data Conditioning (PDCSAP) algorithm (Stumpe et al., 2012; Smith et al., 2012). For our analysis, we used the PDCSAP light curves which have been corrected for the instrumental and systematic errors as well as light dilution effects from nearby stars. The time series in each sector has a $\sim$4-day gap because of the data downlink and telescope re-pointing.

The TESS object of interest (TOI) releases portal²²2https://tess.mit.edu/toi-releases/ announced the transiting planet candidate around HD 18599 (TIC 207141131) as TOI 179.01 (Guerrero et al., 2021). The candidate passed all tests from the Alerts Data Validation Report (Twicken et al., 2018) and is listed on the Exoplanet Follow-up Observing Program (ExoFOP)³³3https://exofop.ipac.caltech.edu/tess/ webpage as having a period of 4.1374 d and a transit depth of about 1097 ppm.

We performed independent transit search in the PDCSAP using transit-least-squares ⁴⁴4https://github.com/hippke/tls (Hippke & Heller, 2019) after applying a biweight filter using wotan ⁵⁵5https://github.com/hippke/wotan with a window length of 0.5 days to detrend the stellar variability in the concatenated light curves in sectors 2 & 3, separately from sectors 29 & 30 due to the long gap. A transit signal with 4.13 d with a signal detection efficiency (SDE; Pope et al., 2016) of 17 and a 1.75 hr duration was detected. No other significant periodic transit signal was detected in further iterations.

We obtained observations of HD 18599 with the intention of following up the candidate planet using the Infrared Array Camera (IRAC, Fazio et al., 1998) 4.5$\mu m$ channel, as part of Spitzer TESSTOO program 14084 (P.I. Crossfield). Observations with Spitzer have a number of advantages over those taken by TESS. First, Spitzer has a smaller pixel scale (1.2 arcsec/pix) than that of TESS’s (21 arcsec/pix). This allows Spitzer observation to localize the signal by excluding the nearby stars that are blended with the target star in the TESS photometric aperture (see Figure 1). Second, the effects of limb darkening in Spitzer is reduced because it operates in the near-infrared, relative to TESS which operates in the optical. Third, Spitzer has better sampling of the transit due to the shorter (2-second) cadence of our Spitzer data compared to the shortest available TESS (20-second) cadence. Hence, more accurate planet parameter estimates are obtained when modeling the Spitzer transit light curve jointly with TESS. In conjunction with the TESS bandpass, the 4.5$\mu m$ IRAC bandpass also provides a relatively broad wavelength baseline which is very useful for planet validation (See Section 4.3).

We used integration time of 2 seconds to keep the detector from saturating and minimize data downlink bandwidth. Following Ingalls et al. (2016), the target was placed on the "sweet spot" of the detector ideal for precise time-series photometry of bright stars like HD 18599. We then extracted the Spitzer light curves following Livingston et al. (2019) which was based on the approach taken by Knutson et al. (2012) and Beichman et al. (2016). In brief, we compute aperture photometry using circular apertures centered on HD 18599, for a range of radii between 2.0 and 5.0 pixels, corresponding to 2.4"–6.0". We used a step size of 0.1 pixel from 2.0 to 3.0, and a step size of 0.5 from 3.0 and 5.0. The optimal aperture was selected by minimizing the photon noise due to sky background and correlated noise due to inter- and intra-pixel gain variations. The resulting light curve is shown in the last panel of Figure 2.

We also conducted ground-based follow-up transit observation of HD 18599 on 2018 Oct. 16 using the SIRIUS camera (Nagayama et al., 2003) on-board the 1.4-m Infrared Survey Facility (IRSF) telescope located in Sutherland, South Africa. The instrument is capable of simultaneous imaging in $J,H,K_{s}$ bands which is ideal for planet validation. We created light curves in three bands using the standard reduction method and aperture photometry following Narita et al. (2013). Although, we were not able to detect the shallow 1.1 ppt (parts per thousand) event on target in all bands, we were able to rule out the deep eclipses from nearby faint stars that could reproduce the TESS detection. This adds further evidence that the signal indeed originates from HD 18599 as shown in Section 2.1.2.

We observed HD 18599 on 2021 Sep 18 with 1,0m LCO-CTIO in the B and $z_{s}$ bands, on 2018 Nov 18 with 1.0m LCO-CPT in the $g^{\prime}$ band, on 2018 Nov 23 with the 1.0m LCO-CTIO in the $u^{\prime}$ band, and on 2018 Dec 22 with the 0.4m LCO-CTIO in the $i^{\prime}$ band. All observations were full transit except that one taken in $i^{\prime}$ band. We scheduled our transit observations using the TESS Transit Finder, which is a customized version of the Tapir software package (Jensen, 2013). The photometric data were calibrated and extracted using AstroImageJ (Collins et al., 2017). Comparison stars of similar brightness were used to produce the final light curves, each of which showed a roughly 2-ppt dip near the expected transit time. The observations are summarized in Table 1 and plotted in Figure 4.

We observed HD 18599 on 2019 Nov 21 using NGTS (Next Generation Transit Survey) based at ESO’s Paranal Observatory in Chile. This array of twelve 0.2m telescopes is equipped with 2K$\times$2K e2V deep-depleted Andor Ikon-L CCD cameras with 13.5$\mu$m pixels, corresponding to an on-sky size of 4.97,̈ each with a custom NGTS (550-927 nm) filter. These observations were performed in the ’multi telescope’ observing mode (Smith et al., 2020), using three of the NGTS telescopes simultaneously. The data were reduced using a custom aperture photometry pipeline (Bryant et al., 2020), which uses the SEP library for both source extraction and photometry (Bertin, E. & Arnouts, S., 1996; Barbary, 2016).

We observed HD 18599 in the $I_{C}$ band from the Perth Exoplanet Survey Telescope (PEST) near Perth, Australia. The 0.3 m telescope is equipped with a 1530$\times$1020 SBIG ST-8XME camera with an image scale of 1.2" pixel^-1 resulting in a 31’$\times$21’ field of view. A custom pipeline based on C-Munipack⁶⁶6https://c-munipack.sourceforge.net was used to calibrate the images and extract the differential photometry.

HD 18599 was observed as part of the Kilodegree Extremely Little Telescope (KELT) survey using a 42mm-aperture telescope located in Sutherland, South Africa (Pepper et al., 2012). The telescope is equipped with Mamiya 645 80mm f/1.9 42mm lens with a 4k$\times$4k Apogee CCD, a pixel scale of $23\arcsec$ and a field of view of $26^{\circ}\times 26^{\circ}$. The target was observed with with a 20-30 minute cadence. The data were reduced in a standard manner using the pipeline described in detail in Siverd et al. (2012) and Kuhn et al. (2016).

We obtained 9 spectra with typical exposure time of 600 s conducted between 2019 Sep 10 and 19 UT using FEROS echelle spectrograph with a resolution R=48000, wavelength coverage between 350 and 920nm, and mounted on the MPG/ESO-2.2m telescope in La Silla, Chile (Kaufer et al., 1999). The spectra have typical SNR of 80. The observations were performed in the Object-Calibration mode to allow precise RV observations. The spectra were then processed with the CERES pipeline (Brahm et al., 2017) to obtain both RVs and activity indicators. The measurements are given in Table 3, and shown in Figure 13

The MINiature Exoplanet Radial Velocity Array (Minerva)-Australis is an observatory located in Queensland, Australia, dedicated to the precise radial-velocity and photometric follow-up of TESS planet candidates (e.g. Addison et al., 2019, 2021b; Wittenmyer et al., 2022) It consists of four 0.7 m robotic telescopes fiber-fed to a KiwiSpec spectrograph with spectral resolution of $\sim$80,000 and wavelength coverage between 480 and 620 nm (Wilson et al., 2019).

We obtained 31 spectra of HD 18599 between 2019 Jan 6 and Jan 29 with a typical exposure time of 20-30 minutes. Radial velocities for the observations are derived for each telescope by cross-correlation, where the template being matched is the mean spectrum of each telescope. The instrumental variations are corrected by using simultaneous Thorium-Argon arc lamp observations. The measurements are given in Table 3, and shown in Figure 13.

We conducted high resolution spectroscopy of HD 18599 using CTIO High Resolution spectrometer (CHIRON) on the 1.5-m SMARTS telescope. It has a spectral resolution of 80,000 with wavelength coverage between 4500 and 8900 Å. Between 2019 Feb and 2020 Dec, we took a total of 6 spectra with typical SNR between 53 and 130.

We conducted high resolution spectroscopy of HD 18599 using Las Cumbres Observatory’s (LCO; Brown et al., 2013) Network of Robotic Echelle Spectrographs (NRES; Siverd et al., 2018). It has a spectral resolution of 53,000 with wavelength coverage between 3800 and 8600 Å. The observations were conducted on several nights with exposure times between 480s and 1200s, resulting in one observation at the unit at the Cerro-Tololo Inter-American Observatory (CTIO) in Chile and two observations at the South African Astronomical Observatory (SAAO). All three observations consisted of three consecutive exposures which were binned on a nightly basis. Due to weather conditions, only two out of the three nightly binned spectra had enough SNR for a confident spectral classification using the SpecMatch-Synth code⁷⁷7https://github.com/petigura/specmatch-syn. The spectrum obtained at the LCO CTIO node has SNR=102 and the spectrum obtained at the LCO SAAO node has SNR=45. We note that the stellar parameters from the spectral classification are in agreement with the ones from the isochrone estimates.

To obtain the physical properties of HD 18599, we utilized the Python package isochrones (Morton, 2015a)⁹⁹9http://github.com/timothydmorton/isochrones; v.2.1 that relies on the MESA Isochrones & Stellar Tracks (MIST; Dotter 2016) grid to infer stellar parameters using a nested sampling scheme given photometric or spectroscopic data and other empirical constraints. In particular, we used 2MASS ($JHKs$) (Skrutskie et al., 2006) and along with Gaia DR2 parallax (Gaia Collaboration et al., 2018) and extinction. We corrected the parallax for the offset found in Stassun & Torres (2018) while quadratically adding 0.1 mas to the uncertainty to account for systematics in the Gaia DR2 data (Luri et al., 2018). Additionally, we used $T_{\mathrm{eff}}$, $\log g$, and $[\mbox{Fe}/\mbox{H}]$ derived from spectroscopy (see Section 2.2) or taken from the literature as additional priors. The results of isochrones are summarized in Table 4.

This section presents an independent method to derive empirical stellar parameters of HD 18599 which also serves to cross-check our results obtained from the isochrones method (Section 3.1.1). Following the procedures described in Stassun & Torres (2016); Stassun et al. (2017); Stassun et al. (2018), we first construct the broadband spectral energy distribution (SED) of HD 18599, and then derive the stellar radius using the Gaia DR2 parallax. The broadband photometry measurements include the FUV and NUV magnitudes from GALEX, the $B_{T}V_{T}$ magnitudes from Tycho-2, the $JHK_{S}$ magnitudes from 2MASS, the W1–W4 magnitudes from WISE, and the $GG_{\rm RP}G_{\rm BP}$ magnitudes from Gaia. Altogether, the available photometry spans the full stellar SED over the wavelength range 0.2–22 $\mu m$. Then, we fit the SED with the NextGen stellar atmosphere models taking into account the effective temperature ($T_{\rm eff}$) and metallicity ([Fe/H]) derived from the spectroscopic analysis. The extinction ($A_{V}$) was set to zero because of the star being very near at d=40 pc.

Figure 6 shows the SED plot with broadband photometry measurements superposed with the best-fit model. The resulting fit is excellent with a reduced $\chi^{2}$ of 1.9 (excluding the GALEX FUV and NUV fluxes, which are consistent with moderate chromospheric activity). Integrating the unreddened SED model yields the bolometric flux at Earth of $F_{\rm bol}=8.013\pm 0.093\times 10^{-9}$ erg s^-1 cm^-2. Taking the $F_{\rm bol}$ and $T_{\rm eff}$ together with the Gaia DR2 parallax, which was adjusted to account for the systematic offset reported by Stassun & Torres (2018), yields the stellar radius as $R_{\star}=0.781\pm 0.016$ R_⊙. In addition, the stellar bolometric luminosity is obtained directly from $F_{\rm bol}$ and the parallax, yields $L_{\rm bol}=0.3709\pm 0.0044$ L_⊙. Finally, we estimate the stellar mass from the empirical relations of Torres (2010) and a 6% error from the empirical relation itself yields $M_{\star}=0.87\pm 0.05$ $M_{\odot}$, whereas the mass estimated empirically from $R_{\star}$ together with the spectroscopic $\log g$ yields $M=0.56\pm 0.13$ M_⊙. This discrepancy suggests that the spectroscopic $\log g$ may be slightly underestimated. In any case, these values are in agreement with the values estimated using isochrones method in Section 3.1.1 which uses high-resolution spectra. The final values are listed in Table 4.

HD 18599 is previously known to be a young field star based on previous studies (Grandjean et al., 2020, 2021). Thus, it comes to no surprise that we did not find HD 18599 to be a member of any of the 1641 clusters catalogued in Cantat-Gaudin & Anders (2020) as well as among the young moving groups identified in (Gagné & Faherty, 2018; Gagné et al., 2018). We double checked for co-moving stars in the neighborhood by querying the kinematics of all sources within 10’ of the target from the Gaia DR2 catalog and found no match as found using the comove ¹⁰¹⁰10https://github.com/adamkraus/Comove code. We found a match in TGv8 catalog (Carrillo et al., 2020) using 4728513943538448512 and found the probability of HD 18599 being in the galactic thin and thick disk to be 98.5% and 1.4%, respectively. This is not surprising given HD 18599’s high galactic latitude of $\sim$53 deg.

The strength of the lithium 6708Å absorption feature has traditionally been used as a youth indicator for Sun-like stars (Soderblom et al., 2014). When mixed down into a layer in a star that is roughly $2.5\times 10^{6}$ K, Li is burned into heavier elements, and thus rapidly destroyed within the first $\sim$500 Myr. Thus, measuring the lithium 6708Å absorption feature is a strong evidence that the star’s age is less 1 Gyr. We inspected Li I $\lambda$6708 in FEROS and CHIRON spectra to confirm the youth of HD 18599. We measured an equivalent width (EW) of 42$\pm$2mÅ where the uncertainty comes from the standard deviation of EW measurement for each spectra. This Li strength generally agrees with those found in stars in Hyades and Praesepe corresponding to an age of 800 Myr. The EW width is narrower than those measured by Bowler et al. (2019) for 58 Li-rich K to M stars (100–600 mÅ) in the solar neighborhood with ages between 10-120 Myr.

We can also estimate the stellar age by taking advantage of the observed chromospheric activity together with empirical age-activity-rotation relations. For example, taking the chromospheric activity indicator, $\log R^{\prime}_{HK}=-4.41\pm 0.02$ from Boro Saikia et al. (2018) and applying the empirical relations of Mamajek & Hillenbrand (2008), gives a predicted age of $0.30\pm 0.05$ Gyr. Whilst this star is an X-ray source based on detection from ROSAT, the X-ray strength is weak ($\log Lx/Lbol=-4.64\pm 0.25$) which corresponds to $1-\sigma$ age range from X-ray of $475^{+734}_{-305}$ Myr. The X-ray count to luminosity calibration is from Fleming et al. (1995), and the X-ray luminosity to age calibration is from Equation A3 in Mamajek & Hillenbrand (2008).

The top panel in Figure 2 clearly shows significant spot-modulated rotational signals in the TESS light curves. To measure the rotation period of HD 18599, we used the light curves from TESS PDCSAP and KELT Oelkers et al. (2018). The top panel in Figure 9 shows the Lomb-Scargle (LS) periodogram of the four sectors of TESS (red) and the 6-year KELT light curves (black). Both show a consistent peak at $P_{\rm{rot}}$ =8.8 d. The bottom panel in Figure 9 shows the spot evolution of HD 18599 over each of four TESS sectors. We note that the light curves from sectors 2 and 3 showed a strong secondary peak at 1/2 rotation period in the LS periodogram, with the corresponding phase folded light curve showing that the star likely has large spot regions on both hemispheres. By sectors 29 and 30, these spots had evolved such that the 8.71-day period is the dominant frequency in the periodogram, which we also confirmed from KELT light curves. We note also that the stellar rotation period does not coincide with the planetary orbital period.

We can further corroborate the activity-based age estimate by also using empirical relations to predict the stellar rotation period from the activity. For example, the empirical relation between $R^{\prime}_{HK}$ and rotation period from Mamajek & Hillenbrand (2008) predicts a rotation period for this star of $9.7\pm 1.3$ d, which is compatible with the rotation periods above as well as with the $P_{\rm{rot}}$ of 8.69 d reported by KELT, and also compatible with the projected rotation period inferred from the spectroscopic $v\sin i$ and $R_{\star}$ which gives $P_{\rm rot}/\sin i=7.2\pm 3.3$ d. Using $P_{\rm{rot}}$ =8.69 d, $R_{\star}$ =0.77 $R_{\odot}$, $v\sin i$ =4.3 km/s, we derive an inclination of 74 deg.

Using the Generalized Lomb-Scargle (GLS) periodogram (Zechmeister & Kürster, 2009) on TESS light curves, we derived $P_{\rm{rot}}$ and amplitude of 8.582549$\pm$0.003679 d and 0.8673$\pm$0.0014%, respectively. Using the age model as a function of "smoothed" rotation amplitude presented in Morris (2020), we derived a poorly constrained age of 1145${}^{+1718}_{-960}$ Myr.

Finally, we can also estimate the age from the observed $P_{\rm{rot}}$ and empirical gyrochronology relations Mamajek & Hillenbrand (2008), we derive a median age of 386 Myr, with 3$\sigma$ range of 261-589 Myr. Using Barnes (2007), we derive a median age of 247 Myr, with 3$\sigma$ range of 185-329 Myr.

Figure 11 shows the summary of the age estimates using the different methods consistent with the age between 0.1 to 1 Gyr. Assuming all the above methods are independent, we compute a weighted mean of 273$\pm$22 Myr. We adopt a median age of 300 Myr for the discussion in Section 5.

A number of astrophysical scenarios can mimic the transit signal detected from TESS photometry, including eclipsing binary (EB) with a grazing geometry, hierarchical EB (HEB), and background/foreground EB (BEB) along the line of sight of the target. Firstly, the scenario that HD 18599 is actually a grazing EB is ruled out based on the small (<1km s^-1) RV variations which disfavors stellar mass companions (Section 2.2). The TESS light curves also do not exhibit secondary eclipses and odd-even transit depth variations. The light curve is also flat-bottomed as compared to V-shaped for EBs. Secondly, the scenario that HD 18599 is actually an NEB is ruled out due to the lack of Gaia sources found within the TESS and Spitzer photometric apertures that are bright enough to reproduce the TESS detection. Thirdly, the scenario that HD 18599 is actually an HEB with component stars of different colors is ruled out based on the achromatic transits from TESS and Spitzer (Section 4.1).

The most plausible HEB scenarios for HD 18599 involve pairs of eclipsing M dwarfs. Eclipses of such stars are deeper than the K-dwarf HD 18599 in longer wavelengths. Limits on whether the transit depth decreases in shorter wavelengths can therefore rule out certain HEB scenarios. Similar to the procedure described in Bouma et al. (2020), we fitted for the observed depths in different bandpasses. To perform the calculation, we assumed that each system was composed of the primary star (HD 18599, Star 1), plus a tertiary companion (Star 3) eclipsing a secondary companion (Star 2) every $P_{\mathrm{orb}}$ d. For a grid of secondary and tertiary star masses ranging from 0.07¹⁵¹⁵15This mass is comparable to the lowest star mass known (von Boetticher et al., 2017). to 0.5$M_{\odot}$, we then calculated the observed maximum eclipse depth caused by Star 3 eclipsing Star 2 in TESS and Spitzer bandpasses using the following procedure. First, we interpolated $L_{\star}$ and $T_{\mathrm{eff}}$ of Star 2 and Star 3 from MIST isochrones given their masses, and the age, metallicity, and mass of Star 1 in Table 4. We then computed the blackbody functions of each stars given their $T_{\mathrm{eff}}$ then convolved it with the transmission functions for each band downloaded from the SVO filter profile service¹⁶¹⁶16http://svo2.cab.inta-csic.es/theory/fps/. We then integrated the result using trapezoidal method and computed the bolometric flux, $F_{\rm bol}$, using the integrated functions above. Using Stefan-Boltzman law and given $T_{\mathrm{eff}}$ and $L_{\star}$, we computed the component radii and luminosity to derive the eclipse depth.

For an HEB system with identical component stars, these stars should be very small ($m_{1}=m_{2}=$0.1$M_{\odot}$) taking into account dilution from the central star to reproduce the observed TESS depth ($\sim$1 ppt). In this contrived scenario, the eclipse depth in Spitzer I2 is about 8 times deeper than in TESS band, because the early M-dwarf blackbody function turns over at much redder wavelengths than the central K star blackbody (Wien’s law). Thus, there is no plausible HEB configuration explored in our simulation above that can reproduce the observed depth in Spitzer I2. Even in the extreme case of a grazing orbit (i.e. $b=0.99$) for the tertiary companion, the resulting decrease in eclipse depth in Spitzer I2 is still twice as large as the observed depth. Altogether, the multiwavelength depth constraint rules out the HEB scenario. Note that the "boxy" shape of the transit signal especially in Spitzer I2 can also help rule out most of the parameter space of the HEB configurations considered above which generally produce V-shaped eclipses.

Finally, the scenario that HD 18599 is actually a BEB is negligibly small and can be completely ruled out. To quantify the probability of chance alignment of a BEB to HD 18599, we use the population synthesis code TRILEGAL¹⁷¹⁷17http://stev.oapd.inaf.it/cgi-bin/trilegal (Girardi et al., 2005), which simulates stellar parameters in any Galactic field. We found there is a 25% chance to find an EB brighter than Tmag=15, within an area equal to the TESS photometric aperture (24 TESS pixels $\approx$ 2.94 arcmin²). We can eliminate the presence of any stellar companion with $\Delta M$<5 up exterior to 0.1" using our high spatial resolution speckle images (Section 2.4)¹⁸¹⁸18Interior to 0.1” however, we could not clearly show that there is no bright enough star in the target’s current position based on the earliest archival image taken during POSS-1 sky survey $\sim$60 years ago when the target would have moved 0.94” relative to the background stars. The brightness of the target essentially masked the view along the line of sight to the target in the archival image. This result is consistent with the observed paucity of close binaries in TESS host stars (Ziegler et al., 2021; Lester et al., 2021). Within this area, we found that there is 7.66$\times 10^{-5}$ chance to find a chance-aligned star. Note that this is a conservative upper limit because this result assumes all stars are binary and preferentially oriented edge-on to produce eclipses consistent with TESS detection.

Although our spectroscopic analysis disfavors the existence of an BEB spectroscopically blended with HD 18599, we cannot rule out all false positive scenarios involving HEB with identical color components. For this reason, we use VESPA to model the relevant EB populations statistically.

We quantify the false positive probability (FPP) of HD 18599.01 using the Python package VESPA ¹⁹¹⁹19https://github.com/timothydmorton/VESPA which was developed as a tool for robust statistical validation of planet candidates identified by the Kepler mission (e.g. Morton, 2012) and its successor K2 (e.g. Crossfield et al., 2016; Livingston et al., 2018b; Mayo et al., 2018). In brief, VESPA compares the likelihood of planetary scenario to the likelihoods of several astrophysical false positive scenarios involving eclipsing binaries (EBs), hierarchical triple systems (HEBs), background eclipsing binaries (BEBs), and the double-period cases of all these scenarios. The likelihoods and priors for each scenario are based on the shape of the transit signal, the star’s location in the Galaxy, and single-, binary-, and triple-star model fits to the observed photometric and spectroscopic properties of the star generated using isochrones.

As additional constraints, we used the available speckle contrast curves described in Section 2.4, the maximum aperture radius (maxrad)–interior to which the transit signal must be produced, and the maximum allowed depth of potential secondary eclipse (secthresh) estimated from the given light curves. Similar to Mayo et al. (2018), we computed secthresh by binning the phase-folded light curves by the measured the transit duration and taking thrice the value of the standard deviation of the mean in each bin. Effectively, we are asserting that we did not detect a secondary eclipse at any phase (not only at phase=0.5) at 3-$\sigma$ level. Given these inputs, we computed a formal FPP=1.48e-11 which robustly qualifies HD 18599.01 as a statistically validated planet.

Additionally, we validated HD 18599.01 using the Python package TRICERATOPS ²⁰²⁰20https://github.com/stevengiacalone/triceratops which is a tool developed to validate TOIs (Giacalone & Dressing, 2020; Giacalone et al., 2021). TRICERATOPS validates TOIs by calculating the Bayesian probabilities of the observed transit originating from several scenarios involving the target star, nearby resolved stars, and hypothetical unresolved stars in the immediate vicinity of the target. These probabilities were then compared to calculate a false positive probability (FPP; the total probability of the transit originating from something other than a planet around target star) and a nearby false positive probability (NFPP; the total probability of the transit originating from a nearby resolved star). As an additional constraint in the analysis, we folded in the speckle imaging follow-up obtained with the Zorro imager on Gemini-South. For the sake of reliability, we performed the calculation 20 times for the planet candidate and found FPP=0.0027$\pm$0.0008 and NFPP=0 (indicating that no nearby resolved stars were found to be capable of producing the observed transit). Giacalone et al. (2021) noted that TOIs with FPP < 0.015 and NFPP < 10^-3 have a high enough probability of being bona fide planets to be considered validated. Our analysis using TRICERATOPS therefore further added evidence to the planetary nature of HD 18599.01. We now refer to the planet as HD 18599 b in the remaining sections.