© The Authors 2023

1 Introduction

Thanks to the constantly increasing number of characterised extrasolar planets, there is a growing number of studies related to the exoplanet demographics, defined by Gaudi et al. (2021) as a research area mainly dedicated to study planets' distribution as a function of the physical parameters that may influence planet formation and evolution, exploring a broad range of the parameter space. Among the many physical parameters involved, a particularly crucial one is the age of the exoplanet systems. Drawing a relationship between the observed properties of the systems and their ages is a necessary step to take in order to build a general picture of how planets form and evolve with time, and to constrain the timescales of the main processes at work within the first hundreds of millions years from planet formation. That is why it is crucial to detect planets at different stages of their evolution (infant: age <20 Myr; young: 20< age <100 Myr; intermediate-age: up to 800 Myr), and determine their orbital and main physical parameters. Among several key aspects, characterising young and close-in super-Earths and sub-Neptunes spotted at different ages offers the opportunity of theoretically studying how the atmospheric mass-loss rate evolves with time, in response to the high stellar high-energy irradiation, after dissipation of the protoplanetary disk.

Currently, the statistical sample of exoplanet systems younger than ~800 Myr is still small if compared to the total number of discoveries achieved so far. This prevents drawing results from comparative analyses with more evolved systems. Detecting infant and young planets is particularly challenging, especially for blind searches carried out with the radial velocity (RV) method because of the very high level of scatter (up to several hundreds of ms⁻¹) mostly due to stellar magnetic activity, which heavily hampers the detection of periodic signals produced by planetary companions. The photometric transit method proved to be a more fruitful detection technique in this case, in particular thanks to the observations of the Kepler/K2 and the Transiting Exoplanet Survey Satellite (TESS) space telescopes. They yielded discoveries of planetary-sized companions, which triggered RV follow-up campaigns to measure their masses and bulk densities. Among these, the long-term Italian project Global Architecture of Planetary Systems (GAPS; Covino et al. 2013; Poretti et al. 2016) is investing a good fraction of the allocated observing time to follow-up stars with an age between ~2–600 Myr using the High Accuracy Radial velocity Planet Searcher (HARPS-N) spectrograph (Cosentino et al. 2012) located in the Northern Hemisphere at Telescopio Nazionale Galileo (TNG). After an initial effort in revising previously announced and debated planet detections (e.g. Carleo et al. 2018, 2020; Damasso et al. 2020), the main focus of the GAPS young planets sub-programme shifted to the measurement of the fundamental properties of transiting planets. GAPS has observed and characterised several systems so far (e.g. Carleo et al. 2021; Suárez Mascareño et al. 2021; Nardiello et al. 2022), and further characterisation studies have been carried out in parallel in the Southern Hemisphere with the HARPS spectrograph (e.g. Benatti et al. 2021; Desidera et al. 2023). The measurement of the planetary dynamical masses remains crucial to investigate the evolution of exo-atmospheres even when it is only possible to provide upper limits, and the current escape rates can be assessed through atmospheric evaporation models (e.g. Benatti et al. 2021; Carleo et al. 2021; Poppenhaeger et al. 2021; Maggio et al. 2022).

In this paper, we present a characterisation study of the two-planet system detected by the TESS space telescope (Ricker et al. 2016) around the bright and 400-Myr old Sun-like star HD 63433, a member of the Ursa Major moving group, also known as TESS object of interest (TOI) 1726. The system was first validated and partly characterised by Mann et al. (2020). In particular, they revealed that the orbit of the innermost planet is prograde, and emphasised that HD 63433 represents a benchmark system for investigating the planet's evolution after the first hundreds millions years. The first results for this system led the GAPS collaboration to select HD 63433 as a promising target for an intense RV follow-up for further characterisation. The architecture of HD 63433 was also investigated by Dai et al. (2020), who measured the sky-projected spin-orbit angle of planet c, consistent with an aligned and prograde orbit. Results from Mann et al. (2020) and Dai et al. (2020) indicate that the system apparently did not experience catastrophic events capable of heavily influencing the dynamical evolution. The presence of a cold Jupiter in the system was examined by Hirsch et al. (2021), who analysed sparse RVs with a long time baseline, and found evidence only for a signal clearly ascribable to stellar activity.

The main purpose of the GAPS follow-up campaign of HD 63433 was measuring the masses of the two known planets. Accomplishing this goal becomes even more compelling to interpret the recent observational results of Zhang et al. (2022). They found evidence of Lyα absorption during a transit of HD 63433 c, suggesting ongoing atmospheric evaporation from this planet, while they suggest that HD 63433 b has already lost its primordial atmosphere: these observations have to be necessarily reconciled with mass measurements within the framework of current theoretical models. The study by Zhang et al. (2022) emphasises a key aspect: the two planets represent a great opportunity for a comparative study within the same system, to examine how their different orbital configurations led to the apparently diverse atmospheric mass-loss timescales.

In this work, we focus on the planetary mass measurements through the exploitation of RVs calculated from HARPS-N spectra, and we investigate the time evolution of the planetary atmospheres through evaporation models and updated stellar evolutionary tracks. The paper is organised as follows. In Sect. 2, we describe the photometric and spectroscopic datasets, and in Sect. 3 we present our measured fundamental stellar parameters. In Sect. 4, we present a frequency content analysis of the spectroscopic dataset, as a preparatory step to a more sophisticated analysis discussed in Sect. 5. In Sect. 6, we use the mass-radius diagram for a general comment on the possible planetary structures constrained by our results, especially in the context of known young exoplanets. Our results for planets b and c are then used for investigating the atmospheric evolutionary history in Sect. 7. Conclusions are summarised in Sect. 8.

2 Observations and data reduction

2.1 TESS photometry

TESS observed HD 63433 in Sector 20 (from 24 December 2019 to 21 January 2020) during Cycle 2 (GO22038, PI: Roettenbacher; GO22203, PI: Ge; GO22032, PI: Metcalfe), and again in Sectors 44–47 (from 12 October 2021 to 28 January 2022) during Cycle 4 (GO4104, PI: Huber; GO4242, PI: Mayo; GO4191, PI: Burt; GO4060, PI: ; GO4039, PI: Mann). In our study, we used the short-cadence (2 min) light curves. We did not make use of the official Presearch Data Conditioning Simple Aperture Photometry (PDCSAP, Smith et al. 2012; Stumpe et al. 2012, 2014) light curves, because they are affected by systematics due to overcorrections and/or injection of spurious signals. We corrected the Simple Aperture Photometry (SAP) data by using cotrending basis vectors obtained with the tools developed by Nardiello et al. (2020). The extracted light curve for all Sectors is shown in Fig. 1. In order to model and remove the clearly visible variability due to the stellar activity, we used the package PyORBIT (Malavolta 2016; Malavolta et al. 2016, 2018). First, we masked all the transits from the light curve, then we modelled the stellar activity by using Gaussian rocess (GP) regression performed with the package celerite2 (Foreman-Mackey et al. 2017; Foreman-Mackey 2018). We interpolated each excluded point within the masked transits using the best-fit GP solution, and we corrected the light curve variability by dividing the observed flux for the best-fit model. The final result is a flattened light curve with all the transits preserved, which we used in the combined light curve-RV analysis (Sect. 5). We already verified in previous works (see, e.g. Nardiello 2020; Nardiello et al. 2021) that this procedure for detrending the light curve does not alter the transit signal, by affecting its shape and depth. At the time of writing, according to the Web TESS Viewing Tool HD 63433 won't be observed by TESS anymore up to September 2023 (Cycle 5).

Fig. 1 TESS light curve of HD 63433 for Sector 20 (upper panel), and Sectors from 44 to 47 (lower panel). Planetary transits and bad-quality points have been removed from the dataset. The red curve represents the best-fit GP model.

2.2 STELLA photometry

