Orbital obliquity of the young planet TOI-5398 b
and the evolutionary history of the system ^††thanks: Table A.1 is only available in electronic form at the CDS via anonymous ftp to https://cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/.^††thanks: Based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated by the Fundación Galileo Galilei (FGG) of the Istituto Nazionale di Astrofisica (INAF) at the Observatorio del Roque de los Muchachos (La Palma, Canary Islands, Spain).

The RV spectroscopic time series will be available in electronic format as supplementary material for the Letter at the CDS.

To properly take into account the errors associated with the orbital parameters, we used Gaussian priors on the orbital period ($P_{\rm b}$), the central time of the first transit ($T_{0,{\rm b}}$), the planet-to-star radius ratio ($R_{\rm b}/R_{\star}$), and the impact parameter ($b$). We imposed a Gaussian prior on the eccentricity following Van Eylen et al. (2019) and included the RV semi-amplitude $K_{\rm b}$ in the modelling. We also imposed Gaussian priors on the host star density ($\rho_{\star}$), rotational period ($P_{\rm rot}$), and projected rotational velocity ($v\sin{i_{\star}}$). In particular, we treated the stellar equatorial velocity ($v_{\rm eq}$, obtained from $P_{\rm rot}$) and inclination ($i_{\star}$) as independent variables. Therefore, we are not subject to the bias described in Masuda & Winn (2020). We treated the stellar limb-darkening (LD) contribution by estimating $\rm u_{1}$ and $\rm u_{2}$ using PyLDTk⁵⁵5https://github.com/hpparvi/ldtk (Husser et al. 2013; Parviainen & Aigrain 2015) and assuming a boxcar filter as the passband in the HARPS-N spectral range. We also modelled and fit the effect of the stellar convective blueshift (CB) on the in-transit RV curve (e.g. Cegla et al. 2016) using a linear law as a function of the limb angle. Short-term stellar activity was included in the model through a jitter term to be added in quadrature to the errors of the RV.

We performed a global optimisation of the parameters by running PyDE (Storn & Price 1997) and performing a Bayesian analysis of the RM signal in the RV time series using emcee Foreman-Mackey et al. (2013). This choice was motivated by the fact that MCMC samplers allow the specification of a prior on derived parameters, as in the case of $v\sin{i_{\star}}$, which is not possible in the case of nested sampling algorithms.

We performed a preliminary analysis including a Gaussian Process with a Matern-3/2 kernel (e.g. Palle et al. 2020) to model the stellar activity, and the simultaneous modelling of the RM signals associated with both planets. The results are perfectly consistent with those presented in Table 1 with comparable error bars, while the obliquity of planet c remains unconstrained. For this reason, we preferred the simpler model presented in Sect. 3. A more complex modelling will be presented after securing more spectroscopic and photometric in-transit data.

We considered a template protoplanetary disc similar to the nominal one from Turrini et al. (2023) (see Table 2 for a summary of its main physical parameters) and sampled 4$\times$10⁴ planet formation tracks. Half of our samples assume centimetre(cm)-sized pebbles and the other half millimetre(mm)-sized pebbles, to explore the impact of faster and slower growth rates of the planetary cores (Johansen et al. 2019). The resulting populations of synthetic planets are shown in Fig. 5. Each formation track is characterised by a different combination of the initial time and orbital distance of the planetary seed, both values being obtained using uniform extractions (see Table 2 for the relevant ranges and see Fig. 4 for illustrative examples of the growth tracks).

Because of the uncertainties of the real characteristics and lifetime of TOI-5398’s now-dispersed disc, we considered as successful all planets between 8 and 80 $M_{\oplus}$, that is roughly between 20% lower than the mass of planet c and 20% higher than the mass of planet b at their current age, and whose final orbits fall within 0.2 AU. This region in the space of the solutions is identified by the black box in Fig. 5. Consistently with Johansen et al. (2019), we find that forming planets in this mass range is more challenging than forming less massive super-Earths ($<5$ $M_{\oplus}$) and more massive gas giants ($>100$ $M_{\oplus}$), and that it is comparatively easier to form synthetic counterparts of planet b than of planet c (see Fig. 5).

To determine the atmospheric mass fraction at present age, we adopted the core-envelope structure by Fortney et al. (2007) and Lopez & Fortney (2014). This core-envelope model has four unknowns: the mass and radius of the core, $M_{\rm core}$ and $R_{\rm core}$, the radius of the envelope, $R_{\rm env}$, and the atmospheric mass fraction, $f_{\rm atm}$. These quantities are linked by four relations and two observational constraints: the measured planet radius and total mass, $R_{\rm p}$ and $M_{\rm p}$; the Lopez & Fortney (2014) relation which links $R_{\rm env}$, $f_{\rm atm}$, and $M_{\rm p}$; and a relation between $R_{\rm core}$, $R_{\rm p}$, and $M_{\rm p}$ at the given age and star–planet distance. For the latter, we adopted the internal structure models by Fortney et al. (2007), where cores can be composed of different ice-rock or rock-iron mixtures.

We note that the relation by Lopez & Fortney (2014) was developed for H–He-dominated atmospheres, and it also takes into account the cooling and contraction of the envelope as a consequence of its thermal evolution (Lopez et al. 2012). For atmospheres with solar metallicity, the radius of the envelope decreases with age as $t^{-0.11}$, while for enhanced opacities the variation is $\propto t^{-0.18}$.

