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

Transiting exoplanets orbiting young nearby stars are ideal laboratories for testing theories of planet formation and evolution. However, to date only a handful of stars with age<1 Gyr have been found to host transiting exoplanets. Here we present the discovery and validation of a sub-Neptune around HD 18599, a young (300 Myr), nearby (d=40 pc) K star. We validate the transiting planet candidate as a bona fide planet using data from the TESS, Spitzer, and Gaia missions, ground-based photometry from IRSF, LCO, PEST, and NGTS, speckle imaging from Gemini, and spectroscopy from CHIRON, NRES, FEROS, and Minerva-Australis. The planet has an orbital period of 4.13 d, and a radius of 2.7Rearth. The RV data yields a 3-sigma mass upper limit of 30.5Mearth which is explained by either a massive companion or the large observed jitter typical for a young star. The brightness of the host star (V~9 mag) makes it conducive to detailed characterization via Doppler mass measurement which will provide a rare view into the interior structure of young planets.


INTRODUCTION
The majority of currently known exoplanets orbit mature host stars.In order to have a complete understanding of how planetary systems evolve from birth to maturity, we need to have a sample of systems at various evolutionary stages.In particular, young exoplanets (<1 Gyr in age) inhabit a very important part of the exoplanet evolution-★ E-mail: jpdeleon@g.ecc.u-tokyo.ac.jp ary timescale, where formation mechanisms, accretion, migration and dynamical interactions can significantly change the shape of observed planetary systems.This is well illustrated by our knowledge of the planetary system for which we have by far the most information -the Solar system.When we study the Solar system's history, it is immediately apparent that the first few hundred million years were an extremely chaotic and volatile time, with the system evolving dramatically on short timescales.Evidence abounds of the giant impacts that shaped both the terrestrial and giant planets, including the big smash which led to Mercury's reduced size and enhanced density (e.g.Benz et al. 1988Benz et al. , 2007;;Chau et al. 2018), the collision that shattered the proto-Earth and formed the Moon (e.g.Benz et al. 1989;Canup 2012;Canup & Asphaug 2001), the potential origin of Mars' hemispheric anomaly (e.g.Andrews-Hanna et al. 2008), the formation of Jupiter's diluted core (e.g.Liu et al. 2019), and Uranus' peculiar axial tilt (e.g.Slattery et al. 1992;Parisi & Brunini 1997;Kegerreis et al. 2018).At the same time, the system's small body populations reveal the scale of the migration of the giant planets.In the case of Uranus and Neptune, those small bodies reveal a significant outward migration 1 (e.g.Lykawka et al. 2010Lykawka et al. , 2011)).In the case of Jupiter and Saturn, the evidence points to significant inward migration, and potentially a period of significant dynamical instability -though the scale, rate, and chaoticity of that migration remains the subject of much debate (e.g. the Nice model: Tsiganis et al. (2005), the Grand Tack: O'Brien et al. (2014) and Walsh et al. (2012); or a smoother, less chaotic migration: Lykawka & Horner (2010) and Pirani et al. (2019).For a detailed overview of our knowledge of the Solar system, tailored towards Exoplanetary Science, we direct the interested reader to Horner et al. (2020).
To date, only a handful of exoplanet host stars have a well constrained age.The youngest known transiting systems so far include K2-33 (age=5-10 Myr, David et al. 2016;Mann et al. 2016), HIP 67522 (17 Myr, Rizzuto et al. 2020), TAP 26 (17 Myr, Yu et al. 2017), Au Mic (23 Myr Plavchan et al. 2020), V1298 Tau (23 Myr, David et al. 2019), DS Tuc A (45 Mr, Newton et al. 2019), TOI 837 (45 Myr, Bouma et al. 2020), andTOI 942 (50 Myr, Zhou et al. 2021;Carleo et al. 2020;Wirth et al. 2021).The rarity of known young planets is attributed to the difficulties involved in their detection-the large intrinsic stellar activity in young stars that induces large photometric variations and radial velocity (RV) jitter often a few orders of magnitude larger than the planet signal (e.g.Heitzmann et al. 2021;Nicholson et al. 2021).Recently, novel methods have been developed to overcome this problem by using a method to detrend both photometry and RV data (e.g.Dai et al. 2017;Collier Cameron et al. 2021).Despite the small sample of planets detected so far, there is tentative evidence that there is a measurable change in the occurrence rates of planets with time (See Figure 14; Mann et al. 2016), although more samples are needed to robustly confirm this trend.More recently, Berger et al. (2020) and Sandoval et al. (2020) found that the ratio of super-Earth to sub-Neptune detections in the California-Kepler Survey (CKS) sample increases with system age between 1-10 Gyr even without accounting for completeness effect.Relative to sub-Neptunes, super-Earths appear to be more common around older stars despite the difficulty of detecting small planets around larger stars.Moreover, David et al. (2020) found that the size distribution of small planets also depends on the age of the planet population and that the precise location of the radius valley changes over gigayear timescales.However, we are still in the early stages to confidently corroborate or refute such trends given the low number of transiting planets orbiting young stars.Therefore, by compiling a statistically significant sample of well-characterized exoplanets with precisely measured ages, we should be able to begin identifying new trends as well as other dominant processes governing the time-evolution of exoplanet systems.
Here we present the discovery and validation of a sub-Neptune around HD 18599 also known as TOI 179, a young (300 Myr), 1 in the case of Neptune, travelling over at least 10 au nearby (d=38.6 pc) K star.In Section 2, we discuss the observations including the discovery data from TESS, follow-up photometry from space by Spitzer and from the ground by several telescopes, archival photometry and kinematics from Gaia, spectroscopy from M --Australis, FEROS, CHIRON, and NRES, and speckle imaging from Gemini.In Section 3, we present our analyses to characterize the host star and in particular establish its youth.In Section 4, we derive the properties of the planet such as its precise radius and mass limit.We also synthesize all available data to validate the planet candidate.In Section 5, we put the planet in context with the known young population.Finally, we summarize our results in Section 6 and also motivate further follow-up efforts of HD 18599 to measure the planet's mass and potentially characterize its interior and atmosphere.

TESS discovery
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 ∼4-day gap because of the data downlink and telescope re-pointing.
The TESS object of interest (TOI) releases portal2 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)3 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-squares4 (Hippke & Heller 2019) after applying a biweight filter using wotan5 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.

