A dense mini-Neptune orbiting the bright young star HD 18599

ABSTRACT

Very little is known about the young planet population because the detection of small planets orbiting young stars is obscured by the effects of stellar activity and fast rotation, which mask planets within radial velocity and transit data sets. The few planets that have been discovered in young clusters generally orbit stars too faint for any detailed follow-up analysis. Here, we present the characterization of a new mini-Neptune planet orbiting the bright (V = 9) and nearby K2 dwarf star, HD 18599. The planet candidate was originally detected in TESS light curves from sectors 2, 3, 29, and 30, with an orbital period of 4.138 d. We then used HARPS and FEROS radial velocities, to find the companion mass to be 25.5|$\pm$|4.6 |$M_{\oplus }$|⁠. When we combine this with the measured radius from TESS of 2.70|$\pm$|0.05 |$R_{\oplus }$|⁠, we find a high planetary density of 7.1|$\pm$|1.4 g cm|$^{-3}$|⁠. The planet exists on the edge of the Neptune Desert and is the first young planet (300 Myr) of its type to inhabit this region. Structure models argue for a bulk composition to consist of 23 per cent H|$_2$|O and 77 per cent Rock and Iron. Future follow-up with large ground- and space-based telescopes can enable us to begin to understand in detail the characteristics of young Neptunes in the galaxy.

techniques: photometric, techniques: radial velocities, stars: activity, stars: individual: HD 18599– Planetary Systems

1 INTRODUCTION

The Neptune Desert, also known as the sub-Jovian Desert, is a region in the period-mass-radius planetary parameter space relatively devoid of planets (Szabó & Kiss 2011). The trapezoidal shape of the desert was outlined by Mazeh, Holczer & Faigler (2016), with two firm edges that were easily identifiable in results drawn from Kepler Space Telescope (Borucki et al. 2010). It seems that planets closest to their host stars come in two flavours, either large gas giants, or small rocky planets. The desert can therefore be explained by the combined effects of tidal migration and planetary photo-evaporation by the host star (Owen et al. 2018). Giant planets are too massive to be significantly affected by the high-energy XUV radiation from the host star, no matter how close to their host star they migrate, and therefore they maintain a dominant gaseous envelope. Planets like Neptune, on the other hand, can not hold on to their gaseous envelopes, losing them early in their evolution, and ending up as left-over hot cores, with minimal residual atmospheres. This mechanism naturally explains the Neptune Desert region.

Although the Neptune Desert is ’dry’, it is not completely devoid of planets. Recently a few planets have been discovered residing close to their stars, yet with masses and radii in the Neptune regime. One of the first planets announced from the Next Generation Planet Search (NGTS; Wheatley et al. 2018), NGTS-4b (West et al. 2019), was a dense planet found to have a period, mass, and radius that placed it in the Neptune Desert. Subsequently, Transiting Exoplanet Survey Satellite (TESS) detected the planet TOI-849b (Armstrong et al. 2020), another dense super-Neptune within the desert. This planet likely hosts no atmosphere at all, probably the result of a previous gas giant that was entirely stripped of its gaseous envelope. A further TESS discovery was planet orbiting LTT 9779 a unique desert Neptune, having an orbital period of 19 h, it maintains an equilibrium temperature above 2000 K, making it the first Ultrahot Neptune (Jenkins et al. 2020). More intriguing is the fact that this planet does host a significant atmosphere (see Crossfield et al. 2020; Dragomir et al. 2020), even though the system is not young. Given that the Neptune desert is well-defined empirically and explained by models, the study of objects within the desert is particularly intriguing though bright examples to test photo-evaporation models are needed.

Since photo-evaporation appears to be the mechanism that gives rise to the desert, we expect the majority of mass-loss to happen very early in the evolution of the system. Stars emit the majority of their high-energy X-ray and UV radiation in the first 100 Myrs of their life (Jackson, Davis & Wheatley 2012), the photons that heat the upper atmospheres of planets sufficiently to cause gas to escape the planet’s gravitational potential. Other models that can explain planets to lose significant amounts of their atmospheres, generally happen on longer time-scales, like core-driven mass-loss (Gupta & Schlichting 2020) or Roche Lobe Overflow (Valsecchi, Rasio & Steffen 2014). Therefore, studying this region as a function of age will not only shed light on the dominant processes, but also allow a more detailed characterization of the timescales these processes operate over.

1.1 Significance of HD 18599

The southern star HD 18599 (TOI-179) has a V magnitude of 8.99 mag and had an estimated spectral type of K2V based on data from the Hipparcos mission (Perryman et al. 1997; van Leeuwen 2007). The Gaia EDR3 parallax gives rise to a distance of 38.632 pc (Gaia Collaboration 2021), and given the brightness, this is consistent with the aboslute magnitude expected for a K-dwarf. The star was screened as part of the first phases of the Calan-Hertfordshire Extrasolar Planet Search project (CHEPS; Jenkins et al. 2009, 2017) through measurements of the stellar chromospheric activity, metallicity, rotational velocity, and kinematics. They found the star to have a log|$R^{\prime }_{\rm {HK}}$| index of −4.39 dex (Jenkins et al. 2008, 2011), vsin(i) of 4.3 km s|$^{-1}$|⁠, and kinematic U, V, and W velocities of −3.49, 8.50, and −6.21 km s|$^{-1}$|⁠, respectively (Murgas et al. 2013). Applying gyrochronological-based relations to these measured values (e.g. see Mamajek & Hillenbrand 2008) provides estimates of the age of the system, coming in at around 300 Myrs.

Although the star was recognized to be too active for inclusion in the CHEPS precision radial velocity (RV) program, (since that study focused only on the most chromospherically quiet stars), the TESS (Ricker et al. 2014) mission does not have the same activity constraints, and a transiting planet candidate was detected as part of that mission. The candidate was subsequently validated as a planet by (de Leon et al. 2022, in preparation), providing another possible planet on the edge of the Neptune Desert. However, unlike most of the other systems in this region, we can confidently say this is younger than typical field stars, and given the brightness of the host, it could allow more stringent constraints to be placed on the time-scales of the photo-evaporative process that is believed to drive the existence of the desert, along with allowing atmospheric characterization of the younger cohort of planets.

2 OBSERVATIONS

This candidate was first detected through the observation of 16 transit events by the TESS mission, yielding a period of 4.137 d from an EXOFASTv2 (Eastman et al. 2019) analysis listed on the Exoplanet Follow-up Observing Program.¹ After scrutinizing the processed light curve publicly available at MAST,² and given that the star is bright, we selected this candidate as a high-priority target to be observed in a HARPS campaign starting in 2018. Upon examination of archival data, we found the star was already observed with sparse sampling between 2014 and 2017. Later, additional observations with FEROS were planned to help constrain the orbital and physical parameters of the planet. Finally, Minerva-Australis observations to help validate the candidate were also acquired. The following sections describe the observations and data acquisition of HD 18599 in detail.

2.1 Photometry

2.1.1 WASP photometry