HD 63433 was observed with the STELLA telescope and its wide-field imager WiFSIP (Strassmeier et al. 2004) between 29 September 2020 and 17 March 2021 in Johnson V and Cousin I bands (82 and 81 epochs, respectively). The follow-up was intended to monitor the photometric evolution in response to a changing activity during one of the spectroscopic observing seasons. Each nightly observing block consisted of seven exposures with an exposure time of 1 second for each band. To avoid saturation of this bright target (V_T ~ 7 mag), the telescope was defocussed to get about 5 arcsec for the full width at half maximum (FWHM) of the object point spread function (PSF). The data reduction followed that detailed description in Mallonn & Strassmeier (2016). We averaged the nightly exposures for each filter. The light curves are shown in Fig. C.1 (left column). The scatter and median uncertainty for the V and I-band light curve are 0.009 mag and 0.01 mag, respectively.

2.3 HARPS-N spectroscopic data

The spectroscopic follow-up of HD 63433 was carried out with the HARPS-N spectrograph (Cosentino et al. 2012) mounted at the Telescopio Nazionale Galileo (TNG) on the island of La Palma (Canary Islands, Spain). The observations consist of 103 spectra collected between 26 February 2020 and 17 April 2022 (781 days), with a typical exposure time of 600 s, and a median S/N of 203 measured at a wavelength of ~550 nm. The spectra have been reduced with the standard Data Reduction Software (DRS) pipeline (version 3.7.1) through the YABI workflow interface (Hunter et al. 2012), which is maintained by the Italian centre for Astronomical Archive (IA2)¹. The RVs and activity diagnostics such as the FWHM and bisector velocity span (BIS) have been derived from the DRS cross-correlation function (CCF), which was calculated by adopting a reference mask fora star of spectral type G2, and a half-window width of 200 km s⁻¹. The RVs have a median internal error σ_{RV, DRS} = 0.81 m s⁻¹, and an rms of 23.3 m s⁻¹.

As a sanity check, we measured the RVs also using the Template Enhanced Radial velocity Reanalysis Application (TERRA) pipeline (v1.8; Anglada-Escudé & Butler 2012), which, for young and active stars, is in principle more effective in providing RVs less affected by stellar activity than those calculated from the CCF. TERRA dataset is characterised by a median internal error σ_RV,TERRA = 1.25 m s⁻¹, and an rms of 24.1 m s⁻¹. In Sect. 5, we investigate the results obtained using both DRS and TERRA RVs.

We calculated two additional chromospheric activity diagnostics from the HARPS-N spectra, namely the log ${R^{\prime }}_{\text{HK}}$ from the Call H&K lines (provided by the DRS), and the index based on the H-alpha line calculated with the tool ACTIN (v1.3.9; Gomes da Silva et al. 2018). RVs and activity diagnostics are listed in Tables D.1 and D.2.

We note that there is very little overlap between the beginning of the last season of the RVs and the final part of the TESS light curve of Sector 47. The lack of simultaneous photometric monitoring with the RVs makes it impossible to use the light curve as an ancillary dataset to filter out the activity term in the RVs. That is especially true for a star such as HD 63433, which shows rapidly changing photospheric active regions in the light curve that will generate RV variations independent of any planet. The light curve observed by STELLA is characterised by sparse sampling, and overlaps only partially with the second season covered by HARPS-N. Therefore, even the STELLA dataset is not of practical use to effectively correct stellar activity in the RVs.

3 Fundamental stellar parameters

HD 63433 is a dwarf star that belongs to the Ursa Major moving group, with an estimated age of roughly 400 Myr (see Mann et al. 2020). Given its relatively young age, the spectroscopic analysis might be hampered by the intense magnetic fields that alter the structure of the stellar photosphere, in particular the upper layers. This alteration affects the formation of spectral lines in these layers for which abundances show a trend with optical depths (see Baratella et al. 2020b for further details). As a consequence, the standard spectroscopic approach, that is the equivalent width (EW) method, based on the use of iron (Fe) lines, could fail. In particular, when deriving the microturbulence velocity (ξ) by imposing that weak and strong iron lines have the same abundance, ξ could be over-estimated, which leads to an under-estimation of the iron abundance ([Fe/H]).

Following the same strategy and the same line list as in Baratella et al. (2020a), we applied a new method that consists of using a combination of Fe and titanium (Ti) lines to derive T_eff (through excitation equilibrium), and using only Ti lines to derive log 𝑔 (through ionisation equilibrium) and ξ (by zeroing the trend between individual abundances and strength, or EW, of the lines). The code MOOG (Sneden 1973) was used to analyse the HARPS-N co-added spectrum. We measured line EWs with the ARES v2 code (Sousa et al. 2015) and discarded lines with errors larger than 10% and with EW > 120 mÅ. The 1D LTE model atmospheres linearly interpolated from the ATLAS9 grid of Castelli & Kurucz (2003), with new opacities (ODFNEW), were adopted.

From the spectral type of the star (G2V), we estimated an input T_eff of 5660 K from the relation by Pecaut & Mamajek (2013). Using different colour indexes (the reddening is negligible as the star is 22 pc distant from the Sun) and different calibrated relations (Casagrande et al. 2010, Mucciarelli et al. 2021), the photometric estimates of T_eff range between 5589 K (from Bp – Rp), 5619 K (from V − K) and 5762 K (from J – K). The derived surface gravity based on the Gaia DR3 parallax (Gaia Collaboration 2022) is between 4.50 and 4.54dex (error of the order of 0.05), depending on the photometric T_eff considered. Finally, using the Dutra-Ferreira et al. (2016) relation, we estimated an initial ξ of 0.095 ± 0.05 km s⁻¹.

The final values obtained with the new spectroscopic approach are: T_eff=5700±75 K, log g=4.54±0.05 and ξ=1.03±0.10 km s⁻¹. The iron abundance is [Fe/H] I=0.02±0.06, [Fe/H] II=0.05±0.05, while the titanium abundance is [Ti/H] I = 0.04±0.07, [Ti/H] II=0.05±0.04. The uncertainties reported are the quadratic sum of the scatter due to the EW measurements and the contribution of the parameter uncertainties to the final abundances. Our results are in perfect agreement with the initial guesses and with the recent analysis of Mann et al. (2020).

Fixing T_eff, log g, ξ, and [Fe/H] I to the values found above, we measured the stellar projected rotational velocity (v sin i_★) using the same MOOG code and applying the spectral synthesis of two regions around 6200 and 6700 Å. We adopted the same grid of model atmosphere and, after fixing the macroturbulence velocity to the value of 3.2 km s⁻¹ from the relationship by Brewer et al. (2016), we find a v sin i_★ of 7.2 ± 0.7 km s⁻¹, consistent with the result by Mann et al. (2020).

Finally, we also derived the lithium abundance log A(Li)_NLTE from the measured lithium EW (=84.5±2.0 mÅ) and considering our stellar parameters previously derived together with the NLTE corrections by Lind et al. (2009). The values of the lithium abundance is 2.56 ± 0.06 and the position of the target in a log A(Li)_NLTE−T_eff diagram is compatible with the UMa group, as expected (see Mann et al. 2020).

We calculated the stellar radius and mass with the EXO-FASTv2 tool (Eastman et al. 2019), by fitting the stellar Spectral Energy Distribution (SED) and providing the stellar luminosity derived from the SED as input to the MIST stellar evolutionary tracks (Dotter 2016). For the SED we considered the Tycho B and V magnitudes (Høg et al. 2000), the 2MASS near-IR J, H and K magnitudes (Cutri et al. 2003), and the WISE mid-IR W1, W2, W3 and W4 magnitudes (Cutri et al. 2021). We imposed Gaussian priors on (i) the stellar effective temperature and metallicity from our analysis of the HARPS-N spectra; (ii) the parallax 44.6848 ± 0.0228 mas from the Gaia EDR3 (Gaia Collaboration 2016, Gaia Collaboration 2021) and (iii) the stellar age 413 ± 23 Myr from the most updated value for the Ursa Major association (Jones et al. 2015). We used uninformative priors for all the other parameters. The best fit of the SED is shown in Fig. 2. We found R_★ = 0.897 ± 0.019 R_⊙ and M_★ = 0.994 ± 0.028 M_⊙, in excellent agreement with the previous determination by Mann et al. (2020).