Assuming cores composed of 50% rocks and 50% ices, and a high-metallicity envelope, we found for planet b a first closed solution with $M_{\rm core,b}=16.9$ $M_{\oplus}$ and $R_{\rm core,b}=2.8$ $R_{\oplus}$, corresponding to the measured planetary mass and optical radius (Fig.6). The alternative assumption of solar metallicity yields a maximum value of the planetary radius that is about 7% lower than the available measurement. The core and envelope densities result around $4.2$ g cm^-3 and $0.55$ g cm^-3, respectively. The mean planetary density is around $0.3$ g cm^-3.

For planet c, we obtained $M_{\rm core,c}=11.7$ $M_{\oplus}$ and $R_{\rm core,c}=2.6$ $R_{\oplus}$, a solution self-consistent with the assumed core-envelope structure and the measured mass and optical radius. These values imply a core density of $\sim 3.6$ g cm^-3, an envelope density $\sim 1.3$ g cm^-3, and a mean planetary density $\sim 1.5$ g cm^-3. These core sizes and masses are larger than predicted by the planetary formation models described in Sect. 4.3, while the densities are lower.

Alternatively, we searched for a second possible configuration assuming a denser Earth-like core composition with 67% rocks and 33% iron. In this case, we found a solution only for planet b with a smaller and less massive core (Fig.6), having $M_{\rm core,b}=2.0$ $M_{\oplus}$ and $R_{\rm core,b}=1.2$ $R_{\oplus}$, more in line with the prediction of the planetary formation models.

We adopted the ATES photoionisation hydrodynamics code by Caldiroli et al. (2021, 2022). To speed up the computation, we resorted to their analytical approximation, which provides the planetary mass-loss rate as a function of planetary mass and radius, star–planet distance, and stellar high-energy flux. This is achieved by evaluating the effective photoevaporation efficiency, $\eta_{\rm eff}$, which enters into the classical energy-limited formulation (Erkaev et al. 2007):

where $F_{\rm XUV}$ is the X-EUV flux at the (average) orbital distance, $\rho_{\rm p}$ is the mean planetary mass density, and the factor K accounts for the host star tidal forces (Erkaev et al. 2007). In all cases, we consider the stellar evolutionary track and the long-term change of the XUV irradiation.

The track followed by TOI-5398 in the theoretical temperature–luminosity diagram is shown in Fig. 7, and was computed using the web-based interpolator⁷⁷7https://waps.cfa.harvard.edu/MIST/interp_tracks.html of the MESA Isochrones and 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_{\star,\,X}=10^{29}$ erg s^-1 in the band 5–100 Å, derived from the rotation–activity relationships by Pizzolato et al. (2003). To evaluate the stellar irradiation in the XUV band, we adopted the scaling law between EUV (100–920 Å) and X-ray (5–100 Å) luminosities derived by Sanz-Forcada et al. (2022):

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

To determine the evolution of the planetary radius associated with the variation of the planetary mass, we adopted the core–envelope structure model described in Appendix D. With this model, we incorporate a time dependence of the envelope size linked to the variation of the surface equilibrium temperature during the evolution.

We then proceeded with the time evolution of the photoevaporation. For each time step of the simulation, we compute the mass-loss rate and update $f_{\rm atm}$ and the planetary mass, obtaining a new value of $R_{\rm env}$ with the relation by Lopez & Fortney (2014). The latter quantity added to the core radius (assumed constant) provides the updated planetary radius.

According to the aforementioned scenario and assuming a core that does not change in size or mass, we followed the planetary evolution back in time. We stopped our simulations in the past at 3 Myr, the end time of planetary formation models, that is, when the circumstellar disc has already disappeared and each planet is in its final, stable orbit. The results of the simulations, described in Sect. 4.4, are shown in Figures 8 and 9, and in Table 3.

The mass-loss rate for planet b closely follows the broken power-law decrease in the XUV irradiation over time (Fig. 7). The variation for planet c appears smoother and faster. The reason for the different behaviour is that the mean density of planet b changes little during the evolution, while it increases substantially for planet c. This also explains why the mass-loss rate of planet c becomes lower than that of planet b at an age of around $530$ Myr.

We checked the deviations of the actual numerical results of the ATES code from the analytic approximation we employed. Our computed mass-loss rates are systematically lower by a factor of 1.5–2 than the more precise values. This difference is negligible for the evolution of planet b but may imply an underestimation of the initial mass and radius for planet c (by less than a factor of 2). The mass-loss rate could be overestimated by a similar factor if the EUV flux were predicted with the scaling law proposed by Johnstone et al. (2021) instead of the adopted Eq. 3.

We stress that for the backward-in-time modelling, we are considering as initial conditions the current planetary core–envelope structure (Appendix D). These initial conditions provide strong constraints on the possible evolutionary histories of the two planets.

The system’s dynamics is currently dominated by the mutual secular perturbations of planets and tides. In the long term, the system is stable even if the outer planet (b) were closer than the nominal orbit (up to a limit of $a_{b}>0.07$ AU) and a high eccentricity (0.2) is assumed. At present, the two planets are very close to a 9:4 mean motion resonance, which may have been broken after the dissipation of the gaseous disc. It is a weak resonance, but it may be responsible for an eccentricity value slightly different from zero for both planets, built during the resonance capture induced by the migration of the planets while embedded in the disc.

If planet b were in an almost circular orbit, the secular perturbations on planet c would cause only a low forced eccentricity (of the order of 0.001), which might trigger some tidal deformation — as in the case of Io and Jupiter. On the other hand, if planet b were on an eccentric orbit with a value equal to the upper boundary (0.12) given by the data analysis (see Table 1) then the forced eccentricity of planet c would be about 0.08 with a sinusoidal oscillation amplitude of 0.16. In the latter case, planet c would be likely to undergo strong tidal deformations, which would possibly represent an additional heat source.