WASP-South, an array of 8 wide-field cameras, was the Southern station of the WASP transit-search project (Pollacco et al. 2006). It observed the field of HD 18599 in 2010 and 2011, when equipped with 200-mm, f/1.8 lenses, and then in 2012, 2013, and 2014, equipped with 85-mm, f/1.2 lenses. Observations were made on each clear night over an observing season of 175 nights per year, accumulating 86 000 photometric data points with a typical 15-min cadence. At V = 9, HD 18599 is by far the brightest star in the 48-arcsec extraction aperture of the 200-mm data. A V = 10.9 star, 88 arcsec away, is, however, within the 113-arcsec extraction aperture of the 85-mm data, and will cause a |$\sim$|10 per cent dilution. We studied the periodogram of the WASP light curves in Section 3.2 to find the rotation period of HD 18599 in order to properly detrend the TESS data.

2.1.2 TESS photometry

HD 18599 was observed by TESS in sectors 2, 3, 29, and 30 with Camera 3 in short-cadence mode (2 min), starting from 2018 August 22 to the 2020 September 22, producing a total time baseline of 2 yr, during which 21 transit features were observed. Observations were analysed and processed with the Science Processing Operations Center (SPOC; Jenkins et al. 2016) by NASA Ames Research Center, after which it was promoted to TOI-179.

We extracted the Simple Aperture Photometry (SAP) and the Presearch Data Conditioning (PDC; Jenkins et al. 2016) light curves that were produced by the SPOC pipeline, removed data points flagged as low quality, and finally we normalized the light curve. We note that even after this procedure, stellar variability still remains in the light curve and must be taken care of to determine the most optimal companion parameters. In the end, the SAP and PDCSAP light curves look almost identical, but we use the latter throughout the paper in any case.

Finally we applied Gaussian Processes (GP) to detrend and remove the stellar variability signal in the light curve (see Section 3.3 ). The detrended light curve can be found in Table 1.

Table 1.

Detrended TESS photometry for HD 18599. The full table is available in a machine-readable format from the online journal. A portion is shown here for guidance.

Time	Flux	Flux
(BJD-2450000)	(normalized)	error
8354.11157	0.9989	0.0004
8354.11296	0.9992	0.0004
8354.11574	0.9993	0.0004
8354.11712	0.9994	0.0004
8354.11851	0.9994	0.0004
8354.11990	0.9999	0.0004
8354.12129	0.9990	0.0004
8354.12268	0.9993	0.0004
8354.12407	0.9993	0.0004
8354.12546	0.9995	0.0004
|$\vdots$|	|$\vdots$|	|$\vdots$|

2.2 Radial velocities

We observed HD 18599 with three fiber-fed high precision échelle spectrographs: High-Accuray Radial velocity Planet Searcher (HARPS; Pepe et al. 2002), FEROS, and Minerva-Australis to obtain high-precision RVs and constrain the planetary mass.

2.2.1 HARPS

HARPS is mounted on the ESO 3.6-m Telescope at the La Silla observatory in Northern Chile, and has a spectral resolution of R = 120 000. Observations include data acquired prior to the fiber updates we refer to as ‘pre’ (PI: Lagrange). Data acquired after the fiber updates in 2015 (PI: Lagrange, Díaz, Berdiñas, Jordán, Brahm, and Benatti)³ are labelled as ‘post’. In total, we obtained 103 observations from HARPS. HARPS observations are processed on site using the standard ESO Data Reduction Software (DRS). The DRS pipeline data products were then reprocessed independently using the HARPS-TERRA code (Anglada-Escudé & Butler2012). HARPS-TERRA creates a high signal-to-noise template from the individual observations. Then, each observation is matched using a |$\chi ^2$| process relative to this template producing an RV for each observed spectrum. The code also computes stellar activity indices for the Calcium ii H & K lines (⁠|$\lambda _{H}$| = 3933.664 Å, |$\lambda _{K}$| = 3968.470 Å), which are not available from the DRS-reduced spectra. For the derivation of S-indices, the code integrates the flux in these lines and compares with the flux on adjacent chunks in the continuum, following the procedure described in Lovis et al. (2011). S-indices from HARPS-TERRA are calibrated to the Mt. Wilson system (⁠|$S_{MW}$|⁠) and they can be used as a direct proxy to monitor the chromospheric activity of the star. Additional activity indices such as the full-width at half-maximum of the cross-correlation function (CCF FWHM) and the bisector inverse slope (BIS) are taken directly from the fits headers of the DRS-reduced data products. Uncertainties in BIS are taken as twice the RV errors and the CCF FWHM uncertainties are 2.35 times the RV errors (Zechmeister et al. 2018).

2.2.2 FEROS

The FEROS spectrograph is mounted on the 2.2-m MPG/ESO telescope also at La Silla observatory and has a spectral resolution of R = 48 000 (Kaufer et al. 1999). Nine observations were made during September 10 and 19 under program 0103.A-9004(A) (PI: Vines) which were reduced with the CERES pipeline (Brahm, Jordán & Espinoza 2017). After the reduction, CERES calculates RVs through the CCF method. We used a binary K5 mask for each epoch and fit a double Gaussian to the CCF in order to find the RV. A double Gaussian was used in order to account for scattered moonlight contamination. In addition to the RVs, CERES also provides the CCF FWHM and the BIS as activity indices. The uncertainties in BIS are calculated internally by CERES while the uncertainties of the CCF FWHM are calculated by dividing the standard deviation of the Gaussian fitted to the CCF, and dividing by the S/N at 5130 Å, a procedure analogous to the RV uncertainty calculation done by CERES. S-indices index are computed from each 1D spectra after correcting to rest-frame, following the procedure outlined by Jenkins et al. (2008) using the Ceres-plusplus⁴ code.

2.2.3 Minerva-Australis

Minerva-Australis is an array of four PlaneWave CDK700 telescopes located in Queensland, Australia, fully dedicated to the precise RV follow-up of TESS candidates (e.g. Jordán et al. 2020; Addison et al. 2020a, b). The four telescopes can be simultaneously fiber-fed to a single KiwiSpec R4-100 high-resolution (R = 80 000) spectrograph (Barnes et al. 2012; Addison et al. 2019).

HD 18599 was monitored by Minerva-Australis in its early commissioning period, with a single telescope between 2019 January 6 and 29. Exposure times were 30 min, with some nights receiving two or three consecutive exposures. We obtained a total of 31 observations, with a mean S/N of 23|$\pm$|5 per pixel. RVs 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.

We show the RVs with the activity indicators in Table 2.

Table 2.

RVs and activity indices for HD 18599.

BJD	RV	RV err	BIS	BIS err	S|$_{\text{M}W}$|	S|$_{\text{M}W}$| err	FWHM	FWHM err	instrument
(-2450000)	(m s|$^{-1}$|⁠)	(m s|$^{-1}$|⁠)	(m s|$^{-1}$|⁠)	(m s|$^{-1}$|⁠)	(dex)	(dex)	(km s|$^{-1}$|⁠)	(km s|$^{-1}$|⁠)
6980.6959	12.7646	1.2472	5.2533	0.8099	|$-$|0.0053	0.0028	|$-$|0.0356	0.0029	HARPS_pre
6980.7020	11.4832	1.2733	5.8890	0.8002	|$-$|0.0034	0.0028	|$-$|0.0337	0.0010	HARPS_pre
6982.6788	16.0365	1.3771	8.2058	1.0231	0.0031	0.0033	0.0113	0.0032	HARPS_pre
6982.6849	14.4035	1.9666	10.6088	1.0709	0.0214	0.0034	0.0067	0.0046	HARPS_pre
6987.6309	|$-$|18.5854	1.9749	55.0689	1.1134	|$-$|0.0055	0.0035	0.0504	0.0046	HARPS_pre
6987.6371	|$-$|16.9837	1.8491	52.7322	1.1613	0.0077	0.0036	0.0564	0.0043	HARPS_pre
6988.6928	|$-$|19.1186	1.2100	44.4527	0.8029	|$-$|0.0180	0.0026	|$-$|0.0556	0.0028	HARPS_pre
|$\vdots$|	|$\vdots$|	|$\vdots$|	|$\vdots$|	|$\vdots$|	|$\vdots$|	|$\vdots$|	|$\vdots$|	|$\vdots$|	|$\vdots$|