From modelling the activity of the TESS light curve (Fig. 1) with a GP regression, we obtained a rotation period P = 6.48 ± 0.08 days. This result is in agreement with the measurement by Mann et al. (2020), and it is confirmed by looking at the peri-odogram of the V-band STELLA light curve calculated with the generalised Lomb-Scargle (GLS, Zechmeister & Kürster 2009) tool (Fig. C.1), which shows the main and significant peak at 6.4 days. The fundamental stellar properties of HD 63433 are summarised in Table 1.

Fig. 2 Spectral energy distribution of the host star HD 63433 with the best-fit model overplotted (solid line). Red and blue points correspond to the observed and predicted values, respectively.

4 Frequency content analysis of radial velocities and activity diagnostics

We analysed the frequency content of RVs and activity indicators to investigate the presence of sinusoidal periodic signals. The time series of the RVs (both DRS and TERRA), and the corresponding GLS periodograms are shown in the panels of the first two columns of Fig. 3. The main peak is statistically significant (bootstrap false alarm probability FAP <1%) and located at the second harmonic of the stellar rotation period. The periodograms of the residuals, after removing the best-fit sinusoid calculated by GLS (third column of Fig. 3), show the main peak at the first harmonic of the rotation period, with low significance in the case of TERRA RVs. No other significant peaks are detected, especially at the orbital frequencies of the planets, indicating that the RVs scatter is dominated by stellar activity.

The time series and GLS periodograms of the spectroscopic activity diagnostics are shown in Fig. 4. The log R'_HK index and BIS contain a quite significant signal related to the stellar rotation period. The main peaks in the periodograms occur at a frequency corresponding to P_rot,★ and to the second harmonic of P_rot,★, respectively. A non-significant (FAP = 41 %) peak at a frequency compatible with the rotational frequency is also detected in the residuals of the H-alpha index, after removing a long-term curvature clearly visible in the time series, for which we cannot assess a periodicity, if actually present. The periodogram of the FWHM time series has the main peak located at ~162 days (FAP = 0.3%), whose origin is not clear. After a pre-whitening, the periodogram of the FWHM residuals shows a peak at 2.15 days with FAP = 51%. Given the very low statistical significance of this signal, which corresponds to the second harmonic of P_rot,★, in the following analysis we do not use the FWHM to correct the RVs for the stellar activity term.

We note that the RV and BIS data are significantly correlated. The Spearman's rank correlation coefficient is ρ = −0.88, calculated using the r_correlate function in IDL.

Fig. 3 Time series of the RVs calculated from HARPS-N spectra with the DRS and TERRA pipelines, and their GLS periodograms. First column: RV timeseries. The mean value has been subtracted from the original DRS data. Second column: The corresponding GLS periodograms of the two RV dataset. FAP levels are of 1% and 10% are indicated by yellow and red horizontal lines, respectively. They are calculated through a bootstrap analysis. Third column: GLS periodograms of the pre-whitened RV data.

5 Joint analysis of transits and RVs

We modelled the RV time series jointly with the detrended and flattened TESS light curve, which includes only the part of the data with the transits, and a sufficiently long out-of-transit baseline for a proper modelling of the planet ingress and egress. A key point in our analysis is the modelling of the stellar activity term which is the dominant signal in the RV time series. We took special care of it by testing different models, mostly based on the GP regression analysis. A schematic description of the test models is provided in Table 2, with more details given in Appendix A.

For all the tested models, we explored the full parameter space using the publicly available Monte Carlo (MC) nested sampler and Bayesian inference tool MultiNest v3.10 (e.g. Feroz et al. 2019), through the pyMultiNest wrapper (Buchner et al. 2014). Our MC set-up included 300 live points, a sampling efficiency of 0.5, and a Bayesian tolerance of 0.3. We used the code batman (Kreidberg 2015) for modelling the photometric transits. In a few cases, when the same dataset has been used, we performed a Bayesian model selection by comparing the values of the Bayesian evidence ln 𝒵 calculated by MultiNest, and using the empirical scale indicated in Table 1 of Feroz et al. (2011) to assess their relative statistical significance². We assigned the same a-priori probability to each model.

Concerning the light curve, we modelled the limb darkening with a quadratic law, and fitted the coefficients LD_c1 and LD_c2 using the formalism and uniform priors given by Kipping (2013, see Eqs. (15) and (16) therein). We also introduced constant jitters σ_jit,TESS added in quadrature to the nominal photometric uncertainties, one jitter term for Sector 20, and one for Sectors 44-47, to take into account the change in the TESS performance over nearly 2 yr. The complete list of priors used for all the test models is provided in Table A.1.

The main goals of our analysis are the improvement of the orbital parameters for HD 63433 b and HD 63433 c, and the measure of the fundamental planetary parameters radius, mass, and average density. Based on the statistical results and mass measurements summarised in Table 2, we notice that:

• the GP quasi-periodic and quasi-periodic with cosine (QPC) kernels perform equally well (models M1 and M8), and better than the double simple harmonic oscillator (dSHO) kernel (model M7), when modelling only the RV time series. It is not statistically advantageous treating each observing season separately (model M2). The stellar rotation period is recovered with high precision (GP hyper-parameter $\theta ={6.390}_{-0.005}^{+0.007}$ days for M1), even with an adopted uniform and large prior 𝒰(0,10) days;

• results for all the test models show that the planetary orbits can be assumed circular. The data do not allow for significant non-zero eccentricities (see Table 3);

• we cannot constrain the mass m_b of HD 63433 b, but only derive upper limits, with values that change depending on the model. We notice that the orbital period P_b and the stellar rotation period P_rot,★ are not dissimilar, and this could make more difficult the identification of a planetary Doppler signal with a small semi-amplitude, despite the very precise transit ephemeris recovered;

• most of the models allow us to retrieve the mass mc of HD 63433 c with a statistical significance (defined as the ratio $m_c/{\sigma }_{{\text{m}}_{\text{c}}}^{-}$) in the range 2.1−2.7σ (with the exception of models M6, M7, and M9, for which the significance is lower);

• using the more complex multi-dimensional GP (quasi-periodic kernel), we find lower upper limits for the masses, with no improvement in their precision.

Despite the fact that we followed several pathways to model the stellar activity in the RV dataset, the data do not allow for a significant and accurate measurement of the mass m_c, but we can conclude with high confidence that the mass of HD 63433 c is less than twice the mass of Neptune. Nonetheless, we note that in a few cases the median values of Gaussian-like posterior distributions are very similar, and suggest that m_c could be actually consistent with the mass of Neptune. Examples of posteriors for mb and mc are shown in Fig. 5, with reference to Table 2. We also note that, in a few cases, masses of young planets have been claimed with a precision of the order of 3−3.5σ . For example, Desidera et al. (2023) found that TOI-179 b (further discussed in Sect. 6) has a mass of ${24.1}_{-7.7}^{+7.1}{M}_{\oplus }$ (3.1σ precision). The host TOI-179 has a similar mass, and an activity-induced RV scatter very similar to that measured for HD 63433 (~ 24m s⁻¹), while TOI-179 b has an orbital period nearly five times shorter than HD 63433 c, making a more precise RV detection of the planet potentially easier. K2-100 b (Barragán et al. 2019b), is a 750 Myr-old planet with an orbital period of 1.7-d, and a measured mass of 21.8±6.2 M_⊕ (3.5σ precision). Similarly, mass measurements have been claimed for TOI-560b (480Myr; $m={10.2}_{-3.1}^{+3.4}{M}_{\oplus }$, i.e. 3.3σ precision; Barragán et al. 2022b), and for Kepler-411 d (200 Myr; m = 15.2 ± 5.1, i.e. 3σ precision, measured through transit timing variations; Sun et al. 2019). Since we determined the mass of HD 63433 c with a significance of 2.7σ at best, we do not consider this result significant enough for claiming a mass measurement, but indeed it deserves further attention to be confirmed with additional data.