Spitzer follow-up
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 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 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.

IRSF follow-up
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 , ,   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.

LCO follow-up
We observed HD 18599 on 2021 Sep 18 with 1,0m LCO-CTIO in the B and   bands, on 2018 Nov 18 with 1.0m LCO-CPT in the  band, on 2018 Nov 23 with the 1.0m LCO-CTIO in the  band, and on 2018 Dec 22 with the 0.4m LCO-CTIO in the  band.All observations were full transit except that one taken in  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.

NGTS follow-up
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×2K e2V deep-depleted Andor Ikon-L CCD cameras with 13.5m 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).

PEST followup
We observed HD 18599 in the   band from the Perth Exoplanet Survey Telescope (PEST) near Perth, Australia.The 0.3 m telescope is equipped with a 1530×1020 SBIG ST-8XME camera with an image scale of 1.2" pixel −1 resulting in a 31'×21' field of view.A custom pipeline based on C-Munipack6 was used to calibrate the images and extract the differential photometry.

KELT archival data
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×4k Apogee CCD, a pixel scale of 23 and a field of view of 26 • × 26 • .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).

Spectroscopy
We conducted several high resolution spectroscopic observations to characterize the host star and to measure the RV variation induced by orbiting companion.In the following, we describe our observations first and then derive an upper mass limit of the companion.

FEROS
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.2mtelescope 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

M -Australis
The MINiature Exoplanet Radial Velocity Array (M )-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. 2019Addison et al. , 2021b;;Wittenmyer et al. 2022) It consists of four 0.7 m robotic telescopes fiber-fed to a KiwiSpec spectrograph with spectral resolution of ∼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 crosscorrelation, 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.

CTIO/CHIRON
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.

LCO/NRES
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 clas- sification using the SpecMatch-Synth code7 .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.