3 ANALYSIS

3.1 Stellar properties

3.1.1 SPECIES

We analysed the 1D stacked HARPS spectra using SPECIES (Soto & Jenkins 2018; Soto, Jones & Jenkins 2021). SPECIES is an automated code to derive stellar parameters for large samples of stars, using high-resolution echelle spectra. It makes use of equivalent widths from a number of neutral and ionized iron lines (measured with EWComputation, Soto et al. 2021) to derive the atmospheric parameters (temperature, metallicity, surface gravity, and microturbulence). Together with ATLAS9 model atmospheres (Castelli & Kurucz 2004), it solves the radiative transfer and hydrostatic equilibrium equations using MOOG (Sneden 1973), imposing local thermodynamic equilibrium (LTE) conditions, as well as excitation and ionization equilibrium. The correct atmospheric parameters are found when the abundance for ionized and neutral iron are the same, the obtained iron abundance is the same as the one used to create the model atmosphere, and there is no correlation between the neutral iron abundances with the excitation potential, and with the reduced equivalent widths (⁠|$W/\lambda$|⁠). Finally, rotational and macro turbulent velocity were derived using spectral line fitting and analytic relations, respectively. SPECIES finds an effective temperature of 5109 |$\pm$| 50 K, a |$\log ~g$| of 4.41|$\pm$|0.07 dex, a [Fe/H] of –0.01 |$\pm$| 0.05 dex and a vsini of 5.2 |$\pm$| 0.2 km s|$^{-1}$|⁠.

3.1.2 ARIADNE

We used the temperature, |$\log ~g$| and [Fe/H], found by SPECIES as Gaussian priors on a Spectral Energy Distribution (SED) analysis performed with ariadne⁵ (Vines & Jenkins 2022) to find the final set of bulk stellar parameters. ariadne is a publicly available python tool designed to fit catalog photometry to different stellar atmosphere models grids, Phoenix V2 (Husser et al. 2013); BT-Settl, BT-Cond, BT-NextGen (Hauschildt, Allard & Baron 1999; Allard, Homeier & Freytag 2012), Castelli & Kurucz (2004) and Kurucz (1993); which have been convolved with the following filter response functions:

• Johnson UBV

• Tycho-2 BtVt

• 2MASS JHK|$_\text{s}$|

• SDSS ugriz

• Gaia DR2v2 G, RP, and BP

• GALEX NUV and FUV

• TESS

The SEDs are modelled by interpolating the model grids in T|$_{\rm eff}$|-|$\log ~g$|-[Fe/H] space, scaling the synthetic flux by |$(R/D)^2$| and accounting for interstellar extinction through the extinction in the V band, A|$_{\rm V}$|⁠. Additionally, we include an excess noise parameter for every photometric observation in order to account for underestimated uncertainties. ARIADNE samples the posterior space using the dynesty (Speagle 2020) implementation of (dynamic) nested sampling (Higson et al. 2019). Nested sampling (Skilling 2004, 2006) algorithms are designed to estimate the Bayesian evidence of a model producing posterior distributions as a by-product.

We modelled the SED with four models, Phoenix V2, BT-Settl, Castelli & Kurucz (2004), and Kurucz (1993) and finally averaged the posterior parameters from each model, weighting them by their respective relative probability that was computed from their Bayesian evidence estimates, resulting in precise parameters that take into account the different micro-physics and geometry of the different models.

Additional priors for distance, radius, and A|$_{\rm v}$| were drawn from the Bailer–Jones distance estimate from Gaia DR2 (Bailer-Jones et al. 2018), a uniform prior from 0.5 to 20 R|$_\odot$| and a uniform prior from 0 to 0.038, the maximum line-of-sight extinction from the re-calibrated SFD galactic dust map (Schlegel, Finkbeiner & Marc 1998; Schlafly & Finkbeiner 2011).

In Fig. 1, we show the SED for HD 18599, and in Table 3 we report the relevant observational properties with our derived properties, along with the method used.

Figure 1.

The best-fitting Castelli & Kurucz SED (black line) for HD 18599 based on the photometric data (cyan points) presented in Table 3 is shown in the top panel. Horizontal errorbars show the bandpass width. The pink diamonds show the synthetic magnitudes at the wavelengths of the photometric data. The lower panel shows the residuals to the best-fitting model.

Table 3.

Stellar properties for HD 18599.

Property	Value	Source
Astrometric properties
RA	|$02^{\rm {h}} 57^{\rm {m}} 02{_{.}^{\rm s}}8835$|	Gaia
Dec.	|$-56^{\circ }11^\prime 30{_{.}^{\prime\prime}} 7297$|	Gaia
2MASS I.D.	J02570294-5611314	2MASS
TIC	207141131	TESS
Gaia DR2 I.D.	4728513943538448512	Gaia
|$\mu _{{\rm R.A.}}$| (mas y|$^{-1}$|⁠)	|$-$|36.676 |$\pm$| 0.042	Gaia
|$\mu _{{\rm Dec.}}$| (mas y|$^{-1}$|⁠)	50.599 |$\pm$| 0.045	Gaia
|$\varpi$| (mas)	25.9011 |$\pm$| 0.0244	Gaia
Photometric properties
V (mag)	8.968 |$\pm$| 0.012	APASS
B (mag)	10.015 |$\pm$| 0.004	APASS
g (mag)	10.20 |$\pm$| 0.05	APASS
V|$_\text{T}$| (mag)	9.084 |$\pm$| 0.016	Tycho2
B|$_\text{T}$| (mag)	10.015 |$\pm$| 0.004	Tycho2
G (mag)	8.7312 |$\pm$| 0.0007	Gaia
RP (mag)	8.1384 |$\pm$| 0.0019	Gaia
BP (mag)	9.2113 |$\pm$| 0.0020	Gaia
TESS (mag)	8.1796 |$\pm$| 0.0060	TESS
J (mag)	7.428 |$\pm$| 0.018	2MASS
H (mag)	7.029 |$\pm$| 0.015	2MASS
Ks (mag)	6.883 |$\pm$| 0.020	2MASS
NUV (mag)	15.684 |$\pm$| 0.012	GALEX
FUV (mag)	21.389 |$\pm$| 0.264	GALEX
Derived properties
T|$_{\rm eff}$| (K)	5083 |$\pm$| 23	ariadne
|$\left[{\rm Fe}/{\rm H}\right]$|	|$-$|0.05 |$\pm$| 0.04	ariadneARIADNE
vsini (km s|$^{-1}$|⁠)	5.2 |$\pm$| 0.2	SPECIES
log g	4.40 |$\pm$| 0.07	ariadne
|${\rm M}_{\star }$|(⁠|$M_{\odot }$|⁠)	0.807|$^{+0.019}_{-0.007}$|	ariadne
|${\rm R}_{\star }$|(⁠|$R_{\odot }$|⁠)	0.798|$^{+0.006}_{-0.007}$|	ariadne
|$\rho$| (g cm|$^{-3}$|⁠)	2.241|$^{+0.081}_{-0.077}$|	This work
Age (Myr)	300	Gyrochronology
Distance (pc)	38.585|$^{+0.110}_{-0.150}$|	ariadne
Av (mag)	0.030|$^{+0.002}_{-0.014}$|	ariadne
|$P_{\rm rot}$|	8.74 |$\pm$| 0.05	This work