Concerning the transit ephemeris and other specific parameters of the transit model, we did not find significant differences among the results obtained for all the tested models (except for the models M5 and M6, where the light curve is not fitted). We provide in Table 3 the results of model M4, selected among the solutions available, which are all equivalent concerning the transit ephemeris. We underline that the election of one of the test models to our reference model is a tricky issue in this case, because there is not a representation of the activity term in the RVs which we can significantly define the best. Therefore, the choice of model M4 should be considered mainly for illustrative purposes. We found that for model M4 the fitted RV jitter σ_jitt,HARPS-N is nearly half that of model M1, for which we got σ_jitt,HARPS-N = 12.1 ± 1.8 ms⁻¹, suggesting that, at some level, the activity modelling benefits from using the time series of the log ${R^{\prime }}_{\text{HK}}$ index to constrain the GP hyper-parameters. Nonetheless, we note that even the RV uncorrelated jitter is lower for model M5 (${\sigma }_{\text{jitt,HARPS}-\text{N=3}{.8}_{-2.1}^{+1.8}\text{m}\hspace{0.17em}{\text{s}}^{-1}}$), but the masses of planet c are different: for model M4, it is similar to that of Neptune at a ~2.7σ level, while for model M5 we find 7 M_⊕ with ~2σ significance. We show the best-fit solution (model M4) for transits and spectroscopic orbits in Fig. 6.

With more data from new TESS Sectors, we improved the precision of the transit ephemeris and other transit-related parameters, like the planetary radii, with respect to the results of Mann et al. (2020), confirming that the planets have co-planar orbits. Nevertheless, some caution is required about the precision of our determined planetary radii, which, for each planet, is calculated by fitting a light curve obtained by combining all the available transits together, and it is likely optimistic. The unspotted level of HD 63433 is unknown and this introduces a systematic error on the planetary radii. Considering that the flux in the TESS passband is modulated at the level of about 2%, we estimate a systematic error at the level of 1% on the planetary radii from transits close to light maxima and minima, respectively. This comes from the fact that the square of the ratio between the planetary and the stellar radii is proportional to the relative flux deficit at the centre of a transit. Therefore, a change in the out-of-transit flux by about 2% between light maxima and minima will affect the planetary radius derived from the corresponding transits by about 1%. This error estimate is indeed a lower limit because starspots uniformly distributed in longitude do not contribute to the amplitude of the photometric modulation, but affect the measurement of the radius ratio.

An additional systematic effect is produced by the occultations of starspots during transits, but its amplitude is probably lower by about 50% than the effect produced by the changes in the out-of-transit light reference level, if we assume that spots in this G dwarf have a temperature deficit of about 1000 K with respect to the unperturbed photosphere and a filling factor of about 2% in the latitude bands occulted by the planets (e.g. Ballerini et al. 2012). We did not find a significant RV acceleration ($\dot{\gamma }$ = 0.004 ± 0.017 m s⁻¹d⁻¹).

Limits on outer Jovian planets around HD 63433 can be placed by the HIPPARCOS-Gaia proper motion anomaly technique (Brandt 2021; Kervella et al. 2022). For example, the sensitivity limits based on this approach would be compatible with the presence of a ≳1.0 M_Jup companion in the "sweet spot" separation range 3−10 au (Kervella et al. 2022). Interestingly, we note that the reported S/N of the proper motion anomaly for HD 63433 is in the range 2–2.5 (with a growing trend from Gaia DR2 to DR3). While low, it is intriguing, as in this regime of orbital separations Jupiter- and super-Jupiter-mass companions might still be missed by existing RV datasets (see e.g. Fig. 9 of Hirsch et al. 2021), particularly if they were to lie on significantly non-coplanar orbits. Given that the star is nearby, there are good prospects for improved sensitivity to such companions based on Gaia DR4 results, which are expected to be released in late 2025.

We did not measure different nor more precise results for the planetary masses when using the RVs extracted with TERRA. As an example, for the case of models M1, M4, and M6 we got $m_{\text{b}}={3.8}_{-2.6}^{+3.9}$ and $m_{\text{c}}={14.8}_{-7.7}^{+8.0}$, $m_{\text{b}}={2.6}_{-1.8}^{+3.2}$ and $m_{\text{c}}={15.1}_{-7.1}^{+7.5}{M}_{{}_{\oplus }}$, and $m_{\text{b}}={1.8}_{-1.3}^{+2.3}$ and $m_{\text{c}}={5.9}_{-3.8}^{+4.6}{M}_{{}_{\oplus }}$ respectively.

Fig. 4 Time series (left column) and periodograms (right column) of spectroscopic activity diagnostics. The periodogram of the H-alpha index refers to pre-whitened data, after removing the long-term curvature visible in the time series. FAPs are calculated through a bootstrap analysis.

6 Considerations on planetary structure based on the mass-radius diagram

The more conservative conclusion that we can draw from the analysis described in Sect. 5 is that we determined upper limits for the planet masses, m_b ≲ M_⊕ and m_c ;≲ 31 M_⊕. The mass of HD 63433 b remains undetected (and we cannot exclude the possibility that this result is mainly influenced by the choice of a specific model). The mass of HD 63433 c remains (very) uncertain, but the clue that it might be Neptune-like would imply a bulk density that, if confirmed, would be interestingly high for a planet with the size of a mini-Neptune. For this reason, in this Section we comment on this and other possibilities for the planets' physical structures which are compatible with our results, by inspecting the location of the two planets in a mass-radius diagram for exoplanets. Specifically, our considerations concern the solutions which correspond to models M4 and M5 of Table 2).

We show in the upper panel of Fig. 7 a mass-radius diagram which includes well-characterised planets (with masses and radii measured with a precision lower than 30% and 10%, respectively). Assuming the lower value for mc (model M5), HD 63433 c is located in a well populated region of the diagram, corresponding to water worlds, i.e. planets containing significant amounts of H₂O-dominated fluid/ice in addition to rock and gas. Instead, if the estimate of m_c from model M4 is the one more representative of the real mass of HD 63433 c, the planet would occupy a position on the diagram which is less populated, but not empty, suggestive of a structure with a rocky core and a water layer < 50% by mass. The lower panel of Fig. 7 shows a mass-radius diagram which includes only planets with age <900 Myr, with no selection based on the precision of their masses and radii. We note that, assuming the higher mass (density) estimate from model M4, HD 63433 c would have a similar structure of TOI-179b (Desidera et al. 2023; Vines et al. 2023) and Kepler-411 b (Sun et al. 2019). TOI-179 is a K2V star with an age similar to that of the system HD 63433 (400±100 Myr; Desidera et al. 2023). TOI-179 b is a dense mini-Neptune that moves on an eccentric orbit ($e={0.34}_{-0.09}^{+0.07}$) with a shorter orbital semi-major axis (a = 0.0481 ± 0.0004 au) and higher equilibrium temperature (T_eq ~ 990 K, in the case of a circular orbit with same semi-major axis a) than HD 63433 c. Desidera et al. (2023) show that the typical composition of TOI-179 b can be assumed 75% rock + 25% water, and conclude that the planet has likely lost most of its primordial atmosphere, and it is stable against hydrodynamic evaporation. Kepler-411 b is a member of a fairly compact four-planet system orbiting a ~200 Myr old K2V star with a period similar to that of TOI-179 b.

Based on our derived upper limits to the mass of HD 63433 b, we can conclude that its composition could be compatible with that of a water world, and we cannot exclude that it is similar to that of the outermost companion, especially if we consider the results of the M5 model. As for the case of TOI-179 b, noticing that the two planets have similar equilibrium temperatures, HD 63433 b could have lost any primordial gaseous H-He envelope (as we discuss in detail in the next Section), and reached its definitive position on the mass-radius diagram.

Fig. 5 Posterior distributions for the planetary masses for some of the models tested in this work (Table 2). Vertical dashed lines indicate the 95% percentiles.

7 Planetary atmosphere mass-loss due to photoevaporation

