The GAPS Programme at TNG

Context. The GAPS collaboration is carrying out a spectroscopic and photometric follow-up of a sample of young stars with planets (age ≲ 600Myr) to characterise planetary systems at the early stages of their evolution. Aims. For more than 2yr, we monitored with the HARPS-N spectrograph the 400Myr-old star HD63433, which hosts two close-in (orbital periods P b ∼ 7 . 1 and P c ∼ 20 . 5 days) sub-Neptunes detected by the TESS space telescope, and it was announced in 2020. Using radial velocities and additional TESS photometry, we aim to provide the first measurement of their masses, improve the measure of their size and orbital parameters, and study the evolution of the atmospheric mass-loss rate due to photoevaporation. Methods. We tested state-of-the-art analysis techniques and different models to mitigate the dominant signals due to stellar activity that are detected in the radial velocity time series. We used a hydro-based analytical description of the atmospheric mass-loss rate, coupled with a core-envelope model and stellar evolutionary tracks, to study the past and future evolution of the planetary masses and radii. Results. We derived new measurements of the planetary orbital periods and radii ( P b = 7 . 10794 ± 0 . 000009 days, r b = 2 . 02 + 0 . 06 − 0 . 05 R ⊕ ; P c = 20 . 54379 ± 0 . 00002 days, r c = 2 . 44 ± 0 . 07 R ⊕ ), and determined mass upper limits ( m b < ∼ 11 M ⊕ ; m c < ∼ 31 M ⊕ ; 95% confidence level), with evidence at a 2.1–2.7 σ significance level that HD63433c might be a dense mini-Neptune with a Neptune-like mass. For a grid of test masses below our derived dynamical upper limits, we found that HD63433b has very likely lost any gaseous H-He envelope, supporting HST-based observations that are indicative of there being no ongoing atmospheric will keep evaporating over the next ∼ 5 Gyr if its current mass is m c < ∼ 15 M ⊕ , while it should be hydrodynamically stable for higher masses.


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' distribu-tion 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 m s −1 ) 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. 2018Carleo et al. , 2020;;Damasso et al. 2020), the main focus of the GAPS young planets subprogramme 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. 2022).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 twoplanet 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.
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. 2016Malavolta et al. , 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).

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 arcseconds 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.

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)1 .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 for a 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 −1 , and an rms of 23.3 m s −1 .
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 −1 , and an rms of 24.1 m s −1 .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 HK from the CaII 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.

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 g (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 B p − R p ), 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.54 dex (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.10km/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

Parameter
Value Ref.
References. 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 et al. 2016, 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 d.This result is in agreement with the measurement by Mann et al. ( 2020

Frequency content analysis of radial velocities and activity diagnostics
We 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.

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 Z calculated by MultiNest, and using the empirical scale indicated in Table 1 of Feroz et al. (2011) to assess their relative statistical significance2 .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 Eq. 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 two years.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 θ = 6.390 +0.007 −0.005 days for M1), even with an adopted uniform and large prior U(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 m c of HD 63433 c with a statistical significance (defined as the ratio m c /σ − m 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 (quasiperiodic 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 m b and m c 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. (2022) found that TOI-179 b (further discussed in Sect.6) has a mass of 24.1 +7.1 −7.7 M ⊕ (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 (∼ 24 m/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. 2019), is a 750 Myrold planet with an orbital period of 1.7-d, and a measured mass of 21.8±6.2M ⊕ (3.5σ precision).Similarly, mass measurements have been claimed for TOI-560 b (480 Myr; m = 10.2 +3.4  −3.1 M ⊕ , 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.8m s −1 , suggesting that, at some level, the activity modelling benefits from using the time series of the log R HK index to constrain the GP hyper-parameters.Nonetheless, we note that even the RV uncorrelated jitter is lower for model M5 (σ jitt, HARPS−N = 3.8 +1.8  −2.1 m 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 (γ = 0.004 ± 0.017 m s −1 d −1 ).
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 b = 3.8 +3.9 M ⊕ , and m b = 1.8 +2.3 −1.3 and m c = 5.9 +4.6 −3.8 M ⊕ respectively.

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 < ∼ 11 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 m c (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 2 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 Article number, page 7 of 22 Table 2. Summary of the different RV (DRS) + transit photometry + activity diagnostics (when specified) models tested in this work.All include two Keplerians for planets b and c, except for models M5 and M6, for which we adopted circular orbits.Multi-dimensional GP framework d 836.6 ± 0.3 m b = 1.2 +1.4 −0.9 (3.9) ρ b = 0.8 +1.0 −0.6 (2.6) using log R HK and BIS, and the quasi-periodic kernel; circular orbits m c = 7.0 +3.2 −3.0 (12.1) Notes. (a) We note that the analysed datasets related to each model are not always the same, therefore model comparison based on the values of the Bayesian evidence ln Z can be performed only when the same datasets are involved. (b) The 95 th percentile is provided in parenthesis. (c) With "trained" we mean that the same GP kernel is used to model the RVs and the activity diagnostic.In the case of a quasi-periodic kernel, three hyperparameters over four, namely w,θ, and λ, are shared by the two datasets, while two distinct amplitudes h are used. (d) To speed up the analysis, we did not fit the photometric transits, and we assumed circular orbits (the results from other models show that the eccentricities are consistent with zero).The priors for the planet ephemeris are defined based on the results of the M4 model.We adopted planet radii calculated from M4 to derive the densities.
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-179 b (Desidera et al. 2022;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. 2022).TOI-179 b is a dense mini-Neptune that moves on an eccentric orbit (e = 0.34 +0.07 −0.09 ) 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 semimajor axis a) than HD 63433 c.Desidera et al. (2022) 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 fourplanet 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.

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 Xray luminosity, L , X = 7.5 × 10 28 erg s −1 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 Article number, page 8 of 22 18.9 +7.0 −6.9 (30.3) mean density, ρ c [g cm −3 ] 7.1 +2.8 −2.6 (11.7)T eq., c [K] 680 Notes. (a) The uncertainties are given as the 16 th and 84 th percentiles of the posterior distributions.For some of the parameters, we provide the 95 th percentile in parenthesis. (b) In place of a/R , we used the stellar density ρ * as a free parameter with Gaussian prior N(1.376,0.078)ρ , from which we derived the a p /R ratios at each step of the MC sampling. (c) A systematic error of at least 1% of the best-fit radius should be added, to take into account effects due to stellar activity (Sect.5).
Article number, page 9 of 22  (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 −1 cm −2 , to be compared with the estimate of 91 ± 27 erg s −1 cm −2 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 Figures 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=10Myr ∼ 12.02 M ⊕ and an atmospheric fraction f atm,t=10Myr = 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 Article number, page 10 of 22 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 photoevaporation, 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=414Myr = 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=10Myr could have been up to about 40 M ⊕ , f atm,t=10Myr = 82%, and the corresponding initial radius was R b, t=10Myr ∼ 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 efolding mass-loss timescales are always > 5 Gyr, even for the case of a planet with M p = 7.3M ⊕ 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.
Article number, page 11 of 22 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 HD 63433 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 H 2 O-gaseous atmosphere (50% and 100% H 2 O by mass, cyan and blue curves, respectively), and for a planet with different percentages of an H 2 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. (2022) as a reference for the mass and radius of TOI-179 b.

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 re-vised the measure and improve the precision of the planetary radii (r b = 2.02 +0.06 −0.05 R ⊕ and r c = 2.44 ± 0.07R ⊕ ), 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.
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.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 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 hydrodynamically stable or not, we use the Jeans escape parameter where G is the gravitational constant, m H 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

Fig. 1 .
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.
), and it is confirmed by looking at the periodogram 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 d.The fundamental stellar properties of HD 63433 are summarised in Tab. 1.

Fig. 2 .
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.

Fig. 3 .
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.

Fig. 4 .
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.

Fig. 5 .
Fig. 5. Posterior distributions for the planetary masses for some of the models tested in this work (Table2).Vertical dashed lines indicate the 95% percentiles.
Notes.(a)Mass estimated byZhang et al. 2022.They assumed r b = 2.15 R ⊕ and r c = 2.67 R ⊕ as the current values of the planetary radii.

Fig. 6 .
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 HK chromospheric activity diagnostic.
(a) U(a, b) denotes a prior drawn from an uninformative distribution in the range (a,b).N(a, b) denotes a prior drawn from a Gaussian distribution with mean a and sigma b.
(2022), namely log L EUV = (0.793 ± 0.058) log L x + (6.53 ± 1.61) (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.
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 5 of the MESA Isochrones and Stellar Tracks (MIST, Choi et al. 2016).

Table 3 .
Best-fit values of the free parameters of model M4 (RVs calculated by the DRS).

Table 4 .
Results of the planetary mass-loss simulations from the current age to 5 Gyr (Sect.7).We assumed r b = 2.02 and r c = 2.44 R ⊕ as the current values of the planetary radii.

Table A .
1. Priors used for the models listed in Table2.