2MASS (Skrutskie et al. 2006); UCAC4 (Zacharias et al. 2013); APASS (Henden & Munari 2014); WISE (Wright et al. 2010); Gaia (Gaia Collaboration 2016, 2018); GALEX (Bianchi et al. 2011); TESS (Stassun et al. 2018).

|$q_{1,2}$| are the Kipping Limb Darkening parameters.

HD 18599 is previously known to be a young field star based on previous studies (Grandjean et al. 2020, 2021). In a more recent study, de Leon et al. 2022 employed various methods to estimate the age of HD 18599. In that paper, we measured an equivalent width (EW) of 42 |$\pm$| 2 mÅ for Li I |$\lambda$|6708. This measurement generally agrees with those found in stars in Hyades and Praesepe corresponding to an age of |$\sim$|800 Myr. We also estimated the stellar age by taking advantage of the observed chromospheric activity together with empirical age–activity–rotation relations. In particular, we used the chromospheric activity indicator, |$\log R^{\prime }_{HK} = -4.41 \pm 0.02$| from (Boro Saikia et al. 2018), which predicts an 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.

HD 18599 exhibits significant spot-modulated rotational signals in the four sectors of TESS observations. Using a periodogram analysis, we measured the rotation period of HD 18599 for each of these sectors of observations, finding that it has a mean rotation period of |$8.71\pm 0.31$| d. We compute a gyrochronological age estimates of 386 Myr, with |$3-\sigma$| range of 261–589 Myr based on Mamajek & Hillenbrand (2008) model, whereas we compute an age of 247 Myr, with 3-|$\sigma$| range of 185–329 Myr based on Barnes (2007) model. 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 rotation period of 8.69 d reported by KELT. We adopt a median age of 300 Myr throughout the rest of the paper.

An asteroseismology study of HD 18599 was also performed, finding that the active nature of the star produces photometric variability much larger than the amplitude of the expected oscillation modes, making the visual detection of the oscillations very tough. The details of the analysis and the results can be found in Samadi-Ghadim et al. (2022 in preparation).

3.2 Stellar rotation

From Kepler data, it has been shown that the analysis of photometric time series can yield good estimates of the stellar rotation period (see, e.g. McQuillan, Mazeh & Aigrain 2013; Giles, Collier Cameron & Haywood 2017). Rotationally modulated spots on the stellar surface can produce a periodic signal that can be detected in the photometry, since these spots are regions of diminished flux. Additionally, rotational modulations can be a source of false-positive exoplanets signatures, for which the TESS photometry, combined with ground-based data sources, can serve as a powerful diagnostic for assessing the properties of known exoplanets (Kane et al. 2021; Simpson et al. 2022).

We searched each season’s WASP light curve for rotational modulations using the methods presented in Maxted et al. (2011). There is a clear and persistent modulation seen in four of the five seasons (Fig. 2). The mean period is 8.74 |$\pm$| 0.06 d, the amplitude ranges from 6 to 10 mmag, and the false-alarm likelihoods are less than 10|$^{-3}$|⁠. In 2013, the modulation is marginally present, but with a much lower amplitude of 1–2 mmag.

Figure 2.

Periodograms of the WASP-South data for HD 18599 (left) along with the modulation profiles from folding the data (right). The red line is at a period of 8.74 d.

We also computed the autocorrelation function (ACF; McQuillan et al. 2013) for each TESS light curve, using the PCDSAP corrected data provided by the SPOC in order to find a possible value for the stellar rotation period. We follow the procedure described in Soto et al. (2018), where we model the ACF as an underdamped Simple Harmonic Oscillator (SHO; López-Morales et al. 2016). We define

$$\begin{eqnarray*} y(t) = e^{-t/\tau _{\rm AR}} \left[ A \cos \left(\frac{2\pi t}{P_{\rm rot}}\right) + B \cos \left(\frac{4\pi t}{P_{\rm rot}} \right) \right] + y_{0} \end{eqnarray*} $$

(1)

Where |$\tau _{\rm AR}$| is the correlation time-scale and |$P_{\rm rot}$| is the rotation period.

We estimated |$P_{\rm rot}$| using the Edelson & Krolik method (Edelson & Krolik 1988) implemented in astroML python package via a least-squares-fitting process (Vanderplas et al. 2012). Results are tabulated in Table 4.

Table 4.

ACF results from the TESS light curve.

Sector	|$P_{\rm rot}$|	|$P_{\rm rot}$| err
	(d)	(d)
2	4.522	0.001
3	4.526	0.002
29	8.605	0.281
30	8.791	0.110

The rotation period from the first two sectors are consistent with each other and consistent with being an alias from the rotation period derived from WASP data. This is further confirmed with the estimated rotation period of the latter sectors being consistent with the one derived by WASP as well.

We ran the Generalized Lomb–Scargle periodogram (GLS; Zechmeister & Kürster 2009) in a grid of 50 000 trial frequencies ranging from 1 d to 2 yr in period space on the TESS data to find the possible sinusoidal signals. We show these in Fig. 3 where the dashed line represent the 0.1 per cent, 5 per cent, and 10 per cent false-alarm probabilities, but they happen to coincide due to the amplitude of the oscillations. The top panel of Fig. 3 shows the GLS periodogram of sectors 2 and 3, where the two most strong signals correspond to 4.41 and 8.78 d. The middle panel shows the periodogram of sectors 29 and 30 where the most dominant signal is 8.76 d, the second most strong signal is 22.89 d. We show a third signal, which while not as significant, corresponds to the same period as the most strong signal found in sectors 2 and 3, 4.40 d. The bottom panel shows the periodogram of the combined three sectors where the most dominant signal is the 8.77 d one and the 4.40 d signal is the second most significant one.

Figure 3.

From top to bottom: the GLS periodogram for the TESS light curves from sectors 2 to 3 (top), sector 29 to 30 (middle), and the three sectors combined (bottom). The orange crosses show the most significant peaks with their respective periods labelled on top. The vertical red line marks the 4.1375-d planet candidate from the TESS data. The dashed horizontal line are the false-alarm probabilities for 0.1 per cent, 5 per cent, and 10 per cent, they turn out to be superimposed.

Using juliet⁶ (Espinoza, Kossakowski & Brahm 2019) and multinest via the python wrapper pymultinest (Feroz et al. 2019; Buchner et al. 2014), we fit a GP using a SHO kernel from celerite (Foreman-Mackey et al. 2017) to the out-of-transit TESS light curve in order to remove the stellar activity signal present. We used the rotation period derived from WASP as a Gaussian prior for the characteristic frequency of the SHO. We found the best-fitting parameters to be |$S_0 = 0.0008\pm 0.0001$|⁠, |$Q = 0.11\pm 0.01$|⁠, and |$\omega _0 = 0.719\pm 0.004$|⁠. A quality factor less than |$1/2$| means the system is overdamped and thus there is no oscillation aside from the rotation of the star, which is captured by the GP as seen in the top panel of Fig. 4. These parameters translate to a rotation period of |$8.74\pm 0.05$| d. The top panel of Fig. 4 shows the raw light curve with the best-fitting GP, and the bottom panel shows the detrended light curve, with the red vertical lines highlighting the transit features.