We evaluated the mass-loss rate of the planetary atmosphere of planet b and c using the hydro-based approximation developed by Kubyshkina et al. (2018a,b), coupled with the planetary core-envelope model by Lopez & Fortney (2014) and the MESA Stellar Tracks (MIST; Choi et al. 2016). For the stellar X-ray emission at different ages, we adopted the analytic description by Penz et al. (2008), anchored to the current value of the X-ray luminosity, L_★,X = 7.5 × 10²⁸erg s⁻¹ in the band 5–100 Å, derived from XMM-Newton spectra (Zhang et al. 2022). The stellar EUV luminosity (100−920 Å) was computed by means of the new scaling law by Sanz-Forcada et al. (2022). More details on our modelling of atmospheric evaporation are provided in Appendix B. Following Maggio et al. (2022), we expect small but non negligible differences in the evaporation efficiency and timescales if a different X-ray to EUV scaling law is assumed (e.g. King & Wheatley 2021, Johnstone et al. 2021), but a detailed comparison of results is not warranted at this stage, given the uncertainties on the planetary masses.

We performed several simulations of the past and future evolution of the planetary atmospheres, assuming different possible values for the planetary masses at the current age. The grids of test masses are defined by taking into account the upper limits presented in Sect. 5. For planet b, we explored the mass range 2−12 M_⊕, while for planet c we limited our analysis to masses in the range 5−15 M_⊕, because we found that for m_c > 15 M_⊕ the planet is stable against photoevaporation.

First, we considered the cases with the planetary parameters adopted by Zhang et al. (2022) in their 3D hydrodynamic modelling of the hydrogen escape and of the possible Lyα and He I 10833 Å absorption features. In particular, they assumed m_b = 5.5 M_⊕ and R_b = 2.15 R_⊕ for planet b, and m_c = 7.3 M_⊕ and R_c = 2.67 R_⊕ for planet c. Using these values, we predict a current atmospheric mass fraction of ~0.5% for planet b and ~1.4% for planet c, which are in good agreement with the corresponding values of 0.6% and 2% derived by Zhang et al. (2022). We also checked that our X-ray to EUV luminosity scaling provides consistent fluxes at 1 au distance, at the current age: we derived a value of 95 erg s⁻¹ cm⁻², to be compared with the estimate of 91 ± 27 erg s⁻¹ cm⁻² by Zhang et al. (2022), who employed a nominal stellar spectrum from 5 Å to 5 μm. Finally, we compared the predictions of the current mass-loss rates, which resulted in a factor of 1.7 higher for planet b and 2.1 lower for planet c, with respect to the values calculated by Zhang et al. (2022). This relatively small difference is not surprising, given the quite different modelling approaches. However, we stress that these planetary properties are just instantaneous snapshots, which do not take into account the past evolutionary history.

In the following, we report the results of our backward and forward in time simulations, which include the time evolution of the XUV irradiation and of the planetary structure in response to the stellar behaviour. We are interested in investigating how the planetary masses, radii, and atmospheric mass-loss rates change with time due to photoevaporation. The results of the simulations are shown in Figs. 8–9 and in Table 4.

Planet b. First, we computed atmospheric mass fractions and mass-loss rates at the present age for the grid of planet masses, and we simulated the evolution in the future. In Table 4 we report the e-folding mass-loss timescales. These evaporation timescales are very short (1–100 Myr) in any case. Hence, we confirm the conclusion of Zhang et al. (2022), based on observations, that the inner planet has probably lost its atmosphere entirely. Then we explored the possible backward evolutionary pathways, in order to assess the initial planetary masses and radii at the age of 10 Myr (see Appendix B for details on the modelling). In practice, we found that it is always possible to find initial conditions such that the atmosphere is depleted within about 400 Myr. For each assumed planetary mass at the current age, we were able to determine what could be the highest initial mass at 10 Myr. Based on our results for the mass-loss timescales discussed above, we adopt the boundary condition that the atmospheric evaporation is completed at the current age of 414 Myr. These evolutionary histories are shown in Fig. 8. For example, if planet b has m_b = 12 M_⊕ at the present age, it could have started with a mass m_{b, t=10 Myr} ~ 12.02 M_⊕ and an atmospheric fraction f_{atm, t=10 Myr} = 0.14% at 10 Myr.

For planetary masses ≲7 M_⊕ we found a remarkable change in evaporation efficiency. This is due to the relative change of the Jeans escape parameter, Λ, with respect to the control parameter e^Σ in the hydro-based formulation by Kubyshkina et al. (2018b, see Appendix B for details). While for relatively high planetary masses the atmospheric escape is mainly driven by the photo-evaporation, for lower planetary masses the mass-loss rate is increasingly larger because the atmospheric escape is more due to the thermal energy and the lower gravity of the planet. Hence, in order to reach the current age with an atmospheric mass fraction f_{atm,t=414 Myr} = 0, the initial configuration of the planet is characterised by a massive and inflated atmospheric envelope. We found that if the current mass of planet b is 7 M_⊕, m_{b,t=10 Myr} could have been up to about 40 M_⊕, f_{atm,t=10 Myr} = 82%, and the corresponding initial radius was R_{b,t=10 Myr} ~ 14.25 R_⊕. We remark that these are maximum mass values, such that complete evaporation occurs in about 400 Myr. Lower initial masses are also possible, with shorter evaporation timescales and identical atmosphere-depleted status at the present age.

We would like to remind readers that for the backward in time modelling we are considering as a boundary condition that the current mass-loss rate of planet b is zero. That explains the difference with the non-null mass-loss rates indicated in Table 4. In particular, for the cases of a planet with 2 and 4.5 M_⊕, the boundary condition implies the runaway values reached in the late stages before the present age.

A further caveat about the simulations for planets with masses ≲ 7 M_⊕ is that a significant decrease in the planetary mass at early ages implies an outward migration of the planet, with a rate depending on the amount of angular momentum driven away from the system by the planetary wind flow (Fujita et al. 2022). An increase in the star-planet distance yields a decrease in XUV irradiation and equilibrium temperature, hence a decrease in mass-loss rate. Exploring this highly non-linear regime is beyond the scope of the present paper.

Planet c. In this case, the evolution is more simple, with no predicted change of regime (Fig. 9), due to the larger distance of the planet from the host star, a lower equilibrium temperature, and lower high-energy irradiation. In most of the cases that we have examined, i.e. for masses in the range 7.5–15 M_⊕ at the current age, the planet mass and radius are scarcely affected by atmospheric evaporation, although the atmospheric mass fraction can change significantly during the system lifetime. The e-folding mass-loss timescales are always > 5 Gyr, even for the case of a planet with M_p = 7.3 M_⊕ considered by Zhang et al. (2022). The mass-loss timescale of 0.9 Gyr predicted by these authors was significantly lower since it was computed assuming a constant XUV irradiation and mass-loss rate. Instead, in our evolutionary models, the XUV luminosity drops to about 1/3 the current value already at 1.3 Gyr, and it keeps decreasing at later times. Only if a current mass of 5 M_⊕ is considered, the mass and radius change substantially starting from an age of 10 Myr.

Fig. 6 First two rows: TESS light curve and best-fit transit models (left panels), and spectroscopic orbits (HARPS-N DRS RVs; right panels) related to planets HD 63433 b and HD 63433 c. Here, we show the solution obtained with model M4. The upper right plot shows the RV residuals, after removing the Doppler signal due to HD 63433 c and the stellar activity term from the original data, phased to the orbital period of 7.01794 d: it appears clear that the spectroscopic orbit corresponding to HD 63433 b remains not characterised. The phase-folded plot for HD 63433 c shows the RV residuals after removing the stellar activity term from the original data. Third row: Zoomed view of the stellar activity term in the HARPS-N RVs, after removing the Doppler signal due to HD 63433 c from the original RVs, as modelled by a GP quasi-periodic kernel. In model M5, the GP quasi-periodic kernel was trained on the time series of the log $R_{\text{HK}}^{\prime }$ chromospheric activity diagnostic.