Gaia astrometry
Between 25 July 2014 and 23 May 2016, the ESA Gaia satellite measured about 300 billion centroid positions of 1.6 billion stars.The positions, proper motions, and parallaxes of the brightest 1.3 billion sources were calculated for the second data release (DR2) (Gaia Collaboration et al. 2018).HD 18599 was assigned the Gaia DR2 identifier 4728513943538448512, and had 265 "good" astrometric observations, indicating a nearby (d=38.6 pc), high-proper motion (-36.68, 50.60 mas/yr) star.
We further leverage Gaia DR2 to search for direct and indirect evidence of potential contaminating sources.In our sample, we can probe Gaia DR2 sources separated from the target as close as 1 .Gaia DR2 can also be useful to look for hints of binarity.Evans (2018) proposed that systems with large Astrometric Goodness of Fit of the astrometric solution for the source in the Along-Scan direction (GOF_AL > 20)) and Astrometric Excess Noise significance (D > 5)8 are plausibly poorly-resolved binaries.Stars that are exceptionally bright or have high proper motion are proposed to explain the large offset due causing difficulties in modelling saturated or fast-moving stars, rather than unresolved binarity.We found GOF_AL = 7.0 and D = 0 for HD 18599 which are well below the aforementioned empirically-motivated cutoffs, indicating the target is indeed single.

Gemini speckle imaging
The presence of multiple unresolved stars in photometric and spectroscopic observations of a transiting planetary system biases measurements of the planet's radius, mass, and atmospheric conditions (e.g.Furlan & Howell 2020;Southworth & Evans 2016).To determine if any fainter point sources existed closer to target inside of Gaia's point-source detection limits and to rule out false positives caused by an eclipsing binary as well as to search for potential (sub-)stellar companions within a few arcseconds from the target, we conducted speckle imaging using the Zorro instrument in b and r cameras centered on 562nm and 832nm, respectively, mounted on the 8.1m Gemini South Telescope at the Cerro Pachon, Chile (Scott et al. 2021).Smooth contrast curves were produced from the reconstructed images by fitting a cubic spline to the 5- sensitivity limits within a series of concentric annuli.The speckle observations with their corresponding contrast curves in Figure 5 illustrate that no companions were detected within a radius of 1.2 down to a contrast level of 7 magnitudes, and no close companion star was detected within angular resolutions of the diffraction limit (0.02 ) out to 1.2 .At the distance of HD 18599 (d=38.6 pc), these limits correspond to spatial limits of 0.8 to 46 AU.These observations sharply reduce the possibility that an unresolved background star is the source of the transits.The contrast curves are also used as additional constraints for false positive calculation in Section 4.3.

HOST STAR PROPERTIES
HD 18599 (HIP 13754, TOI 179, TIC 207141131) is a known young, nearby K2V star included in the SPHERE GTO The SpHere INfrared survey for Exoplanets (SHINE) survey sample (Chauvin et al. 2017).Based on 16 HARPS spectroscopic observations, Grandjean et al. (2020Grandjean et al. ( , 2021) ) did not detect any planetary companion.They also measured a systemic RV=115.9±38.9m/s (rms) and chromospheric flux ratio of log    =-4.310 indicating a high activity index mainly due to stellar starspots.In the following, we characterize the host star in detail.