Figure 4.

The normalized TESS light curve including data from the three sectors, with the best-fitting SHO kernel GP (black line). The bottom panel shows the detrended light curve in ppm. The vertical red lines show the transit features.

3.3 Activity indices

We searched for correlations between the RVs and the activity indices: (BIS), CCF FWHM and S|$_{\rm MW}$|⁠, and we found a strong anticorrelation between the BIS and RVs, no correlation for the CCF FWHM, and a weak correlation for S|$_{\rm MW}$|⁠, quantified by the Pearson r correlation coefficients: −0.62, 0.05, and 0.49, respectively. We show our thresholds for correlation in Table 5.

Table 5.

Adopted thresholds for the Pearson r correlation coefficient. We only show the positive values, but the same thresholds apply for their negative counterparts.

r	Strength
0.00–0.19	Very weak
0.20–0.39	Weak
0.40–0.59	Moderate
0.60–0.79	Strong
0.80–1.00	Very strong

In addition to the Pearson r coefficients, we fit a linear model to each of the correlations using Markov Chain Monte Carlo (MCMC) through emcee (Foreman-Mackey et al. 2013) to probe the posterior parameter space. We employed a Gaussian likelihood that takes into account uncertainties in both axes. We show the best-fitting model parameters in Table 6 and the correlations in Fig. 5. Given that the slopes of the CCF FWHM and S|$_{\rm MW}$| are statistically consistent with zero, we conclude that the only significant correlation is the BIS.

Figure 5.

From left to right: the bisector velocity span with 1α, 2α, and 3α confidence regions, FWHM, and S-index correlation with the RVs after subtracting their respective mean values. HARPS_pre, HARPS_post, and FEROS are denoted by circles, hexagons, and triangles, respectively. Each observation is colour coded by its BJD timestamp. In all cases, the blue line corresponds to the best-fitting linear model. Confidence regions were not added for the FWHM and S index panels due to them being unconstrained.

Table 6.

Intercept and slopes of the best-fitting linear models of each correlation.

Correlation	Slope	Intercept
BIS	|$-0.836^{+0.088}_{-0.096}$|	17.302|$^{+2.222}_{-2.254}$|
CCF FWHM	0.000|$^{+0.005}_{-0.005}$|	|$-0.032^{+0.161}_{-0.162}$|
S|$_{\rm MW}$|	0.001|$^{+0.005}_{-0.005}$|	|$-0.001^{+0.162}_{-0.162}$|

Using the GLS periodogram we searched for periodicities in the RVs, BIS, S|$_{\rm MW}$| indices, and CCF FWHM. Signals at 1, 0.5, and 0.3 d arise from the window function, as well as a broad signal around 800–1000 d, which reappears in the RVs, S|$_{\rm MW}$|⁠, and the FWHM (see Fig. 6). The BIS shows a peak at 4.40 d, which coincides with the Prot alias from TESS and WASP. Given that the period candidate is very close to the rotation period alias tracked by the BIS, a straightforward decorrelation is not possible since it would also inadvertently remove the candidate signal from the data, therefore, we chose to model the raw RVs but with the incorporation of different noise and activity models.

Figure 6.

From top to bottom: the GLS periodogram for the RVs and activity indices from HARPS and FEROS: BIS, S-index, and CCF FWHM. The bottom panel shows the window function. Horizontal lines from top to bottom are the 0.1 per cent, 5 per cent, and 10 per cent false-alarm probabilities. The vertical red line marks the 4.1375-d planet candidate signal from the TESS data. The orange crosses show the most relevant peaks with their values in period written above them.

3.4 Global modelling

We jointly modelled the detrended TESS data with the HARPS and FEROS RVs, in order to obtain the mass, radius, and orbital parameters of TOI-179 b, with emperor.t, an upgraded version of emperor (Pena & Jenkins in prep) that handles both RVs and light curves.

emperor models transits with the transit light-curve calculation code pytransit (Parviainen 2015) using the quadratic limb darkening law from Mandel & Agol (2002), along with the triangular limb-darkening parametrization proposed by Kipping (2016). The transit model also includes an offset and white noise term, as well as a dilution factor (Espinoza et al. 2019), which we fixed to unity due to the expected dilution being |$\lt 1{{\ \rm per\ cent}}$| (de Leon et al. 2022, in preparation). RVs are modelled with the Keplerian model, including a velocity offset |$\gamma$| and an instrumental white noise parameter |$\sigma$| for each instrument and a first-order acceleration term. emperor implements instrumental q order autoregressive and p order moving average models (AR(q)MA(p); Tuomi et al. 2013) as a correlated noise model that deals with noise from, for example, stellar variability. Additionally, emperor has incorporated a linear stellar activity correlation model, following the description in Anglada-Escudé et al. (2016). emperor also allows for the modelling of the stellar density, which has been shown to improve parameter estimation as it deals with the degeneracy between the radius ratio and the impact parameter (Sandford et al. 2019; Vines et al. 2019). The algorithm assumes a Gaussian generative model and thus both the RV jitter and photometric white noise are added in quadrature to their respective observational uncertainties. The described model is given by

$$\begin{eqnarray*} rv_{i, \rm {ins}} &=& k(t_i) + ar_{i, \rm {ins}}(q) + ma_{i, \rm {ins}}(p) + A_{i, \rm {ins}} + \gamma _{\rm ins} + \dot{\gamma }_{\rm ins} \nonumber \\ &&\,+ \sigma _{i, \rm ins}+ \sigma _{i, \rm ins} \end{eqnarray*} $$

(2)

$$\begin{eqnarray*} ar_{i, \rm {ins}}(q) = \sum _{j=1}^q\phi _{j, \rm ins}\exp \left(\frac{t_{i-j} - t_i}{\alpha _{j, \rm ins}}\right)rv_{i-j, \rm ins} \end{eqnarray*} $$

(3)

$$\begin{eqnarray*} ma_{i, \rm {ins}}(p) = \sum _{j=1}^q\omega _{j, \rm ins}\exp \left(\frac{t_{i-j} - t_i}{\beta _{j, \rm ins}}\right)\epsilon _{i-j, \rm ins} \end{eqnarray*} $$

(4)

$$\begin{eqnarray*} A_{i, \rm ins} = \sum _\xi C_{\xi , \rm ins}\xi _{i, \rm ins} \end{eqnarray*} $$

(5)

Where |$k(t_i)$| is the Keplerian model evaluated at time |$t_i$|⁠, |$\dot{\gamma }_{\rm ins}$| is an acceleration term for each instrument, |$\gamma _{\rm ins}$| is the velocity offset of instrument ins, and |$ar_{i, \rm {ins}}(q)$| and |$ma_{i, \rm {ins}}(p)$| are the qth and pth order autoregression and moving average, respectively, with |$\phi$| and |$\omega$| being the strength of the correlated noise, |$\alpha$| and |$\beta$| are the characteristic time associated to the correlated noise, and |$\epsilon _{i-j, \rm ins}$| are the residuals of the model so far. |$A_{i, \rm {ins}}$| is the linear activity indices correlation model where |$\xi$| denotes the activity indices for each instrument and |$C_{\xi , \rm ins}$| are the activity coefficients. Finally, |$\sigma _{i, \rm ins}$| is the white noise term added for each instrument.