Fig. 7 Upper panel: Mass-radius diagram for exoplanets selected from the TEPCAT sample, available at https://www.astro.keele.ac.uk/jkt/tepcat/ (updated to 7 October 2022; Southworth 2011). Black dots represent planets with mass and radius measured with a relative precision lower than 30%and 10%, respectively. Planets of the HD63433 system are indicated by reddish triangles (planet b), which denote mass upper limits, and squares (planet c), with reference to the corresponding models listed in Table 2. Theoretical curves for some planet compositions are overplotted, as calculated by Zeng et al. (2019) assuming 1 milli-bar surface pressure. Models for a planet with an H2O-gaseous atmosphere (50% and 100% H2O by mass, cyan and blue curves, respectively), and for a planet with different percentages of an H2 gaseous envelope over a 50% water-rich layer (magenta curves) are calculated for an isothermal fluid/steam envelope equilibrium temperature of 700 K, similar to that of HD 63433 c (680±12 K). Lower panel: Mass-radius diagram for planets with age < 900 Myr. Triangles are used to identify planets for which only mass upper limits are available. No selection is made based on the precision of mass and radius measurements. We used Desidera et al. (2023) as a reference for the mass and radius of TOI-179 b.

Fig. 8 Evolutionary history of planet HD 63433 b. (a) Mass-loss rate vs. time for different values of the planetary mass at the present age. (b) Fraction of planetary mass lost vs. time. Panels c and d show the evolution of the planetary mass and radius, respectively. Dashed lines indicate backward time evolution from the current age, solid lines forward time evolution. Note that all the simulations were constructed with the boundary condition that the planet completely lost the atmosphere at the present age of 414 Myr.

8 Conclusions

In this paper, we analysed new photometric and spectroscopic data of the 400 Myr-old star HD 63433 (TOI-1726), which hosts a multi-planet system. The presence of two planets in this relatively young system offers the opportunity to examine differences in planetary evolution pathways over the first hundreds of millions years after the system formation, under the influence of the same host star. We provided a characterisation of the system, with special focus on the measurement of the planet radii and masses, and on the investigation of the past and future planetary atmosphere evolution.

Thanks to new TESS photometry and follow-up with HARPS-N, we could improve the planetary ephemeris with respect to the study of Mann et al. (2020), with statistical evidence in favour of circular orbits for the two planets. We revised the measure and improve the precision of the planetary radii ($r_{\text{b}}={2.02}_{-0.05}^{+0.06}{R}_{\oplus }$ and r_c = 2.44 ± 0.07 R_⊕), and provide the first dynamic constraints on the planetary masses. After testing several approaches to mitigate the dominant stellar activity contribution in the RV time series, we can only provide upper limits for the mass m_b of the innermost planet b, conservatively concluding that m_b ≲ 11 M_⊕ at 95% of confidence. As for the larger planet c, not all the models converge to similar mass estimates, but results from a few test models suggest a Neptune-like mass with a significance of 2.1−2.7σ, and that HD 63433 c could be a dense mini-Neptune, with a few known counterparts in the mass-radius diagram. The more conservative conclusion is that m_c ≲ 31 M_⊕ at 95% of confidence.

Based on the constraints to the masses derived from our analysis, we explored the evolution of the atmospheres of both planets, taking into account the decay of stellar activity and XUV irradiation with time. We estimate that the current mass-loss timescale of planet b is sufficiently short to justify the assumption that its atmosphere is already evaporated. This finding is supported by the results of Zhang et al. (2022), who reported no evidence of Lyα or He I absorption due to photoevaporation in dedicated photometric observations with HST.

However, the past history of the planet is quite uncertain, because the evolution of the mass-loss rate is strongly dependent on the assumed mass and current atmospheric mass fraction. For masses ⪞7 M_⊕, the atmospheric mass fraction was already small at early ages (~ 10 Myr), hence the planetary mass and radius changed little in time. Instead, a planet with a lower mass at the present age could have started as an inflated body with a few tens of Earth masses, with a structure dominated by a heavy atmosphere. In any case, the current planetary structure is essentially depleted of a gaseous envelope.

For HD 63433 c, our models indicate that the atmospheric evaporation will keep going over the next ~4.5 Gyr, because the decay of the XUV irradiation determines a steady decrease of the mass-loss rate with time. Our prediction for the mass-loss timescale is different from the shorter value determined by Zhang et al. (2022), who did not take into account the evolution of the stellar activity. Another effect included in our modelling is the time evolution of the stellar bolometric luminosity and planetary equilibrium temperature. Neglecting this effect, the atmosphere of a 13 M_⊕ planet would stop evaporating within 0.6 Gyr, while we found that such a planet keeps losing its atmosphere with a timescale longer than 5 Gyr. Summarising, if the current mass of HD 63433 c is 15 M_⊕, the atmosphere should be stable against evaporation; for masses in the range 7–13 M_⊕ the planet mass and radius do not change appreciably in time, while for the case with m_c = 5 M_⊕ we predict a decrease of 0.8% in mass and 25% in radius over the next ~4.5 Gyr.

We acknowledge that HD 63433 has also been followed-up with the CARMENES spectrograph (Mallorquín Díaz et al. 2023). The GAPS and CARMENES teams have coordinated the submission of two studies, which have been intentionally carried out in an independent way.

Fig. 9 Evolutionary history of planet HD 63433 c. Panel a: mass-loss rate vs. time for different values of the planetary mass at the present age. Panel b: fraction of planetary mass lost vs. time. Panels c and d: time evolution of the planetary mass and radius, respectively. Dashed lines indicate backward time evolution from the current age, solid lines forward time evolution.

Acknowledgements

This work has been supported by the PRIN-INAF 2019 "Planetary systems at young ages (PLATEA)" and ASI-INAF agreement no. 2018-16-HH.0. A.M. and D.L. acknowledge partial financial support from the ASI-INAF agreement no. 2018-16-HH.0 (THE StellaR PAth project), and from the ARIEL ASI-INAF agreement no. 2021-5-HH.0. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Appendix A Description of the GP models and model priors

Hereafter we provide details about the different GP kernels tested in this work to filter out the stellar activity contribution to the observed RV variability.

GP quasi-periodic (QP) kernel. An element of the QP covariance matrix (e.g. Haywood et al. 2014) used in our work is defined as follows:

$$\begin{array}{l}{k}_{QP}\left(t,t^{\prime }\right)=h^2.\exp \left[-\frac{{\left(t-t^{\prime }\right)}^2}{2{\lambda }_{\text{QP}}^2}-\frac{{\sin }^2\left(\pi \left(t-t^{\prime }\right)/\theta \right)}{2{\omega }^2}\right]+\hfill \\ +\left({\sigma }_{\text{RV}}^2\left(t\right)+{\sigma }_{\text{jit}}^2\right)\cdot {\delta }_{t,t^{\prime }}\hfill \end{array}$$(A.1)

Here, t and t' represent two different epochs of observations, σ_RV is the radial velocity uncertainty, and δ_t,t' is the Kronecker delta. Our analysis takes into account other sources of uncorre-lated noise – instrumental and/or astrophysical – by including a constant jitter term σ_jit which is added in quadrature to the formal uncertainties σ_RV. The GP hyper-parameters are h, which denotes the scale amplitude of the correlated signal; θ, which represents the periodic timescale of the correlated signal, and corresponds to the stellar rotation period; w, which describes the "weight" of the rotation period harmonic content within a complete stellar rotation (i.e. a low value of w indicates that the periodic variations contain a significant contribution from the harmonics of the rotation periods); and λ_QP, which represents the decay timescale of the correlations, and is related to the temporal evolution of the magnetically active regions responsible for the correlated signal observed in the RVs.

GP quasi-periodic with cosine (QPC) kernel. The QPC kernel has been introduced by Perger et al. (2021). The covariance matrix element implemented in our work is described by the equation