isochrones
To obtain the physical properties of HD 18599, we utilized the Python package isochrones (Morton 2015a) 9 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 ( ) (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  eff , log , and [Fe/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.

Spectral Energy Distribution
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. (2017Stassun 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     magnitudes from Tycho-2, the   magnitudes from 2MASS, the W1-W4 magnitudes from WISE, and the  RP  BP magnitudes from Gaia.Altogether, the available photometry spans the full stellar SED over the wavelength range 0.2-22 .Then, we fit the SED with the NextGen stellar atmosphere models taking into account the effective temperature ( eff ) and metallicity ([Fe/H]) derived from the spectroscopic analysis.The extinction (  ) 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  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  bol = 8.013 ± 0.093 × 10 −9 erg s −1 cm −2 .Taking the  bol and  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  ★ = 0.781 ± 0.016 R .In addition, the stellar bolometric luminosity is obtained directly from  bol and the parallax, yields  bol = 0.3709 ± 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  ★ = 0.87 ± 0.05  , whereas the mass estimated empirically from  ★ together with the spectroscopic log  yields  = 0.56 ± 0.13 M .This discrepancy suggests that the spectroscopic log  may be slightly underestimated.In any case, these values are in agreement with the values estimated using isochrones method in Section 3.1.1which uses high-resolution spectra.The final values are listed in Table 4. Soderblom et al. (2014) provides a comprehensive review of the techniques to determine approximate stellar ages.To age-date HD 18599, first we search for coeval, phase-space neighbors and compile a sample of candidate siblings to compare with the empirical sequences of young clusters (color-magnitude diagram,    , Li I 6708 equivalent width, stellar rotation period) and gyrochonology relationships.Age dating stars is notoriously difficult except for ensemble of stars.Hence, most of the currently known young stars hosting transiting planets are found in stellar associations where age can be reliably measured.In fact, the youngest transiting host star known so far, K2-33, is found in a 5-10 Myr moving group in Upper Sco.According to MESA grids, the Pre-Main Sequence (PMS) contraction time of a  ★ =0.8 star is ∼67 Myr.Thus, the parameters of the star do not significantly change between 0.1 to 1 Gyr.However, we can look at various youth indicators to derive the age of HD 18599.

Stellar association
HD 18599 is previously known to be a young field star based on previous studies (Grandjean et al. 2020(Grandjean et al. , 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 comove10 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 ∼53 deg.

Li I 6708 Å equivalent width
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 × 10 6 K, Li is burned into heavier elements, and thus rapidly destroyed within the first ∼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 6708 in FEROS and CHIRON spectra to confirm the youth of HD 18599.We measured an equivalent width (EW) of 42±2mÅ 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.

Activity indicators
We can also estimate the stellar age by taking advantage of the observed chromospheric activity together with empirical age-activityrotation relations.For example, taking the chromospheric activity indicator, log    = −4.41± 0.02 from Boro Saikia et al. ( 2018) and applying the empirical relations of Mamajek & Hillenbrand (2008), gives a predicted age of 0.30 ± 0.05 Gyr.Whilst this star is an X-ray source based on detection from ROSAT, the X-ray strength is weak (log / = −4.64 ± 0.25) which corresponds to 1 −  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).

Stellar rotation period and amplitude
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  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

Gyrochronology
Finally, we can also estimate the age from the observed  rot and empirical gyrochronology relations Mamajek & Hillenbrand (2008), we derive a median age of 386 Myr, with 3 range of 261-589 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±22 Myr.We adopt a median age of 300 Myr for the discussion in Section 5.

Transit modeling
In Section 2.1.2,we confirmed that the TESS signal indeed originates from HD 18599.In the following, we model the TESS and Spitzer light curves jointly to robustly measure the the planet's parameters, in particular the transit depth as a function of bandpass.
After removing all flagged cadences and those that are more than 3- above the running mean, we flattened and normalized the raw light curves using a median filter with kernel size of 301 cadences corresponding to ∼5× the transit duration.After the preprocessing steps, we model the TESS and Spitzer light curves using the Python package allesfitter11 detailed in Günther & Daylan (2020).In brief, allesfitter is a tool developed for joint modeling of photometric and RV data with flexible systematics models and was extensively used in related exoplanet studies with TESS (e.g.Huang et al. 2018;Dragomir et al. 2019;Addison et al. 2021a).Parameter estimation is done with nested sampling using dynesty (Speagle 2020) of transit and RV models defined in ellc (Maxted 2016) and Gaussian Processes for systematics models using celerite (Foreman-Mackey et al. 2017).Our motivation to use allesfitter is its flexible capability to model the transit and RV data either separately or jointly including model comparison.Nested sampling is also more efficient than MCMC especially for complex models that require more than 10 model parameters.
We set the following as free parameters: the orbital period  orb , mid-transit time  0 , scaled semi-major axis /  , (cosine of) inclination , and quadratic limb darkening coefficients in q-space ( 1 and  2 ) as prescribed by Kipping (2013).We also fit for the logarithm of the Gaussian flux errors (log ), and the logarithm of the two hyper-parameters of the Gaussian process (GP) model with an approximated Matérn-3/2 kernel for the out-of-transit baseline: where  and  reflects the characteristic amplitude and length scale of the GP, respectively.We also tried other kernels, such as simple harmonic oscillator, but Matérn-3/2 kernel yielded the highest Bayesian evidence.
We imposed uniform priors on the host star's density computed from the stellar parameters reported in Table 4, taking into account uncertainties.This ensures that at each step in the Nested Sampling run, the stellar density is derived from the planet's orbit (/  and  orb ), and the fit gets penalized for deviation from the prior.We also imposed (truncated) Gaussian priors on  1 and  2 , with mean and standard deviation computed using limbdark12 , which is a tool for Monte Carlo sampling an interpolated grid of the theoretical limb darkening coefficients (in the TESS and Spitzer bandpasses) tabulated by Claret et al. (2012) given  eff , [Fe/H], and log  of the host stars.We ran the nested sampler in dynamic mode with 500 live points until it reached the prescribed 1% tolerance criterion for convergence.Based on 35611 posterior samples, we report the median and 68% credible interval of the resulting marginalized posterior distributions in Table 5.The TESS and Spitzer light curves with best-fit transit model is shown in Figure 2 and Figure 3.The comparison of the marginalized posterior distributions of the said parameters from modeling using TESS-only, Spitzer-only, and jointly are shown in Fig. 12.The joint transit modeling of the the TESS and Spitzer lightcurves resulted to tighter constraints on the radius ratio and impact parameter.

Companion mass constraint
To put a limit on the mass of putative HD 18599 b, we fit an RV model with a circular orbit to the RV data from M -Australis and FEROS based on their small RV scatter, more precise than the rest of our RV data.We used the RV model included in the PyTransit python package13 which we simplified to have five free parameters: phase-zero epoch T 0 , period, RV semi-amplitude, RV zero point, and RV jitter term.For the T 0 and the period, we put Gaussian priors using the T 0 and period derived from the transit analysis.For the other parameters we put wide uniform priors.We ran the built-in Differential Evolution optimizer and then sample TESS-only (blue), Spitzer-only (red), and joint (black) dataset.The higher radius ratio obtained from the joint fit relative to the single band fits is due to its larger and more precise impact parameter.the parameters with Markov Chain Monte Carlo (MCMC) using 30 walkers and 10 4 steps.We use the following equation from Cumming et al. (1999) to derive the planet mass, where   is planet mass,  ★ is star mass,  is orbital period,  is RV semi-amplitude,  is eccentricity (fixed to zero), and  is inclination (fixed to 90 • ).To propagate uncertainties, we use the posteriors for  ★ and  from previous analyses.
In Figure 13, we plot Keplerian orbital models correspond- ing to different masses encompassing the 68 th , 95 th , and 99.7 th percentiles of the semi-amplitude posterior distribution.The 3- upper limit is 30.5  ⊕ which places the companion 2 orders of magnitude below the deuterium burning mass limit.The bestfit semi-amplitude is 7.42  ⊕ , which corresponds to a mass of M p =7.42  ⊕ , and the best-fit jitter values are   , =28.3 m s −1 and   , =12.4 m s −1 for M -Australis and FEROS, respectively.
We calculated an expected planetary mass of 6.66 +15.84 −2.82  ⊕ with MRExo14 , which uses a mass-radius relationship calibrated for planets around Kepler stars (Ning et al. 2018).This mass corresponds to a semi-amplitude of 2.99 m s −1 , but the observed RV data exhibits significantly larger variability (  =27m s −1 for M -Australis and   =14m s −1 for FEROS) as expected for young stars.We interpret this variability as being responsible for the large jitter value found by the fit, which suggests it is out-ofphase with HD 18599.01.Grandjean et al. (2021) also identified activity as the source of the stellar jitter.We note that a detailed RV mass measurement of the companion will be presented in an accompanying paper (Vines et al., in prep.).

False positive analysis
In this section, we present the validation of the planet candidate HD 18599.01.First, we describe each constraint to eliminate false positive scenarios constrained by our data.Next, we present an FPP calculation using VESPA (Morton 2015b) and TRICERATOPS (Giacalone & Dressing 2020) to demonstrate that the probability that HD 18599.01 is an astrophysical false positive is small enough to formally validate it as a planet.

Eliminating false positive scenarios
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  orb d.For a grid of secondary and tertiary star masses ranging from 0.0715 to 0.5 , 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  ★ and  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  eff then convolved it with the transmission functions for each band downloaded from the SVO filter profile service 16 .We then integrated the result using trapezoidal method and computed the bolometric flux,  bol , using the integrated functions above.Using Stefan-Boltzman law and given  eff and  ★ , 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 ( 1 =  2 =0.1 ) taking into account dilution from the central star to reproduce the observed TESS depth (∼1 ppt).In this contrived scenario, the eclipse depth in Spitzer I2 is about 8 times deeper than in TESS band, because the early Mdwarf 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. = 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 TRILEGAL17 (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 ≈ 2.94 arcmin 2 ).We can eliminate the presence of any stellar companion with Δ<5 up exterior to 0.1" using our high spatial resolution speckle images (Section 2.4) 18 .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×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.

Statistical validation with VESPA and TRICERATOPS
We quantify the false positive probability (FPP) of HD 18599.01 using the Python package VESPA19 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- level.Given these inputs, we computed a formal FPP=1.48e-11which robustly qualifies HD 18599.01 as a statistically validated planet.
Additionally, we validated HD 18599.01 using the Python package TRICERATOPS20 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±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 TRICER-ATOPS 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.

DISCUSSION
Here, we consider the nature of HD 18599 b (=TOI 179 b) by placing it in context with the population of known exoplanets 21 .Figure 14 shows HD 18599 (marked as star) in the context of known transiting planets orbiting >1 Gyr old field stars (black points).Also shown are transiting planets orbiting young (<100 Myr, blue) and adolescent stars (100<age<1000 Myr, orange).HD 18599 b resides in the parameter space consistent with other planets orbiting adolescent stars.In particular, HD 18599 b is most similar to K2-284 b in terms of orbital period and radius among the known young transiting planets.K2-284 b also orbits a K-type field star with an age of 100-760 Myr, consistent with HD 18599 (David et al. 2018).When compared to transiting planets around older field stars, HD 18599 b appears to fall in the large-radius tail of the size distribution for close-in sub-Neptunes.
The measured period of  orb =4.13 d and radius of  P =2.7  ⊕ places it close to the boundaries of the Neptunian desert defined by Mazeh et al. (2016) as shown in the lower right panel of Figure 14.HD 18599 b joins the sparsely populated region within or near the boundaries of the Neptune desert, with the youngest among planet groups occupying majority of such special region.Among the known young planets, K2-100 b is located the deepest within the Neptune desert with  orb =1.67 d and  P =3.88  ⊕ .This is unsurprising because the planet's atmosphere is still currently experiencing photoevaporative escape (Gaidos et al. 2020).The youngest in the group is K2-33 b at 5-10 Myr with  orb =5.42 d and  P =5.04  ⊕ .
In the two lower panels in Figure 14, it is clear that the young planet population appears inflated relative to the older planet populations, especially in the case of planets orbiting low-mass stars (e.g.Mann et al. 2016;Rizzuto et al. 2020;David et al. 2019;Newton et al. 2019;Bouma et al. 2020).However, most transiting planets found orbiting adolescent stars like HD 18599 do not appear to be clear outliers in the period-radius diagram (e.g.Mann et al. 2020;David et al. 2018;Livingston et al. 2018a).One possible explanation for this observed behavior comes from photoevaporation theory.Young stars emit X-ray and EUV radiation strong enough to drive atmospheric escape.In this picture, atmospheric loss operates in a timescale on the order of hundreds of millions of years for 21 Based on a query of the NASA Exoplanet Archive "Confirmed Planets" table on 2022 January 31, https://exoplanetarchive.ipac.caltech.edu/close-in planets with thick atmospheres to a few billion years for the largest and most massive planetary cores (David et al. 2020).
Another explanation for the anomalously large radii of young planets comes from core-powered mass loss theory (Gupta & Schlichting 2019).In this picture, the cooling of planetary core provides the energy for atmospheric loss.This effect alone effectively reproduces the observed valley in the radius distribution of small close-in planets (Fulton et al. 2017).More recently, Berger et al. (2020) found the first evidence of a stellar age dependence of the planet populations straddling the radius valley.They found that the fraction of super-Earths to sub-Neptunes increases from 0.61 ± 0.09 at young ages (<1 Gyr) to 1.00 ± 0.10 at old ages (>1 Gyr), consistent with the prediction by core-powered mass loss which operates on a gigayear timescales.
For a mini-Neptune with a predominantly H2-He envelope presumably similar to HD 18599 b, its radius will contract to approximately half of its original size in the first 500 Myr due to radiative cooling and XUV-driven mass loss (Howe & Burrows 2015).Given the wide age range of HD 18599 b, it is difficult to ascertain whether enough time has passed to allow HD 18599 b to contract to its final equilibrium radius and whether either one or combination of both effects is responsible for the observed properties of HD 18599 b.At this point, we cannot be sure that the relatively large size of HD 18599 b is due to its young age.However, the properties of HD 18599 b do not appear to be merely a consequence of observational bias.

SUMMARY AND FUTURE PROSPECTS
Young exoplanets inhabit a very important part of the exoplanet evolutionary timescale, where formation mechanisms, accretion, migration and dynamical interactions can significantly change the shape of observed planetary systems.To date, only a handful of planetary systems <1 Gyr old are known.The main reason for the rarity of known young planets is the strong stellar activity of young stars, which makes it hard to find the subtle planetary signal in the face of large stellar variations (e.g.Plavchan et al. 2020;Barragán et al. 2019).Despite the challenges, the young planet population is an emerging field that is expected to yield highly impactful scientific results.Therefore, by compiling a statistically significant sample of well-characterized exoplanets with precisely measured ages, we should be able to begin identifying the dominant processes governing the time-evolution of exoplanet systems.In this light, we present the discovery and validation of a sub-Neptune orbiting the young star HD 18599.Complementary to TESS data, we utilized a suite of follow-up data including photometry with Spitzer/IRAC, IRSF/SIRIUS, & KELT, speckle imaging with Gemini/Zorro, and high resolution spectroscopy from ESO 2.2m/FEROS, LCO 2m/NRES, SMARTS 1.5m/CHIRON, & MKO 4x0.7m/.By implementing a similar validation framework and analyses presented previously, we found that the planet has an orbital period of 4.13 d, a radius of 2.7  ⊕ , a mass of <0.4 JUP , and an equilibrium temperature of 934±10 K.When compared to transiting planets around older field stars, HD 18599 b appears to fall in the large-radius tail of the size distribution for close-in sub-Neptunes.
A comparison between the typical densities of young and old planets may be more insightful than simply comparing radii and periods.However, planets around young stars are challenging for RV observations and only a handful have been successful to obtain mass measurements (e.g.Barragán et al. 2019Barragán et al. , 2022;;Zicher et al. 2022).Fortunately, the proximity and brightness of HD 18599 b would allow mass measurement with high resolution Doppler spectroscopy, unlike the majority of young K2 planets orbiting faint host stars.Measuring the planet's mass provides a rare opportunity to derive the planet's density and model its interior composition given sufficient precision.This would then allow direct comparison of bulk parameters such as density with similar sub-Neptunes orbiting older host stars.Our estimate of mass upper limit from Section 2.2 implies an RV amplitude of at most 300 m s −1 .Thus, a dedicated campaign designed to observe and model the RV systematics may allow a precise mass measurement of HD 18599 b.A potential caveat is that HD 18599 is rotating faster than 15 days which results in rotational broadening at a level that limits RV precision and/or a large stellar RV jitter caused by strong stellar variability expected for young stars.
Moreover, a measurement of the planet's obliquity or the skyprojected angle between the stellar spin axis and the planet's orbital axis via observation of the Rossiter-McLaughlin (RM) effect is also a good avenue to explore.The planet's obliquity is a tracer for any dynamical processes, such as tidal interaction with stellar companions (e.g.Batygin 2012), that the system might have undergone to produce a misaligned orbit, assuming well-aligned initial condition.
Obliquity measurements were conducted for a growing number of young planets (e.g.Johnson et al. 2022;Wirth et al. 2021;Gaidos et al. 2020).These planets were found to have well-aligned orbits which indicates that misalignments may be generated over timescales of longer than tens of Myr.We estimate the RV amplitude due to the RM effect to be ≈140 m s −1 using the relation ΔRV RM =  sin  (  /  ) 2 √ 1 −  2 (Gaudi & Winn 2007), making this one of the most amenable targets for RM follow-up.This makes obliquity measurements feasible for HD 18599 b, following only DS Tau Ab and AU Mic b to have obliquity measurement among young systems so far (Zhou et al. 2020;Martioli et al. 2020).Constraining the star-planet obliquity of this young system may shed light on possible migration mechanisms (e.g.Kozai, secular, tidal interactions) to explain its current architecture such as its seemingly short period as compared to other young planets (e.g.Montet et al. 2020).
Following Kempton et al. (2018), we compute the transmission spectroscopy metric (TSM) which is a general metric useful for prioritizing of transmission spectroscopy targets for future infrared observations e.g. with JWST (McElwain et al. 2020).It is defined as where  is a scale factor equal to 1.28 appropriate for HD 18599 b,   is the planet radius in Earth radii,   is the planet mass equivalent to 1.436 1.7  appropriate for the size of HD 18599 b,  * is the star radius in solar radii,   is the apparent magnitude of the host star in the J band, and  eq is the planet's equilibrium temperature in K, where  is the orbital semimajor axis in solar radii, and assuming zero albedo and full day-night heat redistribution.We computed TSM=21 which can be interpreted as a transmission spectrum measurement with signal to noise of 21 in a 10 hour window with the NIRISS instrument (Louie et al. 2018).Combined with its young age, this makes HD 18599 a compelling target for atmospheric characterization.

Figure 1 .
Figure 1.3" × 3" DSS2 (blue filter) image showing the target (green square) on the center and nearby Gaia sources (orange circles), superposed with the Sector 2 SPOC photometric aperture (blue polygon).We establish a dilution of <1% for HD 18599 based on the flux contributions of all sources within or near the aperture perimeter.

Figure 2 .
Figure 2. Raw light curves used in transit modeling superposed with best-fit transit model with baseline trend.The top row shows the TESS PDCSAP light curves (cadence=2 min) in sectors 2 & 3 where the individual transits are highlighted in blue and numbered relative to the first (n=0).Panels 2 to 12 show the zoomed-in view of individual transits with best-fit model.The last panel shows the Spitzer light curve.

Figure 3 .
Figure 3. Same as Figure 2 but for TESS sectors 29 and 30.Note the significantly higher sampling rate due to the shorter (20-second) cadence.

Figure 4 .
Figure 4. Light curves used in this work.The top row shows the phase-folded data with best-fit transit model (See Section 4.1).The blue data points show binned data every 2-minutes and red lines show a fit using random samples from the MCMC posterior.The bottom row shows the residuals between the data and the transit+systematics model.

Figure 5 .
Figure 5.The contrast curves in 562 and 832nm taken by Gemini/Zorro speckle taken by the Gemini-South telescope.The inset is a reconstructed speckle image.

Figure 6 .
Figure 6.Spectral energy distribution (SED) of HD 18599.Red symbols represent the observed photometric measurements, where the horizontal bars represent the effective width of the passband.Blue symbols are the model fluxes from the best-fit NextGen atmosphere model (black).

Figure 7 .
Figure 7. FEROS spectra zoomed around the Li I doublet absorption at 6708 Å (shaded blue region) which is a known strong indicator of youth.The black line is the mean of the all spectra taken at different epochs.

Figure 9 .
Figure 9. Rotation period estimate.The top panel shows the Lomb-Scargle periodogram of the TESS light curves (red) and KELT observations (black).The rotation period is consistently detected in both TESS and the six years of monitoring from KELT.The bottom panel shows the TESS light curves folded to the rotation period.Each rotation period is over-plotted with a slightly different color gradient, where the lighter/darker color represents more recent observations.

Figure 10 .
Figure 10. eff vs. rotation period relation of HD 18599 (white star) compared to representative star cluster members with known ages.The black dashed line shows the negative empirical relation between rotation period and  eff .

Figure 11 .
Figure 11.Summary of age estimates for HD 18599 using various methods.

Figure 12 .
Figure 12.Comparison of marginalized posterior distributions of fits using

Figure 13 .
Figure 13.Phase-folded RVs with Keplerian models corresponding to the 1-, 2-, and 3- mass upper limits fitted to the M -Australis (green) and FEROS (orange) data.Error bars with dark shades show the errors estimated from our spectroscopic analyses.The error bars with lighter shades show the original errors + jitter term value (added in quadrature) from the best-fit RV model (black line, best M p = 7.4 ⊕ ).

Figure 14 .
Figure 14.HD 18599 (marked as star) in the context of known transiting planets orbiting older field stars (gray points).Also shown in orange and blue circles are transiting planets with age<100 Myr and 100<age<1000 Myr, respectively.The contours represent the number density of the known transiting planets orbiting older host stars.HD 18599 b appears to close to the boundaries of the Neptunian desert (solid black lines, Mazeh et al. (2016)) in the period-radius plane in the lower right panel.The dashed lines refer to the boundaries' uncertainty regions.

Table 1 .
Summary of our follow up photometric observations.
NOTE-This table is published in its entirety in a machine-readable format.A few entries are shown for guidance regarding form and content.
† 0.1 mas was added in quadrature.

Table 5 .
Results of joint transit modeling of TESS and Spitzer light curves.