Sampling is carried out through emcee version 2.2.1 and its Parallel Tempering MCMC module, which is capable of sampling multimodal phase spaces without getting stuck in local maxima. This is achieved by modifying the posterior space, as shown in equation (6), where |$\beta = 1/T$| with T a temperature, |$\mathcal {L}$| the original likelihood, and |$\Pi$| is the prior density. This effectively flattens the posterior space with increasing temperature, allowing high temperature chains to sample the posterior more easily.

$$\begin{eqnarray*} \mathcal {L}_\beta (x) = \mathcal {L}^\beta (x)\Pi (x) \end{eqnarray*} $$

(6)

Given the active nature of TOI-179 and the significant amount of scatter in the data, we decided to test six different configurations for modelling the system: an AR(1)MA(1) model; activity correlations with the BIS; AR(1)MA(1) plus activities; AR(1)MA(0) model; AR(0)MA(1) model; and a pure Keplerian model, hereafter we will refer to those as runs 1, 2, 3, 4, 5, and 6. Additionally, we ran a noise model with activity correlations as a baseline for the model comparison. For the sampling, we used five temperatures, 1500 walkers, and 6000 steps, with a burn in period of half the number of steps, totaling 45 million samples for the search phase, and 1500 walkers with 6000 steps for the sampling stage, ending the fitting process with 9 million samples for parameter inference. We chose the temperature ladder to decrease as |$1/\sqrt{5}^i$| where |$i=0, 1, 2,... T$| and T the number of temperatures minus one. We used the same configuration for the walkers, temperatures and priors for all configurations mentioned previously.

3.4.1 Prior selection

We chose the period prior to be a Jeffreys prior due to it being uninformative in period space, and we chose the bounds to be from 0.1 to 6 d in order to avoid possible aliases at 8 d. We applied the transformation |$s = \sqrt{e}\sin (\omega)$|⁠, |$c = \sqrt{e}\cos (\omega)$|⁠, both bounded from −1 and 1 to fit for the eccentricity, e, and argument of periastron, |$\omega$|⁠. At each step of the MCMC sampling emperor makes sure e and |$\omega$| are physically plausible, i.e. |$0\le e\lt 1$| and |$0\le \omega \le 2\pi$|⁠. Additionally, an external prior is applied to the eccentricity in order to penalize high values while still allowing them if the data argue for it. The rest of the planetary parameters have uninformative, flat priors. For the stellar density, we chose a prior drawn from the mass and radius estimate derived with ariadne using the symmetrized 3|$\sigma$| credible interval values as the standard deviation to allow for a more conservative prior. We chose uninformative priors for the instrumental parameters as well, with the exception of the RV jitter, where a Normal prior was chosen in order to penalize high values, and the photometric offset, which also has a Normal prior centred around zero since there is no significant offset in the detrended light curve. We summarize the priors in Table 7.

Table 7.

Prior choices used in this work.

Parameter	Prior
Orbital parameters
P (d)	|$\ln \mathcal {U}(0.1, 6)$|
K (ms|$^{-1}$|⁠)	|$\mathcal {U}(0, 50)$|
T|$_c$| (JD)	|$\mathcal {U}(min(t^\dagger), max(t))$|
|$\sqrt{e}\sin (\omega)$|	|$\mathcal {U}(-1, 1)$|
|$\sqrt{e}\cos (\omega)$|	|$\mathcal {U}(-1, 1)$|
|$\omega$| (rads)	|$\mathcal {U}(0, 2\pi)$|
e	|$\mathcal {N}(0, 0.3^2)$|
|$R_p/R_{*}$|	|$\mathcal {U}(0.01, 0.5)$|
b	|$\mathcal {U}(0, 1)$|
Stellar parameters
|$\rho _{*}$| (gcm|$^{-3}$|⁠)	|$\mathcal {N}(2.241, 0.479^2)$|
|$C^\ddagger _\xi$|	|$\mathcal {U}(-max(C_\xi), max(C_\upsilon))$|
|$q^{*}_1$|	|$\mathcal {U}(0, 1)$|
|$q^{*}_2$|	|$\mathcal {U}(0, 1)$|
RV noise parameters
|$\gamma$| (ms|$^{-1}$|⁠)	|$\mathcal {U}(0, 3max(|rv|))$|
|$\sigma$| (ms|$^{-1}$|⁠)	|$\mathcal {N}(5, 5^2)$|
MA coefficient |$\omega$|	|$\mathcal {U}(-1, 1)$|
MA time-scale |$\beta$| (d)	|$\mathcal {U}(0, 10)$|
AR coefficient |$\phi$|	|$\mathcal {U}(-1, 1)$|
AR time-scale |$\alpha$| (d)	|$\mathcal {U}(0, 10)$|
Transit noise parameters
Offset (ppm)	|$\mathcal {N}(0, 0.1^2)$|
Jitter (ppm)	|$\ln {\mathcal {U}}(0.1, 10000)$|
dilution	fixed (1)
Acceleration parameter
|$\dot{\gamma }$| (ms|$^{-1}$|/yr)	|$\mathcal {U}(-1, 1)$|

|$^\dagger t=$|Time baseline of TESS data.

|$^\ddagger$| Activity indices were mean subtracted and normalized to their RMS.

|$^{*}$| Kipping LD parameters.

3.5 Model selection

We use the log posterior probability, Bayes Information Criterion (BIC) and Akaike Information Criterion (AIC), defined in equations (7) and (8), respectively, where n is the number of data points, k the number of parameters, and |$\mathcal {L}$| the maximum likelihood of the model, to do model comparison, using the log posterior to select the final model. The most probable model is Run 2, followed by Run 6, which suggests that the data are dominated by white noise instead of correlated noise. The semi-amplitude of all runs is consistent, with the exception of Run 1, which produces a smaller semi-amplitude by a factor of 2, but still within the 3|$\sigma$| credability intervals of the other runs.

$$\begin{eqnarray*} {\rm BIC} = k\ln (n) - 2\ln \left(\mathcal {L}\right) \end{eqnarray*} $$

(7)

$$\begin{eqnarray*} {\rm AIC} = 2k - 2\ln \left(\mathcal {L}\right) \end{eqnarray*} $$

(8)

We also did a seventh run that consists of the most probable configuration but with the eccentricity fixed to zero. This run had a resulting log posterior probability of 392697.59, significantly lower than the 392716.16 from Run 2. It is well-known that planets in such close orbits tend to show low eccentricity orbits, or even circular orbits, but since this signal is buried deep within the noise, it is not possible to fully constrain the orbital eccentricity, resulting in an eccentric orbit.

As an additional check, we ran the same configurations including the Minerva data and found that in all three runs, while the signal is still recovered, significant noise is added to the fit, resulting in lower posterior probabilities. Thus we choose to exclude the Minerva RVs from the final analysis. Table 8 show the aforementioned indicators for each run.

Table 8.

Model statistics for each run, compared against the baseline run. The row in bold face indicates the adopted run. The baseline run is a linear fit to the data.

run	Posterior	BIC	AIC
Run 1	25.89	100.38	9.79
Run 2	35.48	|$-$|13.86	|$-$|22.92
Run 3	27.45	128.78	11.01
Run 4	28.25	32.93	|$-$|3.30
Run 5	25.35	|$-$|197.23	219.49
Run 6	23.45	38.21	1.98
baseline	0	0	0