$$\begin{array}{l}{k}_{QPC}\left(t,t^{\prime }\right)=\exp \left(-2\frac{{\left(t-t^{\prime }\right)}^2}{{\lambda }_{QPC}^2}\right)\cdot [h_1^2\exp \left(-\frac{1}{2{\omega }^2}{\sin }^2\left(\frac{\pi \left(t-t^{\prime }\right)}{\theta }\right)\right)+\hfill \\ +h_2^2\cos \left(\frac{4\pi \left(t-t^{\prime }\right)}{\theta }\right)]+\left({\sigma }_{\text{RV}}^2\left(t\right)+{\sigma }_{\text{jit}}^2\right)\cdot {\delta }_{t,t^{\prime }}\hfill \end{array}$$(A.2)

Again, t and t' represent two different epochs of observations; h₁ and h₂ are scale amplitudes; θ still represents the periodic time-scale of the modelled signal, and corresponds to the stellar rotation period; w still describes the weight of the rotation period harmonic content within a complete stellar rotation; λ_QPC is defined as 2 λ_QP, better representing the average lifetime of the activity-related features responsible for the stellar correlated signal in the RVs according to Perger et al. (2021); σ_RV(t) and σ_jit are the radial velocity uncertainty at epoch t and the uncorrelated jitter, respectively, and δ_t,t' is the Kronecker delta.

GP rotational or double simple harmonic oscillator (dSHO) kernel. The "rotational" kernel is defined by a mixture of two stochastically driven, damped simple harmonic oscillators (SHOs) with undamped periods of P_★,rot and P_★,rot/2. This can be obtained by combining two SHOTerm kernels included in the package celerite (Foreman-Mackey et al. 2017)³. The power spectral density corresponding to this kernel is

$$\begin{array}{l}S\left(\omega \right)=\sqrt{\frac{2}{\pi }}\frac{S_1{\omega }_1^4}{{\left({\omega }^2-{\omega }_1^2\right)}^2+2{\omega }_1^2{\omega }^2}+\hfill \\ +\sqrt{\frac{2}{\pi }}\frac{S_2{\omega }_2^4}{{\left({\omega }^2+{\omega }_1^2\right)}^2+2{\omega }_2^2/Q^2},\hfill \end{array}$$(A.3)

where

$$S_1=\frac{A^2}{{\omega }_1{Q}_1\left(1+f\right)},$$(A.4)

$$S_2=\frac{A^2}{{\omega }_2{Q}_2\left(1+f\right)}\cdot f,$$(A.5)

$${\omega }_1=\frac{4\pi Q_1}{P_{\text{rot}}\sqrt{4Q_1^2-1}},$$(A.6)

$${\omega }_2=\frac{8\pi Q_2}{P_{\text{rot}}\sqrt{4Q_2^2-1}},$$(A.7)

$$Q_1=\frac{1}{2}+Q_0+\text{Δ}Q,$$(A.8)

$$Q_2=\frac{1}{2}+Q_0.$$(A.9)

The parameters in [A.4-A.9], where the subscripts 1 and 2 refer to the primary (P_★,rot) and secondary P_★,rot/2) modes, represent the inputs to the SHOTerm kernels. However, instead of using them directly, we adopt a different parametrization using the following variables as free hyper-parameters in the MC analysis, from which S_i, Q_i, and ω_i are derived through Eq. [A.4-A.9]: the variability amplitude A, the stellar rotation period P_★,rot, the quality factor Q₀, the difference ΔQ between the quality factors of the first and second terms, and the fractional amplitude f of the secondary mode relative to the primary.

Multi-dimensional GP framework. This framework has been first introduced by Rajpaul et al. (2015). A python version of the code is implemented in the package pyaneti (Barragán et al. 2019a)⁴. In our work, we included the module pyaneti in our own code to make it work with the nested sampler MultiNest. In brief, the multi-dimensional GP framework uses ancillary data sensitive to stellar activity in a combined fit with the RVs (in our case, we used the activity diagnostics BIS and log $R_{\text{HK}}^{\prime }$), with each dataset represented in the time domain by linear combinations of the same GP-drawn function G(t) and its derivative $\dot{G}\left(t\right)$. For all the details of the multi-dimensional GP implementation see Barragán et al. (2022a), especially Sect. 3.1.2 and Eq. (8) therein. For our test, we adopted a QP kernel, as described above, and we fixed to zero the coefficient of $\dot{G}\left(t\right)$ for the case of log $R_{\text{HK}}^{\prime }$.

The prior distributions for the free parameters of the models tested in this work are provided in Table A.1.

Appendix B Details on atmospheric evaporation modelling

The hydro-based approximation developed by Kubyshkina et al. (2018a,b) provides an analytic expression of the planetary mass-loss rate as a function of planetary mass and radius, star-planet distance, and stellar high-energy flux. We employed this approximation to model the time evolution of the planet characteristics, taking into account both the stellar evolutionary track and the long-term change of the XUV irradiation.

For the X-ray evolution as a function of time, we used the description by Penz et al. (2008), based on open cluster and field G-type stars. It describes the saturation and decay phases of stellar activity with simple analytic power laws. Maggio et al. (2022) showed that they represent a good approximation of the more complex behaviour modelled by Johnstone et al. (2021) for late-type stars. In order to evaluate the stellar irradiation in an extreme ultraviolet band we adopted a scaling law between EUV (100-920 Å) and X-ray (5–100 Å) luminosities, calibrated on stars observed in both bands, derived by Sanz-Forcada et al. (2022), namely

$$\log L_{\text{EUV}}=\left(0.793\pm 0.058\right)\log L_{\text{x}}+\left(6.53\pm 1.61\right)$$(B.1)

which is an updated version of the better-known and widely used relationship by Sanz-Forcada et al. (2011). The predicted evolution of the X-ray and EUV stellar luminosities is shown in Fig. C.4.

We also took into account the evolution of the planetary radius according to Lopez & Fortney (2014). This relation was developed for H-He dominated atmospheres, and provides the envelope radius R_env, as a function of the planetary mass, the atmospheric mass fraction f_atm, the bolometric flux received, and the age of the system. It takes into account also the cooling and contraction of the envelope as a consequence of its thermal evolution, based on the grid of models presented in Lopez et al. (2012).

The track followed by HD 63433 in the theoretical temperature-luminosity diagram is shown in Fig. C.3. It was computed for a star with a mass 0.994M_⊙ and metallicity [Fe/H] = +0.02 (Sect. 3) using the web-based interpolator⁵ of the MESA Isochrones and Stellar Tracks (MIST, Choi et al. 2016).

The Lopez & Fortney (2014) relation can be inverted to estimate f_atm at the current age, starting with the known (measured) planetary radius and with a core radius, R_core, computed as in Benatti et al. (2021), assuming that all the planetary mass is concentrated in the core. Different assumptions (or future measurements) of the planetary mass, yield different R_core and f_atm values.

We then proceeded with the time evolution. For each time step of the simulation, we computed the mass-loss rate and we updated f_atm and the planetary mass, thus obtaining a new value of R_env with the relation by Lopez & Fortney (2014). The latter quantity added to the core radius (assumed constant) provides the updated planetary radius. For each starting planetary mass, we let the system evolve from the current age (414 Myr for HD 63433) until 5 Gyr. In order to check if the atmosphere is hydro-dynamically stable or not, we use the Jeans escape parameter

$$\text{Λ}=\frac{G{m}_{\text{H}}{M}_{\text{p}}}{k_{\text{B}}{T}_{\text{eq}}{R}_p}$$(B.2)

where G is the gravitational constant, mH is the hydrogen mass, and k_B is the Boltzman constant, while M_p and R_p are the planet mass and radius, and T_eq is the planet equilibrium temperature, computed as

$$T_{\text{eq}}=T_{\text{eff}}{\left[f_{\text{p}}\left(1-A_{\text{B}}\right)\right]}^{1/4}{\left(\frac{R_{\ast }}{2d}\right)}^{1/2}$$(B.3)

where T_eff and R_* are the stellar effective temperature and radius, respectively, A_B is the Bond albedo, d the star-planet distance, f_p is a parameter that takes into account if the planet is tidally locked or not. We assumed A_B = 0.5, and $f_{\text{p}}=\frac{2}{3}$, because tidal locking is assumed as in (Zhang et al. 2022).