3.5.1 HD 18599 b

The modelling results show a radius ratio of |$R_\text{p}/R_*=0.0311\pm 0.0008$|⁠, resulting in a radius of 2.70 |$\pm$| 0.05 |$R_{\oplus }$|⁠, making HD 18599 b a sub-Neptune planet, while the semi-amplitude, K, of 11 |$\pm$| 3 m s|$^{-1}$|⁠, resulting in a derived mass of M|$_\text{p}=25.5\pm 4.6$| |$M_{\oplus }$|⁠, meaning HD 18599 b has a density of |$\rho _\text{p}=7.1\pm 1.4$| g cm|$^{-3}$|⁠. Such a density therefore is consistent with the planet maintaining an atmosphere that has a 24 per cent H|$_2$|O composition (see Section 4). We summarize the model parameters in Table 9 and show the phase folded RVs with the best-fitting model and 1–3|$\sigma$| credible intervals in Fig. 7, and the phase folded light curve with best-fitting model and 1-3σ credible intervals in Fig. 8.

Figure 7.

HD 18599 RVs folded at 4.137437 d. The black line represents the best-fitting model while the purple-shaded regions show the 1α, 2α, and 3|$\sigma$| confidence regions of the model. The lower panel shows the residuals of the fit. Circles are HARPS_pre data, upside down triangles show HARPS_post data, and upright triangles show the FEROS data. Green points show the data binned to 10 points in phase space. The colourbar encodes the observing time of each observation.

Figure 8.

TESS light curve folded at 4.137437 d in light blue points with pink points showing the light curved binned to 50 points in phase space. The dark blue line represents the best-fitting model and the pink shaded regions show the 1, 2, and 3|$\sigma$| confidence regions of the model. The lower panel shows the residuals of the fit in ppm.

Table 9.

Fitted and derived parameters.

Parameter	Prior
Orbital parameters
Fitted
P (d)	4.137 437 |$\pm$| 0.000 004
K (ms|$^{-1}$|⁠)	11 |$\pm$| 2
T|$_c$| (JD)	2458726.9576 |$\pm$| 0.0004
|$\sqrt{e}\sin (\omega)$|	|$-$|0.1 |$\pm$| 0.3
|$\sqrt{e}\cos (\omega)$|	0.5|$^{+0.1}_{-0.3}$|
|$R_p/R_{*}$|	0.0311 |$\pm$| 0.0008
b	0.58 |$\pm$| 0.11
Derived
|$\omega$| (rads)	0.2|$^{+0.6}_{-0.2}$|
e	0.2|$^{+0.1}_{-0.2}$|
i (deg)	87.7|$^{+0.5}_{-0.7}$|
|${\rm M_p}$| (⁠|$M_{\oplus }$|⁠)	25.5 |$\pm$| 4.6
|${\rm R_p}$| (⁠|$R_{\oplus }$|⁠)	2.70 |$\pm$| 0.05
|$\rho _\text{p}$| (g cm|$^{-3}$|⁠)	7.1 |$\pm$| 1.4
T|$_p$| (JD)	2458726.41|$^{+0.27}_{-0.23}$|
|$a/R_{*}$|	13.78|$^{+0.96}_{-1.08}$|
a (au)	0.05 |$\pm$| 0.02
T|$_\text{eq}$| (K)	863|$^{+21}_{-18}$|
Insolation (⁠|$S_\oplus$|⁠)	145|$^{+16}_{-12}$|
Stellar parameters
|$\rho _{*}$| (g cm|$^{-3}$|⁠)	2.87|$^{+0.62}_{-0.66}$|
|$C_\text{FEROS}$|	0.2 |$\pm$| 0.8
|$C_\text{HARPS_pre}$|	1.0|$^{+0.3}_{-0.8}$|
|$C_\text{HARPS_post}$|	|$-2.1^{+0.7}_{-0.0}$|
|$q_1$|	0.68 |$\pm$| 0.17
|$q_2$|	0.32 |$\pm$| 0.07
RV noise parameters
|$\gamma _\text{FEROS}$| (m s|$^{-1}$|⁠)	|$-86.3^{+16.3}_{-14.5}$|
|$\sigma _\text{FEROS}$| (m s|$^{-1}$|⁠)	10.2|$^{+5.0}_{-4.4}$|
|$\gamma _\text{HARPS_pre}$| (m s|$^{-1}$|⁠)	|$-2.7^{+9.8}_{-9.6}$|
|$\sigma _\text{HARPS_pre}$| (m s|$^{-1}$|⁠)	15.4|$^{+4.5}_{-3.6}$|
|$\gamma _\text{HARPS_post}$| (m s|$^{-1}$|⁠)	|$-68.9^{+11.2}_{-11.7}$|
MA coefficient |$\omega$|	|$-$|-
MA Time-scale |$\beta$| (d)	|$-$|-
AR coefficient |$\phi$|	|$-$|-
AR time-scale |$\alpha$| (d)	|$-$|-
Transit noise parameters
offset (ppm)	|$-$|0.000009 |$\pm$| 0.000003
jitter (ppm)	291 |$\pm$| 4
dilution	1 (fixed)
Acceleration parameter
|$\dot{\gamma }$| (ms|$^{-1}$|/yr)	17 |$\pm$| 3

The ARMA correlated noise model was not included in the selected model as described in the text.

We searched for additional signals in the RV residuals, but nothing significant was found. More sophisticated activity models are needed to further disentangle Keplerian signals from activity induced ones in noisy RV data.

3.5.2 TTV search

The Kepler mission was responsible for the detection of a few thousand exoplanets, from which a handful of multiplanet systems were validated through the transit timing variation (TTV) technique (Holczer et al. 2016). Since the majority of Kepler stars were too faint for precise RV measurements, TTVs became key to derive planetary masses in many cases, thus validating the planetary systems (Lithwick, Xie & Wu 2012). Moreover, this technique supported the idea that hot Jupiters are not part of multiplanet systems (Steffen et al. 2012), therefore setting major constraints on planetary migration models of giant planets. The TTV method (Agol et al. 2005) relies on the comparison between the measured and expected central transit times (⁠|$T_c$|⁠) from a linear model given by |$T_n = T_0 + n \times P$|⁠, where |$T_n$|⁠, P, and n are the central transits, period, and transit number (⁠|$n=1,2,3,...$|⁠), respectively. Therefore, a deviation from the linear model may indicate the occurrence of dynamical interactions, where the most frequent cases are mean motion resonance planet–planet interaction and planet–star interaction leading to orbital tidal decays (Yee et al. 2019).

The TTV analysis was performed using allesfitter⁷ (Günther & Daylan 2020) on the detrended light curve. Stellar and planetary parameters were fixed to the median posteriors from Tables 3 and 9, respectively, except for |$T_n$| which used |$\mathcal {U}(T_n - 0.05, T_n + 0.05)$| for each nth transit. The O-C transit timing plot (Fig. 9) agrees with the linear model thus indicating no dynamical interaction with a second body.

Figure 9.

Observed minus computed mid-transit times for HD 18599b. The residuals (TTV) are shown considering the proposed ephemerides. The dashed cyan line shows zero variation and the shaded cyan region shows the 1|$\sigma$| uncertainties on the linear ephemeris from EMPEROR. The plot shows no significant variation as all observations are within 1|$\sigma$| of the best-fitting ephemerides.