Note that Λ depends on the stellar T_eff through T_eq, and Λ enters also in the analytic expression for the atmospheric mass-loss rate derived by Kubyshkina et al. (2018b), which is also a function of the XUV flux, F_XUV, at the planetary orbital position. According to Fossati et al. (2017), photoevaporation of planetary atmospheres occurs when the Jeans escape parameter is smaller than a critical value (Λ < 80), but two different evaporation regimes are possible, depending on how Λ compares with a control parameter e^Σ, determined by d and by time-dependent values of R_p and F_XUV (Kubyshkina et al. 2018b): if Λ > e^Σ, the evaporation is controlled by the XUV flux, while for smaller values of Λ the atmospheric escape is mainly driven by the planetary intrinsic thermal energy and low gravity.

In this work, assuming an X-ray luminosity of 5 × 10²⁸ erg/s at the current age of the system, and allowing for the evolution of the bolometric luminosity and the effective temperature of the star as described above (Fig. C.3), we investigated the hydrodynamic atmospheric evolution of the two planets of HD 63433 in the future. Furthermore, we studied the system's evolution back in time, as described by Georgieva et al. (2021). According to the aforementioned scenario and assuming a core which does not change in size or mass, we created a synthetic population of planets with different initial masses, and we assigned to them the core radius and core mass fixed by the total mass assumed for the planet at its present age. This framework leaves f_atm to dictate the total initial mass, while the total radius is again based on the analytic fit by Lopez & Fortney (2014). Then we searched for the planet in the synthetic set which ends up with mass, radius and f_atm most similar to the values at the present age, and we selected its predicted evolutionary history. We started all the simulations in the past at an age of 10 Myr, when we suppose that the circumstellar disk is already disappeared and all planets are in their final, stable orbit. These simulations always end at an age of 414 Myr.

Appendix C Additional plots

Fig. C.1 Time series and GLS periodograms of the light curves collected by STELLA in V and I bands. The horizontal dashed line denotes the FAP level of 10%, the vertical dashed line identifies the stellar rotation period.

Fig. C.2 RVs as a function of the CCF-BIS, showing the significant correlation between the two dataset.

Fig. C.3 Evolutionary track of HD 63433 up to an age of 5 Gyr in the effective temperature-bolometric luminosity plane. The current location of the star on the track is marked by a black cross.

Fig. C.4 Time evolution of X-ray (5–100 Å), EUV (100–920 Å), and total XUV luminosity of HD 63433, according to Penz et al. (2008) and the X-ray/EUV scaling by Sanz-Forcada et al. (2022). The grey area is the observed 1σ spread around the median (dark grey line) of the X-ray luminosity distributions for dG stars in the Hyades and Pleiades open clusters. Uncertainties on the age and X-ray luminosity of HD 63433 are within the size of the square symbol.

Appendix D Dataset

References

All Tables

All Figures

	Fig. 1 TESS light curve of HD 63433 for Sector 20 (upper panel), and Sectors from 44 to 47 (lower panel). Planetary transits and bad-quality points have been removed from the dataset. The red curve represents the best-fit GP model.
In the text

	Fig. 2 Spectral energy distribution of the host star HD 63433 with the best-fit model overplotted (solid line). Red and blue points correspond to the observed and predicted values, respectively.
In the text

Fig. 3 Time series of the RVs calculated from HARPS-N spectra with the DRS and TERRA pipelines, and their GLS periodograms. First column: RV timeseries. The mean value has been subtracted from the original DRS data. Second column: The corresponding GLS periodograms of the two RV dataset. FAP levels are of 1% and 10% are indicated by yellow and red horizontal lines, respectively. They are calculated through a bootstrap analysis. Third column: GLS periodograms of the pre-whitened RV data.

In the text

	Fig. 4 Time series (left column) and periodograms (right column) of spectroscopic activity diagnostics. The periodogram of the H-alpha index refers to pre-whitened data, after removing the long-term curvature visible in the time series. FAPs are calculated through a bootstrap analysis.
In the text

	Fig. 5 Posterior distributions for the planetary masses for some of the models tested in this work (Table 2). Vertical dashed lines indicate the 95% percentiles.
In the text

Fig. 6 First two rows: TESS light curve and best-fit transit models (left panels), and spectroscopic orbits (HARPS-N DRS RVs; right panels) related to planets HD 63433 b and HD 63433 c. Here, we show the solution obtained with model M4. The upper right plot shows the RV residuals, after removing the Doppler signal due to HD 63433 c and the stellar activity term from the original data, phased to the orbital period of 7.01794 d: it appears clear that the spectroscopic orbit corresponding to HD 63433 b remains not characterised. The phase-folded plot for HD 63433 c shows the RV residuals after removing the stellar activity term from the original data. Third row: Zoomed view of the stellar activity term in the HARPS-N RVs, after removing the Doppler signal due to HD 63433 c from the original RVs, as modelled by a GP quasi-periodic kernel. In model M5, the GP quasi-periodic kernel was trained on the time series of the log $R_{\text{HK}}^{\prime }$ chromospheric activity diagnostic.

In the text

Fig. 7 Upper panel: Mass-radius diagram for exoplanets selected from the TEPCAT sample, available at https://www.astro.keele.ac.uk/jkt/tepcat/ (updated to 7 October 2022; Southworth 2011). Black dots represent planets with mass and radius measured with a relative precision lower than 30%and 10%, respectively. Planets of the HD63433 system are indicated by reddish triangles (planet b), which denote mass upper limits, and squares (planet c), with reference to the corresponding models listed in Table 2. Theoretical curves for some planet compositions are overplotted, as calculated by Zeng et al. (2019) assuming 1 milli-bar surface pressure. Models for a planet with an H2O-gaseous atmosphere (50% and 100% H2O by mass, cyan and blue curves, respectively), and for a planet with different percentages of an H2 gaseous envelope over a 50% water-rich layer (magenta curves) are calculated for an isothermal fluid/steam envelope equilibrium temperature of 700 K, similar to that of HD 63433 c (680±12 K). Lower panel: Mass-radius diagram for planets with age < 900 Myr. Triangles are used to identify planets for which only mass upper limits are available. No selection is made based on the precision of mass and radius measurements. We used Desidera et al. (2023) as a reference for the mass and radius of TOI-179 b.

In the text

Fig. 8 Evolutionary history of planet HD 63433 b. (a) Mass-loss rate vs. time for different values of the planetary mass at the present age. (b) Fraction of planetary mass lost vs. time. Panels c and d show the evolution of the planetary mass and radius, respectively. Dashed lines indicate backward time evolution from the current age, solid lines forward time evolution. Note that all the simulations were constructed with the boundary condition that the planet completely lost the atmosphere at the present age of 414 Myr.

In the text

	Fig. 9 Evolutionary history of planet HD 63433 c. Panel a: mass-loss rate vs. time for different values of the planetary mass at the present age. Panel b: fraction of planetary mass lost vs. time. Panels c and d: time evolution of the planetary mass and radius, respectively. Dashed lines indicate backward time evolution from the current age, solid lines forward time evolution.
In the text

	Fig. C.1 Time series and GLS periodograms of the light curves collected by STELLA in V and I bands. The horizontal dashed line denotes the FAP level of 10%, the vertical dashed line identifies the stellar rotation period.
In the text

	Fig. C.2 RVs as a function of the CCF-BIS, showing the significant correlation between the two dataset.
In the text

	Fig. C.3 Evolutionary track of HD 63433 up to an age of 5 Gyr in the effective temperature-bolometric luminosity plane. The current location of the star on the track is marked by a black cross.
In the text

Fig. C.4 Time evolution of X-ray (5–100 Å), EUV (100–920 Å), and total XUV luminosity of HD 63433, according to Penz et al. (2008) and the X-ray/EUV scaling by Sanz-Forcada et al. (2022). The grey area is the observed 1σ spread around the median (dark grey line) of the X-ray luminosity distributions for dG stars in the Hyades and Pleiades open clusters. Uncertainties on the age and X-ray luminosity of HD 63433 are within the size of the square symbol.

In the text