4 DISCUSSION

We have confirmed the transit parameters found by de Leon et al. (2022) and have characterized the orbit of the planet. Fig. 10 places HD 18599 b in the period-radius diagram, where we find it at the edge the Neptune desert. While the edge of the desert is fairly populated, the figure shows that HD 18599 b is the youngest mini-Neptune discovered to date within this region, making this planet a key addition to the population.

Figure 10.

Period-radius diagram for planets with radii measured to a precision better than 5 per cent. Relevant discoveries are named, and HD 18599 is shown with a big yellow star at the edge of the desert. The pink shaded area shows the Neptune desert, with the edges defined by Mazeh et al. (2016). The colours show the ages of the systems and black shows systems without an estimated age. We note that ages have not all been measured in an uniform way.

4.1 Possible compositions

We measured the density of HD 18599 b to be 7.1 |$\pm$| 1.4 g cm|$^{-3}$|⁠, and therefore we can compare this value to other known planets with similar masses and radii by placing it on a mass-radius diagram (Fig. 11). HD 18599 b is shown by the red star, and it is placed against several two-layer composition models from Zeng, Sasselov & Jacobsen (2016) and five two-layer envelope models from Zeng et al. (2019). Its position in the mass-radius diagram suggests the planet can be composed of 23 per cent H|$_2$|O, and 77 per cent rock (MgSiO|$_3$|⁠) and iron core, while the 1|$\sigma$| uncertainties allow the composition to vary from 10 per cent to 41 per cent H|$_2$|O.

Figure 11.

Mass-radius diagram for planets with masses and radii measured with a precision better than 20 per cent (grey circles) in the range R|$_{\rm p}\lt 5\,$|R|$_\oplus$| and M|$_{\rm p}\lt 30\,$|M|$_\oplus$|⁠, retrieved from the transiting exoplanets catalog TEPCat (Southworth 2011). Black points show relevant discoveries. HD 18599 b is shown with a red star. Solid, coloured lines show different bi-layer composition models from Zeng et al. (2016), ranging from 100 per cent iron core planet to 100 per cent H|$_2$|O planet. Also, five two-layer models from Zeng et al. (2019) are shown for 2 per cent H|$_2$| envelopes at different temperatures and core compositions.

Due to the inherent degeneracies in three-layer bulk composition models with H/He envelopes on the sub-Neptune regime (⁠|$2\le$| R |$\le 4$| R|$_\oplus$|⁠; Lopez & Fortney 2014), we chose to further investigate the internal composition of HD 18599 b using the public tool smint⁸ (Structure Model INTerpolator) and its interpolation and envelope mass fraction fitting package, first introduced by Piaulet et al. (2021). This code uses models from Lopez & Fortney (2014), Zeng et al. (2016), and Aguichine et al. (2021), along with emcee to perform an MCMC fit of the H|$_2$|O or H/He mass fraction based on the planetary mass, radius, age, and insolation flux. We set up the MCMC with 100 walkers and 10 000 steps for all of the analyses. First, we studied the H|$_2$|O mass fraction (WMF) using the Aguichine et al. (2021) irradiated ocean worlds mass–radius relationships, where we found a WMF of 0.3 |$\pm$| 0.1. The iron core mass fraction remains widely unconstrained with a value of 0.4 |$\pm$| 0.3. The gas-to-core mass ratio analysis show that HD 18599 b has a 1|$^{+14}_{-1}$| per cent mass fraction in H/He. The lack of an H/He envelope and the amount of water in the bulk composition argues for a steam atmosphere surrounding the planet.

The high-mass regime in which HD 18599 b is found is poorly sampled, making the planet an ideal window into the physical and dynamical processes behind the evolution and formation of young Neptunes and their interactions with active stars. Further work to better characterize the atmosphere of the planet are warranted, which would allow a better understanding of the internal composition of the planet and the physics of young planetary atmospheres. Additional precision RV campaigns would also allow the possibility to discover additional companions in the system, in particular since it is known that small planet systems commonly harbour more than one planet (for example Scarsdale et al. 2021; Vach et al. 2021).

4.2 HD 18599 as a stellar activity model testbed

We tested different correlated noise models in order to account for the stellar activity signals appearing in the RV measurements. These models included AR, MA, linear activity correlations, and a combination of all of these. We found that in this case, linear activity correlations gave rise to a higher posterior probability model, while ARMA models fell short of explaining the correlated noise, likely mainly due to the low observational cadence, since these models perform best under high cadence conditions. HD 18599 is an excellent testbed to apply more sophisticated activity models, also further encouraging more RV follow up work to be pursued.

5 CONCLUSIONS

Using TESS photometry and RVs from HARPS and FEROS, we were able to characterize a dense planet orbiting the young star HD 18599. The world is found to have an orbital period of 4.137437 |$\pm$| 0.000004 d, radius of 2.70 |$\pm$| 0.05 |$R_{\oplus }$|⁠, (making it a mini-Neptune planet in size), and unusually high mass of 25.5 |$\pm$| 4.6 |$M_{\oplus }$| resulting in a planet with a density of 7.1 |$\pm$| 1.4 g cm|$^{-3}$|⁠. The mass-radius diagram and bulk composition models suggests an ice rich core with a range of 10–41 per cent H|$_2$|O contribution, and a 1|$^{+13}_{-1}$| per cent H/He envelope mass fraction. Further analysis shows a WMF of 33 |$\pm$| 10 per cent and an unconstrained iron core mass fraction ranging from 10 per cent to 70 per cent. HD 18599 b is the first mini-Neptune to orbit such a young (< 300 Myr) and bright (V|$\sim$|8.9) star, making it an ideal system to dive deep into evolutionary studies and atmospheric characterization.

SUPPORTING INFORMATION

Table 1. Detrended TESS photometry for HD 18599.

Table 2. RVs and Activity Indices for HD 18599.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

MINERVA-Australis is supported by Australian Research Council LIEF Grant LE160100001, Discovery Grants DP180100972 and DP220100365, Mount Cuba Astronomical Foundation, and institutional partners University of Southern Queensland, UNSW Sydney, MIT, Nanjing University, George Mason University, University of Louisville, University of California Riverside, University of Florida, and The University of Texas at Austin.

We respectfully acknowledge the traditional custodians of all lands throughout Australia, and recognize their continued cultural and spiritual connection to the land, waterways, cosmos, and community. We pay our deepest respects to all elders, ancestors, and descendants of the Giabal, Jarowair, and Kambuwal nations, upon whose lands the Minerva-Australis facility at Mt Kent is situated.

JIV acknowledges support of CONICYT-PFCHA/Doctorado Nacional-21191829. JSJ greatfully acknowledges support by FONDECYT grant 1201371 and from the ANID BASAL projects ACE210002 and FB210003. HRAJ acknowledges support from STFC grant ST/T007311/1 This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

DATA AVAILABILITY

The data underlying this article are available in the article and in its online supplementary material.

Footnotes

1

https://exofop.ipac.caltech.edu/tess/target.php?id = 207141131

2

https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html

3

ESO programs 192.C-0224, 0102.C-0525, 0102.D-0483, 0101.C-0510, 0102.C-0451, 0103.C-0759.

4

https://github.com/jvines/Ceres-plusplus

5

https://github.com/jvines/astroARIADNE

6

https://juliet.readthedocs.io/en/latest/

7

https://www.allesfitter.com/

8

https://github.com/cpiaulet/smint