Issue 
A&A
Volume 656, December 2021



Article Number  A157  
Number of page(s)  36  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202142010  
Published online  16 December 2021 
Constraining stellar rotation and planetary atmospheric evolution of a dozen systems hosting subNeptunes and superEarths
^{1}
Space Research Institute, Austrian Academy of Sciences, Schmiedlstrasse 6, 8042 Graz, Austria
email: andrea.bonfanti@oeaw.ac.at
^{2}
School of Physics, Trinity College Dublin, the University of Dublin, College Green, Dublin2, Ireland
Received:
13
August
2021
Accepted:
11
October
2021
Context. Planetary atmospheric evolution modelling is a prime tool for understanding the observed exoplanet population and constraining formation and migration mechanisms, but it can also be used to study the evolution of the activity level of planet hosts.
Aims. We constrain the planetary atmospheric mass fraction at the time of the dispersal of the protoplanetary disk and the evolution of the stellar rotation rate for a dozen multiplanet systems that host subNeptunes and/or superEarths.
Methods. We employ a customdeveloped PYTHON code that we have dubbed PASTA (Planetary Atmospheres and Stellar RoTation RAtes), which runs within a Bayesian framework to model the atmospheric evolution of exoplanets. The code combines MESA stellar evolutionary tracks, a model describing planetary structures, a model relating stellar rotation and activity level, and a model predicting planetary atmospheric massloss rates based on the results of hydrodynamic simulations.
Results. Through a Markov chain Monte Carlo scheme, we retrieved the posterior probability density functions of all considered parameters. For ages older than about 2 Gyr, we find a median spindown (i.e. P(t)∝t^{y}) of ȳ = 0.38_{−0.27}^{+0.38}, indicating a rotation decay slightly slower than classical literature values (≈0.5), though still within 1σ. At younger ages, we find a median spindown (i.e. P(t)∝t^{x}) of x̄ = 0.26_{−0.19}^{+0.42}, which is below what is observed in young open clusters, though within 1σ. Furthermore, we find that the x probability distribution we derived is skewed towards lower spindown rates. However, these two results are likely due to a selection bias as the systems suitable to be analysed by PASTA contain at least one planet with a hydrogendominated atmosphere, implying that the host star has more likely evolved as a slow rotator. We further look for correlations between the initial atmospheric mass fraction of the considered planets and system parameters (i.e. semimajor axis, stellar mass, and planetary mass) that would constrain planetary atmospheric accretion models, but without finding any.
Conclusions. PASTA has the potential to provide constraints to planetary atmospheric accretion models, particularly when considering warm subNeptunes that are less susceptible to mass loss compared to hotter and/or lowermass planets. The TESS, CHEOPS, and PLATO missions are going to be instrumental in identifying and precisely measuring systems amenable to PASTA’s analysis and can thus potentially constrain planet formation and stellar evolution.
Key words: planets and satellites: atmospheres / stars: activity / stars: rotation
© ESO 2021
1. Introduction
A variety of planetfinding facilities have led to the detection of more than 4000 exoplanets and an almost equal amount of candidates^{1}. Such a large number of planets has allowed us to form a picture of the galactic planet population and of the underlying structures, which are believed to be caused by a combination of formation, migration, and atmospheric evolution processes. For example, the socalled subJovian desert (e.g. Davis & Wheatley 2009; Szabó & Kiss 2011; Mazeh et al. 2016) and subNeptune radius gap (e.g. Fulton et al. 2017; Fulton & Petigura 2018; Van Eylen et al. 2018; MacDonald 2019) are among the most noticeable structures in the observed exoplanet population.
The radius gap and the lower edge of the subJovian desert are currently believed to be the consequence of planetary atmospheric evolution, and in particular of atmospheric escape, though in case of the radius gap the nature of the main escapedriving mechanism is still unclear (e.g. Owen & Wu 2017; Owen & Lai 2018; Jin & Mordasini 2018; Ginzburg et al. 2018; Gupta & Schlichting 2019, 2020; Loyd et al. 2020; Sandoval et al. 2021). In general, it is believed that atmospheric escape plays a pivotal role in shaping the longterm evolution of planetary atmospheres. Therefore, adequately tracking planetary atmospheric evolution gives us the possibility to further understand the observed planet population as well as gather critical information about planet formation (e.g. Jin et al. 2014; Kubyshkina et al. 2019b).
Atmospheric escape is particularly strong when planets are young and is driven by the highenergy – Xray and extreme UV (EUV; together XUV) – emission of the host star (blowoff; e.g. Watson et al. 1981; Lammer et al. 2003) and/or by the internal atmospheric energy and low planetary gravity (boiloff; e.g. Stökl et al. 2015; Ginzburg et al. 2016; Owen & Wu 2017; Fossati et al. 2017). For computational reasons, in atmospheric evolution calculations, escape is computed employing analytical formulas (e.g. the energylimited approximation; Watson et al. 1981; Erkaev et al. 2007) that, however, are unreliable, particularly in the regime that most matters for atmospheric evolution (e.g. Kubyshkina et al. 2018a,b; Krenn et al. 2021). To circumvent this problem, Kubyshkina et al. (2018a) present an analytical expression that enables massloss rates to be computed for a variety of planets on the basis of the large grid of hydrodynamic upper atmosphere models published by Kubyshkina et al. (2018b), thus enabling both boiloff and blowoff to be simultaneously accounted for.
Estimating the evolution of a planetary atmosphere is further complicated by the fact that the evolution of the rotation rate, and thus of the highenergy emission, of latetype stars is not unique. As a matter of fact, stellar rotation rate and XUV emission are tightly bound, with faster rotating stars being XUV brighter (e.g. Pallavicini et al. 1981; Pizzolato et al. 2003; Johnstone et al. 2015a). Along the main sequence, both rotation rate and XUV emission decrease with time, but their evolutionary path is not unique: After the saturation regime, stars born with the same properties (i.e. mass and metallicity) can have widely different initial rotation rates and XUV fluxes (e.g. Mamajek & Hillenbrand 2008; Denissenkov et al. 2010; Johnstone et al. 2015a, 2021). The evolutionary paths of both rotation rate and XUV emission merge at an age of about 2 Gyr (e.g. Tu et al. 2015), meaning that for older stars it is not possible to infer the past rotation rate and XUV emission from observations alone. To overcome this problem, Kubyshkina et al. (2019a,b) developed an algorithm that enables the evolution of a hydrogendominated planetary atmosphere to be tracked accounting for atmospheric escape using massloss rates extracted from the grid of hydrodynamic models from Kubyshkina et al. (2018b). The algorithm employs a Bayesian framework and the currently observed system parameters to simultaneously constrain the past evolution of the planetary atmosphere and of the stellar rotation rate, and hence XUV emission. In practice, the algorithm returns a posterior probability density function (PDF) for the planetary initial atmospheric mass fraction and for the rotation rate of the host star at any desired age.
In its original implementation the algorithm had been applied to the HD 3167, K232, Kepler11, and ν^{2} Lup systems, constraining their past evolution (Kubyshkina et al. 2019a,b; Delrez et al. 2021). However, the code relied on conversion parameters (e.g. between the stellar rotation rate and Xray emission, and between Xray and EUV flux) taken from the literature, which suffer from significant uncertainties. We present here a major update to the code, which we have dubbed PASTA (Planetary Atmospheres and Stellar RoTation RAtes), where we also treat these conversion parameters as free parameters and refine the gyrochronological relation. We then apply the algorithm to a large number of systems for which the masses and radii for at least one planet are known to better than 25% and 8%, respectively. Our code is also able to deal with multiplanet systems, which are the best targets as the multiple fitting points provided by the different planets orbiting the same star increase the constraining power of the algorithm.
We aim at studying the stellar spindown after carefully checking the consistency of our adopted framework with the multiple and diverse constraints coming from observations, theory, and literature. We also search for differences between the distribution of rotation rates of young planet hosts with that of young open cluster (OC) stars of comparable mass and age that may hint at the presence of starplanet interactions (SPIs) early on in the evolution of planetary systems. Finally, we look for correlations between the derived planetary initial atmospheric mass fractions and system parameters (e.g. initial atmospheric mass fraction vs. either planetary or stellar mass and semimajor axis) that would enable planetary accretion models, and thus planet formation to be empirically constrained (e.g. Lozovsky et al. 2021).
This paper is organised as follows. Section 2 describes the whole framework of the PASTA code, while Sect. 3 presents the sample of exoplanetary systems analysed in this work. Our results are outlined in Sect. 4 and then thoroughly discussed in Sect. 5. Finally, our conclusions are drawn in Sect. 6.
2. The PASTA algorithm
We modelled the evolutionary history of both planetary atmospheres and stellar rotation rates within a Bayesian framework following a significant upgrade of PASTA, first presented by Kubyshkina et al. (2019a,b). Below, we thoroughly describe the entire algorithm and the implemented upgrades.
2.1. Stellar and planetary models
Modelling the evolution of an exoplanetary system requires the use of several different modelling tools and results. The various ingredients are stellar evolutionary tracks, a model of the stellar XUV flux evolution, a planetary structure model relating planetary parameters and atmospheric mass, and a model computing planetary atmospheric escape.
The stellar evolutionary tracks are those present in the MESA Isochrones and Stellar Tracks (MIST; Choi et al. 2016) grid, which has been computed employing the MESA stellar evolutionary code. Within the MIST grid, we consider tracks computed for solar metallicity stars with a mass ranging between 0.4 and 1.3 M_{⊙}, with steps of 0.05 M_{⊙} up to 0.9 M_{⊙} and steps of 0.02 M_{⊙} for highermass stars. Each evolutionary track lists the theoretically expected values of the stellar bolometric luminosity L_{bol}, effective temperature T_{eff}, and radius R_{⋆} at any given epoch t. These stellar parameters are fundamental to track the equilibrium temperature of the hosted exoplanets over time (Sect. 2.2).
The planetary structure model relating planetary mass M_{p}, radius R_{p}, equilibrium temperature T_{eq}, and atmospheric mass M_{atm} is that presented by Johnstone et al. (2015b). In this work, we employ the model grid described by Kubyshkina et al. (2018b), which spans the 1–40 M_{⊕} range and thus covering superEarths and subNeptunes.
Variations in M_{atm} over time are computed considering the atmospheric massloss rate Ṁ_{atm}, which depends on the stellar XUV emission and planetary properties. PASTA uses the stellar rotation period P_{rot} as proxy for the stellar highenergy flux; therefore, it is essential to reconstruct the evolution of P_{rot} with time t: P_{rot}(t), where the variable of temporal evolution t may span the range from the initial reference epoch (e.g. the time of protoplanetary disk dispersal, t_{disk}) up to the stellar age t_{⋆}. For young stars, that is stars with an age t_{⋆} < 2 Gyr, we modelled P_{rot}(t) as a power law of the form
$$\begin{array}{c}\hfill {P}_{\mathrm{rot}}(t)={P}_{\mathrm{rot},\star}{\left(\frac{t}{{t}_{\star}}\right)}^{x}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}{t}_{\star}<2\phantom{\rule{0.166667em}{0ex}}\mathrm{Gyr}\end{array}$$(1)
where P_{rot,⋆} is the presentday rotation period of the star, while x is a free parameter varying within the [0, 2] range to account for the different rotation rates observed for young latetype stars (see e.g. Tu et al. 2015, Fig. 2).
Instead, for stars older than 2 Gyr, the picture is more complicated. Despite several gyrochronological relations being available in the literature (e.g. Barnes 2003, 2010; Mamajek & Hillenbrand 2008; Collier Cameron et al. 2009), all the periodage relations are calibrated upon OCs, whose age has been previously estimated through a global isochrone fitting as all cluster members are coeval. However, there is a lack of OCs for calibrating gyrochronological laws in the intermediate and oldage domain and quite some doubts have been raised when applying gyrochronology to field stars. For example, Barnes (2009), Brown (2014), and Kovács (2015) found that for field stars isochronal ages are considerably greater than gyrochronological ages. In particular, Kovács (2015) stresses that field stars appear to have significantly lower slowdown rates compared to their cluster counterparts and van Saders et al. (2016) shows that fast rotators are unexpectedly found among stars more evolved than the Sun. A complete review of the reliable domains of application of classical gyrochronological relations is beyond the scope of this work, but the above considerations brought us to adopt a P_{rot}(t) relation, which assume the following form when a star is older than 2 Gyr (ages in Gyr):
$$\begin{array}{c}\hfill {P}_{\mathrm{rot}}(t)=\{\begin{array}{cc}{P}_{\mathrm{rot},\star}{\left(\frac{t}{2}\right)}^{x}{\left(\frac{2}{{t}_{\star}}\right)}^{y}\hfill & t<2\phantom{\rule{0.166667em}{0ex}}\mathrm{Gyr}\hfill \\ {P}_{\mathrm{rot},\star}{\left(\frac{t}{{t}_{\star}}\right)}^{y}\hfill & t\ge 2\phantom{\rule{0.166667em}{0ex}}\mathrm{Gyr}\hfill \end{array}\phantom{\rule{1em}{0ex}}{t}_{\star}\ge 2\phantom{\rule{0.166667em}{0ex}}\mathrm{Gyr},\end{array}$$(2)
where the broken powerlaw enables us to differentiate between two evolutionary regimes, whose break is set at t = 2 Gyr. In fact, stars even of the same spectral type may have different P_{rot} values at t = t_{disk} and hence evolve along different evolutionary rotation tracks, which then tend to converge towards similar P_{rot} values at t ≈ 2 Gyr, where P_{rot}(t = 2 Gyr) depends on stellar mass. As a star may be characterised by highly different spindown rates in the first part of its life (t < 2 Gyr), the first regime of Eq. (2) is governed by the x exponent spanning the [0, 2] range as it happens for Eq. (1). Instead, the factor ${\left(\frac{2}{{t}_{\star}}\right)}^{y}$ allows the function to be continuous, as there is not dependence on the variable t. The following stellar rotation evolution (t > 2 Gyr) is generally quieter; thus, we introduced an additional y exponent, which varies within the [0.01, 1] range, to model the stellar spindown rate. The value of the y exponent is centred around $\frac{1}{2}$, which is the value adopted by the classical Skumanichlaw (Skumanich 1972). Also, we remark that this value is very close to 0.566, which is the value proposed by Mamajek & Hillenbrand (2008) in their powerlaw gyrochronological relation. Equation (2) is set up so that P_{rot}(t = t_{⋆}) = P_{rot,⋆} and it guarantees continuity of the function at t = 2 Gyr. To further ensure that P_{rot}(t) is differentiable at t = 2 Gyr, within the range t ∈ [1.5, t_{lim}], where 2 ≤ t_{lim} = min{2.5, t_{⋆}}, we replaced the broken powerlaw with a 3rddegree single spline fitted onto the powerlaw.
The P_{rot}(t) function does not explicitly include tidal effects that may be induced by the planets on the host star. However the flat and broad priors on the gyrochronological exponent allow us to span between slow and fast spindown rates. The occurring of SPI is mentioned as a possible scenario in Sect. 5.2.
Given P_{rot} at any given time, the Xray stellar luminosity L_{X} (λ = 5–100 Å) can be estimated using scaling relations. Indeed, L_{X}, which traces stellar magnetic activity, can be expressed as a function of Rossby number
$$\begin{array}{c}\hfill \mathrm{Ro}=\frac{{P}_{\mathrm{rot}}}{\tau},\end{array}$$(3)
which probes the efficiency of the stellar magnetic dynamo (see e.g. Noyes et al. 1984). τ is the convective turnover time, which is massdependent and we expressed it as (Wright et al. 2011)
$$\begin{array}{c}\hfill log\tau =1.161.49log\frac{M}{{M}_{\odot}}0.54{log}^{2}\frac{M}{{M}_{\odot}}\xb7\end{array}$$(4)
For increasing stellar rotational velocities, L_{X} increases monotonically (Pallavicini et al. 1981) following different massdependent tracks (McDonald et al. 2019). This behaviour breaks down for very fast rotators (e.g. Micela et al. 1985), so that below a given threshold value of P_{rot} (hence of Ro, say Ro_{sat}) L_{X} saturates, that is L_{X} becomes Roindependent, while keeping the dependence on M_{⋆}: L_{X}(Ro < Ro_{sat}) = L_{X, sat}(M_{⋆}). Therefore, following Wright et al. (2011), we modelled the ratio between L_{X} and the stellar bolometric luminosity L_{bol} as
$$\begin{array}{c}\hfill {R}_{\mathrm{X}}\equiv \frac{{L}_{\mathrm{X}}}{{L}_{\mathrm{bol}}}=\{\begin{array}{cc}{R}_{\mathrm{X},\mathrm{sat}}\hfill & \mathrm{Ro}\le {\mathrm{Ro}}_{\mathit{sat}}\hfill \\ C{\mathrm{Ro}}^{\beta}\hfill & \mathrm{Ro}>{\mathrm{Ro}}_{\mathit{sat}}\hfill \end{array}\end{array}$$(5)
where Ro_{sat} = 0.13 is the saturation threshold of Ro as quantified by Wright et al. (2011), while β = −2.70 ± 0.13 has been estimated by Wright et al. (2011) from an unbiased sample of latetype stars for which both P_{rot} and L_{X} had been measured. As R_{X, sat} ≡ L_{X, sat}/L_{bol} depends on stellar mass, rather than fixing it to the bestfit value as in Wright et al. (2011), we evaluated R_{X, sat} for different mass ranges from McDonald et al. (2019, Fig. 1), adopting the values that are reported in Table 1. For each mass range, the C coefficient in Eq. (5) is then computed accordingly to guarantee the continuity of R_{X} at Ro = Ro_{sat}. For reference, the C values for the average β = −2.70 are reported in Table 1, but we remark that the default behaviour of our tool is that β varies according to a Gaussian prior whose σ = 0.13. Therefore, C depends also on β.
Saturation values (R_{X, sat}) adopted in our work for different mass ranges, as inferred from McDonald et al. (2019).
From L_{X}, we infer the EUV (λ = 100–920 Å) stellar luminosity (L_{EUV}) via the relation reported by SanzForcada et al. (2011)
$$\begin{array}{c}\hfill \begin{array}{cc}\hfill log{L}_{\mathrm{EUV}}& =({q}_{\mathrm{SF}}\pm {\sigma}_{q})+({m}_{\mathrm{SF}}\pm {\sigma}_{m})log{L}_{\mathrm{X}}\hfill \\ \hfill & =(4.80\pm 1.99)+(0.860\pm 0.073)log{L}_{\mathrm{X}},\hfill \end{array}\end{array}$$(6)
where both Xray and EUV luminosities are expressed in erg s^{−1}. The term L_{EUV} is the keyingredient necessary to finally compute Ṁ_{atm}. To this end, we employed the hydrobased approximation (HBA) presented by Kubyshkina et al. (2018a), which is an analytical formulation of Ṁ_{atm} extracted by fitting the results of the large grid of upper planetary atmosphere hydrodynamic models presented by Kubyshkina et al. (2018b). The HBA speeds up decisively the computation of Ṁ_{atm} compared to using hydrodynamic simulations and, at the same time, it removes the physical limitations of the energylimited approximations (e.g. Krenn et al. 2021). In detail, Kubyshkina et al. (2018a) computes Ṁ_{atm} following two different functional behaviours, which depend on the restricted Jeans escape parameter of a planet (Jeans 1925; Fossati et al. 2017)
$$\begin{array}{c}\hfill \mathrm{\Lambda}=\frac{G{M}_{\mathrm{p}}{m}_{\mathrm{H}}}{{R}_{\mathrm{p}}{k}_{\mathrm{B}}{T}_{\mathrm{eq}}},\end{array}$$(7)
where G is the universal gravitational constant, k_{B} is the Boltzmann constant, and m_{H} is the hydrogen mass. In particular, Kubyshkina et al. (2018a) express the atmospheric massloss rate (in g s^{−1}) of a given planet as a function of semimajor axis a and R_{p} as
$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{atm}}={e}^{{\alpha}_{0}}{\left(\frac{{F}_{\mathrm{EUV}}}{\phantom{\rule{0.333333em}{0ex}}\text{erg}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.333333em}{0ex}}{\text{cm}}^{2}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.333333em}{0ex}}{\text{s}}^{1}}\right)}^{{\alpha}_{1}}{\left(\frac{a}{\mathrm{AU}}\right)}^{{\alpha}_{2}}{\left(\frac{{R}_{\mathrm{p}}}{{R}_{\oplus}}\right)}^{{\alpha}_{3}}{\mathrm{\Lambda}}^{K},\end{array}$$(8)
where e is Euler’s number and ${F}_{\mathrm{EUV}}=\frac{{L}_{\mathrm{EUV}}}{4\pi {a}^{2}}$ is the EUV stellar flux at the distance of the planet, while the exponent K is given by
$$\begin{array}{c}\hfill lnK=\zeta +\theta ln\frac{a}{\mathrm{AU}}\xb7\end{array}$$(9)
The values of the exponents α_{0}, α_{1}, α_{2}, α_{3}, and K (through ζ and θ) are listed in Table 2 for convenience and vary according to whether lnΛ is smaller or larger than the threshold value e^{Σ}, where
$$\begin{array}{c}\hfill \mathrm{\Sigma}=\frac{15.6110.578ln\frac{{F}_{\mathrm{EUV}}}{\phantom{\rule{0.333333em}{0ex}}\text{erg}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.333333em}{0ex}}{\text{cm}}^{2}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.333333em}{0ex}}{\text{s}}^{1}}+1.537ln\frac{a}{\mathrm{AU}}+1.018ln\frac{{R}_{p}}{{R}_{\oplus}}}{5.564+0.894ln\frac{a}{\mathrm{AU}}}\xb7\end{array}$$(10)
Coefficients of Eq. (8) taken from Kubyshkina et al. (2018a, Table 1).
Equation (8) enables the evolution of the atmospheric mass with time to then be computed.
2.2. The statistical framework
PASTA works within a Bayesian framework and can take the following jump parameters: (i) the y exponent of the gyrochronological relation; (ii) the stellar rotation period at an age of 150 Myr P_{rot,150}, where 150 Myr is a representative early stage of stellar evolution that enables us to directly compare the posterior PDF obtained for a given star with the respective distribution gathered from measurements of the rotation period of OC stars (e.g. Johnstone et al. 2015a; see below for the reason for which we decided to use P_{rot,150} instead of x as jump parameter); (iii) the β exponent of Eq. (5); (iv) the q_{SF} and m_{SF} coefficients of Eq. (6); (v) the epoch of dispersal of the protoplanetary disk t_{disk}; (vi) t_{⋆}, P_{rot, ⋆}, and M_{⋆}; and (vii) for each planet, a, M_{p}, and ${f}_{\text{atm}}^{\text{start}}={f}_{\mathrm{atm}}({t}_{\mathrm{disk}})={\frac{{M}_{\mathrm{atm}}}{{M}_{\mathrm{p}}}}_{t={t}_{\mathrm{disk}}}$, which is the atmospheric mass fraction evaluated at the beginning of the evolution, that is, for t = t_{disk}.
PASTA starts the evolution after the dispersal of the protoplanetary disk, whose timescale is estimated to be a few megayears (see e.g. Montmerle et al. 2010; Alexander et al. 2014; Kimura et al. 2016; Gorti et al. 2016, and references therein). Therefore, although t_{disk} can be set as a jump parameter, for this work we fix it at t_{disk} = 5 Myr. We will explore the capability of PASTA to also constrain t_{disk} in a future work.
The default setup is imposing Normal (Gaussian) priors 𝒩(μ, σ) on t_{⋆}, P_{rot, ⋆}, M_{⋆}, a, and M_{p}, which are obtained through observations and stellar evolutionary models. Normal priors are imposed also on β, q_{SF}, and m_{SF} having 𝒩(−2.70, 0.13), 𝒩(4.80, 1.99), and 𝒩(0.860, 0.073), respectively.
Instead, the parameters y, P_{rot,150}, and ${f}_{\text{atm}}^{\text{start}}$ are set having uniform priors. In particular, a flat prior on P_{rot,150} implies a uniform sampling of values within the specified P_{rot,150}range and avoids biasing the sampling towards low P_{rot,150} values, which would happen if a uniform prior on x was assumed. In practice, the output of the code are posterior PDFs for P_{rot,150} and ${f}_{\text{atm}}^{\text{start}}$ that are constrained by the currently observed system parameters and according to the considered input models.
As a statistical framework, PASTA uses a Markov chain Monte Carlo (MCMC) scheme implemented through the MC3 tool (Cubillos et al. 2017). At each chain step PASTA samples a set of input parameters from the priors and it first selects the proper grid of stellar models according to the step value of M_{⋆}, so to infer L_{bol}, T_{eff}, and R_{⋆} by interpolation within the stellar evolutionary tracks at the starting epoch of evolution t_{disk}. Then, for each considered planet in the system, PASTA computes the equilibrium temperature as
$$\begin{array}{c}\hfill {T}_{\mathrm{eq}}={T}_{\mathrm{eff}}\sqrt{\frac{{R}_{\star}}{2a}},\end{array}$$(11)
which assumes zero albedo and a moderately rotating planet, so that its entire spherical surface reirradiates the flux received by the host star. Then, T_{eq} and the step values of M_{p} and ${f}_{\text{atm}}^{\text{start}}$ enable interpolation within the planetary structure model grid to obtain the initial planetary radius ${\widehat{R}}_{p}({t}_{\mathrm{disk}})$. At this point, from the step value of P_{rot} computed from Eq. (1) or (2), PASTA estimates F_{EUV} using the scaling relations given in Eqs. (5) and (6), so that Ṁ_{atm} can be computed through Eq. (8).
The code continuously increases the evolutionary age by Δt, which is adjusted such that the atmospheric mass loss is less than 5% of the entire atmospheric mass M_{atm}. At the end of each MCMC step the code has generated a planetary evolutionary track by updating M_{atm} (from Ṁ_{atm} and the time step Δt), and hence the planetary radius ${\widehat{R}}_{p}({t}_{\mathrm{disk}}+\mathrm{\Delta}t)$. This loop terminates (i.e. the planetary evolutionary track reaches its end) when either the evolution reaches the present day stellar age t_{⋆} or the atmosphere is fully lost. The theoretical presentday planetary radius ${\widehat{R}}_{p}({t}_{\star})$ is finally compared with the observed value R_{p} and that specific MCMC step is accepted according to the MetropolisHastings acceptance rule discussed in Cubillos et al. (2017).
The algorithm uses the planetary radii as observational constraints to drive the chains’ convergence. An alternative option implemented in the code is to track the atmospheric mass fraction ${\widehat{f}}_{\mathrm{atm}}(t)$, rather than ${\widehat{R}}_{p}(t)$, to finally compare ${\widehat{f}}_{\mathrm{atm}}({t}_{\star})$ with its observational counterpart f_{atm}(t_{⋆}) that may be obtained through a planetary atmospheric structure model and the planetary basic parameters (e.g. Rogers et al. 2011; Lopez & Fortney 2014; Dorn et al. 2017; Delrez et al. 2021).
The final results are posterior PDFs for each considered parameter. It is crucial to check the priorposterior consistency of those parameters for which normal priors were imposed. If those priors are respected, the posterior PDFs of the parameters for which uniform priors were set (P_{rot,150} and ${f}_{\text{atm}}^{\text{start}}$ in this case) may be considered physically reliable or, at least, consistent with the adopted framework.
As already highlighted by Kubyshkina et al. (2019b), the constraining power of the code increases with increasing number of planets in the system. This is because the code has multiple fitting points (i.e. the radius of each planet in the system), but at the same time all planets in the system orbit around the same star and thus simultaneously constrain P_{rot,150}. Furthermore, the code can also be employed to constrain planetary masses in case the values are either poorly constrained or unavailable, as shown for example by Kubyshkina et al. (2019a) and Bonfanti et al. (2021).
The reliability of PASTA’s results depend on the reliability and accuracy of the input system parameters, of the background models, and of the assumptions. The scheme described above is based on two main assumptions: (i) the considered planets either hosted in the past or still host a hydrogendominated atmosphere (see Owen et al. 2020, for a thorough discussion on the general validity of this assumption); and (ii) planet migration happened inside the disk and the planetary orbital separation does not change following disk dispersal (i.e. a is constant over time). Furthermore, the planetary structure models we consider assume that planetary cores have an Earthlike density, regardless of the planetary mass. Lopez & Fortney (2014) and Petigura et al. (2016) validated this assumption showing that, for planets with a hydrogendominated atmosphere, the planetary radius is generally independent of the assumed core composition. However, a result of this assumption is that the algorithm works appropriately only for planets with an average density ρ_{p} ≤ 1ρ_{⊕}.
3. The sample
We compiled the sample of systems analysed in this work on the basis of the catalogue presented by Otegi et al. (2020), which contains systems for which for at least one planet the uncertainties on M_{p} and R_{p} are smaller than 25% and 8%, respectively. Within this sample, we further selected just systems hosting planets with 5 ≲ M_{p}/M_{⊕} ≲ 30 in order to remain within the boundaries of the planetary structure and escape model grids we have available. Furthermore, we considered exclusively multiplanet systems, to benefit from the simultaneous multiple constrain on P_{rot,150}.
Following the aforementioned criteria, we ended up with a dozen systems. The planetary input parameters required to run PASTA, and their sources, are listed in Table 3 (first five columns). In most cases, we took the same sources as chosen by Otegi et al. (2020), further complementing them with additional sources in case of missing information. For Kepler11, given the several and often not consistent M_{p} values present in the literature, we decided to combine the solutions provided by different authors, in an attempt to increase the robustness of the considered planetary masses. In practice, we treated each estimated ${M}_{p{\sigma}_{\text{low}}}^{+{\sigma}_{\text{up}}}$ as a Normal PDF with possibly asymmetric tails and we summed all PDFs referring to the same planet. We then extracted the mode and estimated the asymmetric uncertainties from the resulting distribution.
Planetary radii R_{p}, masses M_{p}, and semimajor axes a taken from the literature and adopted in this work, followed by the main results obtained with PASTA (median values and 1σ confidence interval), that is the stellar rotation period at 150 Myr (P_{rot,150}), the theoreticallyestimated planetary mass (${\widehat{M}}_{\mathrm{p}}$), and the atmospheric mass fraction at t = t_{disk} = 5 Myr (${f}_{\text{atm}}^{\text{start}}$).
Besides the planetary parameters, PASTA also requires M_{⋆}, t_{⋆}, and P_{rot,⋆} as input. We obtained both stellar mass and age by homogeneously analysing the host stars using the isochrone placement algorithm described in Bonfanti et al. (2015, 2016). Basic input parameters of the isochrone placement were the stellar effective temperature T_{eff,⋆}, metallicity [Fe/H]_{⋆}, gravity log g_{⋆}, Gaia magnitude G_{⋆}, and the Gaia parallax π_{⋆}. When available from observations, the stellar rotation period P_{rot,⋆}, the projected rotational velocity v sin i_{⋆} and the chromospheric activity index log ${R}_{\text{HK}}^{\prime}$ were further added to the basic input set of the isochrone placement to possibly remove isochronal degeneracies, to improve convergence. Following Bonfanti et al. (2021), we finally enlarged the internal uncertainties on M_{⋆} and t_{⋆} by adding in quadrature 4% and 20%, respectively, to account for isochrones’ systematics.
Any observational P_{rot,⋆} value imposes a further constraint during the preliminary stellar age derivation process within the isochrone placement algorithm, but it is an optional parameter. Conversely, the following application of our PASTA algorithm always requires an estimate for P_{rot,⋆}. Among the considered sample of systems, just Kepler11, Kepler411, and WASP47 have published P_{rot} values. For the other stars, we estimated their current rotation periods by numerically inverting the gyrochronological relation of Barnes (2010) as their ages are known from the evolutionary isochrones and tracks via the isochrone placement technique.
All stellar input parameters are listed in Table 4. Except for the M_{⋆}, t_{⋆}, and most of the P_{rot,⋆} values, which we computed as described above, all the other listed parameters were taken from the same sources specified in Table 3.
Stellar parameters.
4. Results
PASTA returns posterior distributions for each jump parameter. As discussed above, the main results are the posterior distributions of P_{rot,150} and ${f}_{\text{atm}}^{\text{start}}$, which were set with uniform priors, but to assess their reliability it is necessary to first check the priorposterior agreement for those parameters for which a Normal prior based on observations was imposed.
As an example illustrating the capabilities of PASTA, we present in Figs. 1 and 2 the priors and posteriors of all considered jump parameters for the K2285 planetary system, which contains four planets. Similar plots, but for all other considered systems, are shown in Appendix A.
Fig. 1. Stellar properties of the K2285 system given by PASTA. First row: P_{rot,150}, t_{⋆}, P_{rot, ⋆}, and M_{⋆}. Second row: L_{X} at 150 Myr, and the x and y exponents adopted by the gyrochronological relations given by Eqs. (1) and (2). Third row: q_{SF} and m_{SF} coefficients of Eq. (6), and the β exponent of Eq. (5). The output posterior PDFs are shown as blue histograms, with the green area representing the 68.3% HPD credible interval. The red curves are the Gaussian or flat (uninformative) imposed priors. When absent (second row), no specific prior was imposed on the parameter; in particular, x and L_{X} are mathematically derived from the PDFs of the jump parameters. The black histogram in the topleftmost panel represents the P_{rot,150} distribution extracted from the sample of Johnstone et al. (2015a) considering stars with a mass differing by less than 0.1 M_{⊙} from the mass of K2285. 
Fig. 2. PDFs given by PASTA for the considered planetary parameters in the K2285 system: semimajor axis (first row), mass (second row), and ${f}_{\text{atm}}^{\text{start}}$ (third row). The black curves in each panel of the third row (evident in the second panel, and just barely visible in the others as they are squeezed towards zero) indicate the predicted presentday atmospheric mass fraction. It follows that K2285 b, d, and e have almost entirely lost their atmospheres at some point in time along the system’s evolution. 
4.1. The exoplanetary system K2285
In the following, we aim to guide the reader through the interpretation of the plots shown in Figs. 1 and 2, and in turn in Appendix A. Starting from Fig. 1, the posterior distributions (blue histograms) of t_{⋆}, P_{rot, ⋆}, and M_{⋆} (first row, second to fourth panel) nicely agree with their respective priors (red Gaussians). The very slight tension between priors and posteriors of q_{SF} and m_{SF} (third row, first and second panel) highlights the importance of setting them as jump parameters rather than fixing them (see also Sect. 5.1). In fact, giving enough degrees of freedom to the algorithm allows it to find the best match between theoretical and observational constraints.
The global priorposterior agreement, despite the large number of interacting models composing the algorithm, suggests that the PDF of P_{rot,150} unveiling the past rotation history of K2285 (first row, leftmost panel) can be considered to be reliable within the given assumptions and gives a median value of ${\overline{P}}_{\mathrm{rot},150}=15.{3}_{9.4}^{+13.5}$ days. A comparison of the PDF of P_{rot,150} with the P_{rot,150} distribution obtained from considering stars having M_{J15} − M_{⋆, K2285}< 0.1 M_{⊙} within the sample of Johnstone et al. (2015a) indicates that, when it was young, K2285 was a slow rotator.
The second row of Fig. 1 contains three panels. The leftmost panel displays the estimated L_{X} at 150 Myr, which is derived through Eq. (5) and considering P_{rot,150}, while the other two panels show the PDFs of the gyrochronological exponents x and y. Although characterised by a well defined peak, the PDF of L_{X} presents a little bump at log L_{X} > 29.5, which results because L_{X} saturates at low P_{rot} values (saturation regime). As a matter of fact, L_{X} and P_{rot} are anticorrelated, but in the regime of very fast rotators (low P_{rot}), L_{X} stops increasing and becomes equal to L_{X, sat}, which depends only on stellar mass (see e.g. Wright et al. 2011, Fig. 2). Therefore, towards large values, the PDF of L_{X} stops decreasing smoothly and all L_{X, sat} values are lumped in the last bins towards high L_{X} values. In some cases, this bin is so pronounced to enter the highest posterior density (HPD) credible interval, but the HPD credible interval for L_{X} should be continuous. This means that in these cases the high L_{X} solution should be considered to be a bias driven by the lumping of high L_{X} values in a narrow range independently of P_{rot}.
The median of the x distribution $\overline{x}=0.{23}_{0.16}^{+0.33}$ points towards a slow stellar spindown if compared to the spindown rates theoretically expected by Tu et al. (2015, Fig. 2). In fact, the median $\overline{x}$ we derived for K2285 corresponds to the 10th percentile of the theoretical distribution. Moreover, the HPD credible interval of y spans values below 0.50 and the median of the distribution is $\overline{y}=0.{32}_{0.24}^{+0.37}$, which is smaller (though within 1σ) than the classical Skumanich exponent that is equal to $\frac{1}{2}$. We refer the reader to Sects. 5.1 and 5.2 for a deeper discussion about the spindown rate of our targets.
Figure 2 shows the PDFs of the semimajor axis a (first row), mass M_{p} (second row), and initial atmospheric mass fraction ${f}_{\text{atm}}^{\text{start}}$ (third row) for each planet detected in the K2285 system. The Normal priors imposed on a are respected by the posteriors. Normal priors were also imposed on the masses of planets b and c (second row, two leftmost panels) as derived from radial velocity, while uniform priors were imposed on planets d and f, which have just upper mass limits (Palle et al. 2019).
Figure 2 shows that PASTA is unable to constrain ${f}_{\text{atm}}^{\text{start}}$ for planets b, d, and e, whose posteriors are essentially flat. In other words, any atmospheric mass fraction at t = t_{disk} is compatible with their presentday atmospheric content, which is negligible. Instead, for planet c PASTA gives a current atmospheric mass fraction of ${f}_{\mathrm{atm},\mathrm{c}}^{\mathrm{now}}$ ≈ 0.025 and returns a well constrained initial atmospheric mass fraction of ${\overline{f}}_{\mathrm{atm},\mathrm{c}}^{\mathrm{start}}=0.{0384}_{0.0053}^{+0.0052}$. Therefore, it is the constraint imposed by planet c that mainly drives the atmospheric evolution of the K2285 system, defining the activity level of the host star over time. As a positive side effect of studying multiplanet system, despite the uniform priors and as a result of planet c constraining the evolution of the stellar rotation rate, PASTA returns a rather tight constraint on the masses of planets d and e, for which we obtain $3.{59}_{0.34}^{+0.33}\phantom{\rule{0.166667em}{0ex}}{M}_{\oplus}$ and $2.{04}_{0.18}^{+0.20}\phantom{\rule{0.166667em}{0ex}}{M}_{\oplus}$, respectively.
The K2285 corner plot of all the relevant jump parameters within our framework is shown in Fig. 3. It emphasises the overall lack of correlation between the majority of the jump parameters with a few exceptions. For example, q_{SF} and m_{SF} are clearly anticorrelated. The reason is that the data point scatter quantitatively described by Eq. (6) occupies a broad strip in the log L_{EUV}log L_{X} plane, which is quantified by the high σ_{q}. By decreasing the intercept q_{SF} and increasing the slope m_{SF} at the same time, the bestfit line still spans the same strip.
Fig. 3. Corner plot of the relevant jump parameters returned by PASTA for the K2285 system. 
Another anticorrelation involves ${f}_{\mathrm{atm},\mathrm{c}}^{\mathrm{start}}$ and M_{c}. Planet c is the only one characterised by a tight constraint on its ${f}_{\text{atm}}^{\text{start}}$. The higher M_{c}, the higher the gravitational potential that retains the atmosphere; therefore, PASTA estimates a lower value for ${f}_{\mathrm{atm},\mathrm{c}}^{\mathrm{start}}$ to still match the presentday atmospheric content given the lower massloss rate. Finally there are positive correlations involving the masses of planets b, d, and e, which are the ones for which PASTA was not able to constrain their respective ${f}_{\text{atm}}^{\text{start}}$. As planet d and e have no observational constraints on their masses, it is likely that the mass correlations reflect correlations present within our planetary structure models, when only T_{eq} drives the mass posterior PDFs.
5. Discussion
We combined the results obtained from all 12 systems to preliminary check whether the results of our framework agree with the considered empirical relations and then derive the distributions of the gyrochronological exponents x and y (Sect. 5.1). Furthermore, as Johnstone et al. (2015a) report the observed rotation periods of OC stars of 150 and 600 Myr, we also compared the average distributions of P_{rot,150} and P_{rot,600} with those derived from stars member of young OCs to look for possible traces of SPI occurring during the first stages of evolution (Sect. 5.2). Finally, we looked for correlations between the derived ${f}_{\text{atm}}^{\text{start}}$ values and system parameters to constrain planetary atmospheric accretion models (Sect. 5.3).
5.1. The gyrochronological exponents
Kubyshkina et al. (2019a,b) presented the results of several tests performed on synthetic systems aiming at validating PASTA’s framework, but they did not enable a validation of the physical results. Because of the nature of the results provided by PASTA, it is not possible to use any real system to validate the output P_{rot,150} and ${f}_{\text{atm}}^{\text{start}}$ distributions. However, as our framework is built on several empirical relations and theoretical models, we checked a posteriori whether the results obtained for the different systems respect all observational, empirical, and theoretical constraints at the same time. For example, the posterior PDFs on M_{⋆} and t_{⋆} are driven by MESA theoretical models, while the priors imposed on q_{SF}, m_{SF}, and β have an empirical root.
The three panels of Fig. 4 show the merged posterior PDFs (blue histograms) for q_{SF}, m_{SF}, and β, which have been obtained by combining the posterior distributions coming from the analysis of each planetary system. The agreement between posteriors and priors (black Gaussians) of these parameters in addition to the overall priorposterior agreement of all the jump parameters seen in each analysis tells us that both the theoretical and empirical components of our framework lead to a consistent picture.
Fig. 4. PDFs (blue histograms) of those jump parameters entering as coefficients in Eqs. (5) and (6). Left and middle panels: intercept q_{SF} and slope m_{SF} of Eq. (6). Right panel: β exponent of Eq. (5). The median and the ±1σ bounds of each PDF are marked as vertical solid and dashed red lines, respectively. The black Gaussians represent the priors inferred from observations for each empirical parameter. 
This preliminary check strengthens the statistics we derived for the gyrochronological exponents. In fact, we also merged the x and y posterior distributions of all our selected targets producing the histograms shown in Fig. 5. We recall that the x exponent controls the rotation rate of the star during its first 2 Gyr of life. This rate may vary significantly from star to star, in fact rotational evolutionary tracks differing by as much as ∼1.5 orders of magnitude tend to converge towards the same P_{rot} after a couple of Gyr (Tu et al. 2015, their Fig. 2). Instead, the y exponent controls the decay of the stellar rotation rate at older ages.
Fig. 5. Top: PDF of the x exponent of our adopted gyrochronological relations expressed in Eqs. (1) and (2) built upon all the stellar targets analysed in this work. The vertical solid red line marks the median of the distribution, while the two vertical dashed red lines marks the σ levels corresponding to the 15.87th and 84.14th percentiles. Bottom: same as the top panel, but for the y exponent of Eq. (2). We note that the top panel is in log scale. 
The classical Skumanich (1972) law models the decay of P_{rot} of main sequence stars as a function of time through a powerlaw whose exponent is $\frac{1}{2}$. Other examples of gyrochronological relations in the form of powerlaws may be found in, for example, Barnes (2003, 2007), Mamajek & Hillenbrand (2008), and Collier Cameron et al. (2009). These works emphasise that P_{rot} may be predicted from knowledge of both stellar age and colour index, where the latter is taken as a proxy of stellar mass. The different authors still suggest an age dependence of the form P_{rot} ∝ t^{α}, where α = 0.5 (Barnes 2003), α = 0.5189 (Barnes 2007), α = 0.566 (Mamajek & Hillenbrand 2008), and α = 0.56 (Collier Cameron et al. 2009). If, on the one hand, the majority of literature sources tend to prefer α values slightly greater than 0.5, on the other hand, Meibom et al. (2009) confirm the empirical P_{rot} ∝ ${t}^{\frac{1}{2}}$ dependence for Gtype stars, but find a slower spindown rate (i.e. compatible with α < 0.5) for K stars. All in all, an exponent of $\sim \frac{1}{2}$ appears to be generally adequate to describe stellar spindowns for population studies, but slight variations may be expected on a star by star basis, especially for stars with a mass different from that of Gtype stars.
Figure 5 indicates that the reference (i.e. median) gyrochronological exponent for our host stars can be estimated as $\overline{y}=0.{38}_{0.27}^{+0.38}$ with a general preference for slow spindown rates, even if the 84.14th percentile ${\overline{y}}_{+1\sigma}=0.76$ demonstrates the heaviness of the right tail. The difference between our inferred $\overline{y}$ value and the y values reported in the literature is well below 1σ.
The preference towards lower spindown rates shown by our stellar sample may be the result of a selection bias. As a matter of fact, exoplanets are preferentially found around quiet stars, hence slow rotators. This is supported also by the x distribution with a median of $\overline{x}=0.{26}_{0.19}^{+0.42}$, which strongly favours a quiet evolution during the first stages of stellar evolution given its median. We thoroughly discuss this in Sect. 5.2. With reference to Tu et al. (2015, Fig. 2, left panel), the range of our x values defined by the 68.3% confidence interval is compatible with the set of rotational evolutionary tracks spanning the same percentile range, even if our x distribution is skewed towards lower values with $\overline{x}$ being approximately the rate displayed by the red track of Tu et al. (2015, 10th percentile of the predicted rotational distribution).
Finally, we stress that gyrochronological relations in the literature are usually calibrated on OC stars (as they are a coeval), which are generally young. Instead, although dependent by the stellar and planetary models in use, our approach – constrained by the presentday M_{atm} – may infer the rate of rotational decay for stars of any age. The lower spindown rates we found are consistent with the conclusions independently drawn by Kovács (2015) or van Saders et al. (2016) about noncluster field stars and give an original contribution to the gyrochronological studies of planethosting stars.
5.2. Rotation rates at young ages
We combined all the P_{rot,150} distributions we obtained from each analysed system finally obtaining the distribution shown by the blue histogram in Fig. 6. This distribution has a median value of ${\overline{P}}_{\mathrm{rot},150}=7.{9}_{5.3}^{+9.8}$ days, where the error bars are 1σ uncertainties. We further used the results of Johnstone et al. (2015a) to construct a P_{rot J, 150} distribution considering stars member of ∼150 Myr old OCs and having the same mass distribution as that of the planet hosts analysed here. This distribution is shown by a black line in Fig. 6 and has a median value of ${\overline{P}}_{\mathrm{rot}\phantom{\rule{0.166667em}{0ex}}\mathrm{J},150}=2.{2}_{1.4}^{+1.5}$ days, which is ∼1σ smaller than the median value derived for our sample of planethosting stars.
Fig. 6. Top: comparison between the P_{rot,150} PDF obtained by combining those distributions extracted from the analysed sample (blue line) and the P_{rot,150} distribution extracted from Johnstone et al. (2015a), considering stars member of ∼150 Myr old OCs and having the same mass distribution as the planethosting stars analysed in this work. The vertical solid red line marks the median value of the distribution shown by the blue line, while the two vertical dashed red lines indicate the 15.87th and 84.14th percentiles. Bottom: same as the top panel, but for the P_{rot,600} PDF. 
Figure 6 suggests that the rotation rate at an age of 150 Myr of the sample of planethosting stars analysed in this work is on average slower than what observed for the majority of stars. Because of the small number statistics it is not possible to conclude whether this is the result of SPI or not, but the most likely explanation for this difference is a selection bias. By construction, PASTA’s capability of constraining P_{rot,150} of a given star hinges on it hosting at least one planet with a hydrogendominated atmosphere that has been and/or is still affected by atmospheric escape. On average, it is easier to find such planets orbiting stars that when young were slow rotators, because if they had been fast rotators the planets would have more likely lost their hydrogen envelope. Similar considerations hold for P_{rot,600}, where ${\overline{P}}_{\mathrm{rotJ},600}=6.{3}_{1.5}^{+1.6}$ days is smaller than ${\overline{P}}_{\mathrm{rot},600}=11.{5}_{5.6}^{+10.9}$ by ∼1σ.
5.3. Constraining planetary atmospheric accretion models
One of the main aims of this work, and more in general of PASTA, is to enable one to look for correlations between the initial atmospheric mass fraction of closein planets and system parameters to constrain planetary formation and atmospheric accretion models. In particular, we focus here on analysing the possible correlations present between initial atmospheric mass fraction and semimajor axis, planetary mass, and stellar mass.
Figure 7 shows the median ${f}_{\text{atm}}^{\text{start}}$ values provided by PASTA as a function of semimajor axis, of planetary mass, and of stellar mass for the analysed systems. There seems to be a trend between ${f}_{\text{atm}}^{\text{start}}$ and a, in which planets with a ≲ 0.25 AU have a larger ${f}_{\text{atm}}^{\text{start}}$ compared to planets at larger a, a possible link between ${f}_{\text{atm}}^{\text{start}}$ and M_{p}, in which ${f}_{\text{atm}}^{\text{start}}$ decreases with increasing M_{p}, and a possible link between ${f}_{\text{atm}}^{\text{start}}$ and M_{⋆}, in which ${f}_{\text{atm}}^{\text{start}}$ increases with increasing M_{⋆}, but none of these correlations is statistically significant. In detail, because of small number statistics, particularly at large orbital separations, the trend between ${f}_{\text{atm}}^{\text{start}}$ and a has a coefficient of determination R^{2} = 0.25. However, analyses of additional systems detected and/or measured for example by TESS (Ricker et al. 2015), CHEOPS (Benz et al. 2021), or PLATO (Rauer et al. 2014) may shed light on the possible existence of such a correlation. This trend is weak also because of the large uncertainties of the ${f}_{\text{atm}}^{\text{start}}$ values obtained for the planets at short orbital separation. Indeed, for most of the planets with a ≲ 0.25 AU we obtained (almost) flat PDFs on ${f}_{\text{atm}}^{\text{start}}$. This may be a result of planetary atmospheric evolution as planets farther away from their host star are less subject to atmospheric escape, and thus have a slower and more unique evolution compared to closerin planets. The link between ${f}_{\text{atm}}^{\text{start}}$ and M_{p} is not significant (R^{2} = 0.12), but it shows the potential capability of PASTA to constrain atmospheric accretion processes. This plot further highlights that PASTA allows one to better constrain ${f}_{\text{atm}}^{\text{start}}$ for larger mass planets, as the vast majority of flat ${f}_{\text{atm}}^{\text{start}}$ posteriors is obtained for lowermass planets. This is again possibly the result of planetary atmospheric evolution, as for larger mass planets the atmospheric escape is weaker and the evolution is slower. Even excluding the ${f}_{\text{atm}}^{\text{start}}$ data points characterised by a flat posterior, also the link between ${f}_{\text{atm}}^{\text{start}}$ and M_{⋆} is not significant (R^{2} = 0.068). If statistically significant, such a correlation would go in the direction suggested by Lozovsky et al. (2021) in which atmospheric accretion of hydrogendominated atmospheres might be more efficient for planets orbiting more massive stars (see also Kennedy & Kenyon 2008a,b, 2009). The analysis of additional systems may lead one to constrain this correlation.
Fig. 7. Initial atmospheric mass fractions as a function of planetary semimajor axis (top), planetary mass (middle), and stellar mass (bottom) for the sample of planetary systems considered here. For flat ${f}_{\text{atm}}^{\text{start}}$ PDFs, the uncertainty would span the entire axis and is therefore not displayed for better visualisation; those data points are indicated by empty markers. 
In Fig. 8, we show the possible presence of multidimensional links among the parameters considered above. In particular, the top panel of Fig. 8, which shows how ${f}_{\text{atm}}^{\text{start}}$ varies as a function of a and M_{p}, highlights that PASTA is capable of better constraining ${f}_{\text{atm}}^{\text{start}}$ for highermass planets with a larger orbital separation. In fact, atmospheric loss is controlled by Ṁ_{atm}, which is positively correlated with R_{p}. Because of their low gravitational potential, lowmass planets can be subject of extremely large massloss rates and thus there is a degeneracy between atmospheric evolutionary tracks characterised by large ${f}_{\text{atm}}^{\text{start}}$ and high Ṁ_{atm} or small ${f}_{\text{atm}}^{\text{start}}$ and low Ṁ_{atm}. For highermass planets, instead, the high gravitational potential limits the massloss rates and thus there is less degeneracy among the possible atmospheric evolutionary tracks, which leads to a more defined ${f}_{\text{atm}}^{\text{start}}$. At equal planetary mass, the orbital separation plays the same role as planetary mass in controlling the escape and thus the degeneracy among evolutionary tracks. Therefore, on average, planets orbiting farther away from their host star lead to a more defined ${f}_{\text{atm}}^{\text{start}}$. However, in general, Fig. 8 leads to the same conclusions drawn by Fig. 7, without further particular additions.
Fig. 8. Planetary semimajor axis as a function of planetary mass (top) and stellar mass (middle), and stellar mass as a function of planetary mass (bottom) for the planetary systems considered here. In each panel, the initial atmospheric mass fraction (${f}_{\text{atm}}^{\text{start}}$) is colourcoded and the size of the data points is proportional to the uncertainty (at the σ level) inferred from the ${f}_{\text{atm}}^{\text{start}}$ posterior PDFs. Squares are for flat ${f}_{\text{atm}}^{\text{start}}$ posterior PDFs. 
6. Conclusions
We have characterised the atmospheric evolution of planets composing a dozen exoplanetary systems. In particular, we constrained the initial atmospheric mass fraction (i.e. fraction of planetary atmospheric mass at the dispersal of the protoplanetary nebula) of the planets and the evolution of the rotation rate of the host stars. To this end, we employed a custom algorithm called PASTA, which is a major update to the code described by Kubyshkina et al. (2019a,b). Compared to the previous version of the algorithm, the major upgrade is the refinement of the gyrochronological relation and the treatment of the coefficients entering the empirical relations between stellar EUV emission, L_{X} emission, and rotation period as jump parameters. The excellent priorposterior agreement of all observational, empirically driven, and theoretically driven parameters suggests that the results on the unconstrained parameters are robust, particularly given the heterogeneity of the components making up our framework.
We found that the average PDF of the stellar rotation period at an age of 150 Myr has a median value of ${P}_{\text{rot,150}}=7.{9}_{5.3}^{+9.8}$, which is ∼1σ higher than what was found by Johnstone et al. (2015a) for coeval OC stars. This result, together with the spindown rate distribution at early stages skewed towards low values, may be due to a selection bias. In fact, PASTA requires systems that host at least one closein planet with a hydrogendominated atmosphere. Therefore, the presence of such a planet is favoured around stars that have evolved as slow rotators since the intense XUV emission of fast rotators would have more likely totally removed the planetary primary atmosphere.
We also found that the rotation rate distribution at an age older than 2 Gyr for the planet hosts considered in this work is compatible with that found in the literature, though our distribution tends to be skewed towards lower spindown rates. This is consistent with the conclusions independently drawn by Kovács (2015) and van Saders et al. (2016), who noticed that field stars exhibit lower spindown rates than their cluster counterparts. Similar conclusions were previously reached by Brown (2014) and Maxted et al. (2015) when specifically comparing cluster stars with noncluster planethosting stars. As suggested by Kovács (2015), the differences arising between field stars and OC members may be due to environmental effects (e.g. a different strength of the magnetic field, density of the interstellar medium, or SPI).
We employed PASTA’s results to look for correlations between the planetary initial atmospheric mass fraction and semimajor axis, planetary mass, and stellar mass. We did not find any statistically significant correlation, partially due to the small number of analysed systems but mostly due to the large uncertainties on the values obtained for the initial atmospheric mass fractions. However, the results indicate that PASTA is better at constraining the initial atmospheric mass fraction of highermass planets orbiting farther away from the host star, which is likely due to the slower atmospheric evolution of these planets compared to that of lowermass and closerin planets that are subject to more intense stellar irradiation. Although not statistically significant, our results hint at the possible presence of a negative correlation of the initial atmospheric mass fraction with planetary mass and of a positive correlation with stellar mass. The latter, in particular, would agree with predictions by planetary atmospheric accretion models (Kennedy & Kenyon 2008a,b, 2009; Lozovsky et al. 2021).
This work represents a first step towards extracting information about the initial stages of the evolution of single planetary systems on the basis of the currently observed system parameters, further appropriately accounting for the observational uncertainties. Primarily, the algorithm requires precise and accurate parameters for a large number of systems. This is the goal of a number of ground and spacebased facilities. In particular, the TESS and CHEOPS missions, in conjunction with groundbased highresolution spectrographs, are providing measurements of the required quality for a number of systems that will enable us in the near future to significantly enlarge the sample size. In the future, the PLATO mission will enable us to significantly increase the sample size, still providing the required accuracy on the system parameters. In parallel, we will aim at further improving the models behind PASTA by adding physics (e.g. the inclusion of elements in addition to hydrogen, a selfconsistent calculation of the heating efficiency in the escape models, nonsolar composition stellar evolutionary tracks, and the evolution of the orbital separation) and thus reliability to the results.
Acknowledgments
DK received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 817540, ASTROFLOW). We thank contributors to the Python Programming Language and the free and opensource community, including: MC3 (Cubillos et al. 2017), NUMPY (van der Walt et al. 2011), SCIPY (Jones et al. 2001), MATPLOTLIB (Hunter 2007), CORNER (ForemanMackey 2016).
References
 Alexander, R., Pascucci, I., Andrews, S., Armitage, P., & Cieza, L. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 475 [Google Scholar]
 Barnes, S. A. 2003, ApJ, 586, 464 [Google Scholar]
 Barnes, S. A. 2007, ApJ, 669, 1167 [Google Scholar]
 Barnes, S. A. 2009, in IAU Symp., eds. E. E. Mamajek, D. R. Soderblom, & R. F. G. Wyse, 345 [Google Scholar]
 Barnes, S. A. 2010, ApJ, 722, 222 [Google Scholar]
 Benz, W., Broeg, C., Fortier, A., et al. 2021, Exp. Astron., 51, 109 [Google Scholar]
 Bonfanti, A., Ortolani, S., Piotto, G., & Nascimbeni, V. 2015, A&A, 575, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bonfanti, A., Ortolani, S., & Nascimbeni, V. 2016, A&A, 585, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bonfanti, A., Delrez, L., Hooton, M. J., et al. 2021, A&A, 646, A157 [EDP Sciences] [Google Scholar]
 Borsato, L., Marzari, F., Nascimbeni, V., et al. 2014, A&A, 571, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brown, D. J. A. 2014, MNRAS, 442, 1844 [NASA ADS] [CrossRef] [Google Scholar]
 Buchhave, L. A., Dressing, C. D., Dumusque, X., et al. 2016, AJ, 152, 160 [NASA ADS] [CrossRef] [Google Scholar]
 Carter, J. A., Agol, E., Chaplin, W. J., et al. 2012, Science, 337, 556 [Google Scholar]
 Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102 [Google Scholar]
 Christiansen, J. L., Vanderburg, A., Burt, J., et al. 2017, AJ, 154, 122 [NASA ADS] [CrossRef] [Google Scholar]
 Cochran, W. D., Fabrycky, D. C., Torres, G., et al. 2011, ApJS, 197, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Collier Cameron, A., Davidson, V. A., Hebb, L., et al. 2009, MNRAS, 400, 451 [NASA ADS] [CrossRef] [Google Scholar]
 Cubillos, P., Harrington, J., Loredo, T. J., et al. 2017, AJ, 153, 3 [Google Scholar]
 Davis, T. A., & Wheatley, P. J. 2009, MNRAS, 396, 1012 [CrossRef] [Google Scholar]
 Delrez, L., Ehrenreich, D., Alibert, Y., et al. 2021, Nat. Astron., 5, 775 [NASA ADS] [CrossRef] [Google Scholar]
 Denissenkov, P. A., Pinsonneault, M., Terndrup, D. M., & Newsham, G. 2010, ApJ, 716, 1269 [Google Scholar]
 Dorn, C., Venturini, J., Khan, A., et al. 2017, A&A, 597, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Erkaev, N. V., Kulikov, Y. N., Lammer, H., et al. 2007, A&A, 472, 329 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 ForemanMackey, D. 2016, J. Open Source Softw., 1, 24 [Google Scholar]
 Fossati, L., Erkaev, N. V., Lammer, H., et al. 2017, A&A, 598, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fressin, F., Torres, G., Rowe, J. F., et al. 2012, Nature, 482, 195 [NASA ADS] [CrossRef] [Google Scholar]
 Fulton, B. J., & Petigura, E. A. 2018, AJ, 156, 264 [Google Scholar]
 Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109 [Google Scholar]
 Ginzburg, S., Schlichting, H. E., & Sari, R. 2016, ApJ, 825, 29 [Google Scholar]
 Ginzburg, S., Schlichting, H. E., & Sari, R. 2018, MNRAS, 476, 759 [Google Scholar]
 Gorti, U., Liseau, R., Sándor, Z., & Clarke, C. 2016, Space Sci. Rev., 205, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Gupta, A., & Schlichting, H. E. 2019, MNRAS, 487, 24 [Google Scholar]
 Gupta, A., & Schlichting, H. E. 2020, MNRAS, 493, 792 [Google Scholar]
 Hadden, S., & Lithwick, Y. 2014, ApJ, 787, 80 [Google Scholar]
 Hadden, S., & Lithwick, Y. 2017, AJ, 154, 5 [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Jeans, J. H. 1925, The Dynamical Theory of Gases (Cambridge: Cambridge At The University Press) [Google Scholar]
 Jin, S., & Mordasini, C. 2018, ApJ, 853, 163 [Google Scholar]
 Jin, S., Mordasini, C., Parmentier, V., et al. 2014, ApJ, 795, 65 [Google Scholar]
 Johnstone, C. P., Güdel, M., Brott, I., & Lüftinger, T. 2015a, A&A, 577, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Johnstone, C. P., Güdel, M., Stökl, A., et al. 2015b, ApJ, 815, L12 [Google Scholar]
 Johnstone, C. P., Bartel, M., & Güdel, M. 2021, A&A, 649, A96 [EDP Sciences] [Google Scholar]
 Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open Source Scientific Tools for Python [Google Scholar]
 Kennedy, G. M., & Kenyon, S. J. 2008a, ApJ, 682, 1264 [NASA ADS] [CrossRef] [Google Scholar]
 Kennedy, G. M., & Kenyon, S. J. 2008b, ApJ, 673, 502 [Google Scholar]
 Kennedy, G. M., & Kenyon, S. J. 2009, ApJ, 695, 1210 [NASA ADS] [CrossRef] [Google Scholar]
 Kimura, S. S., Kunitomo, M., & Takahashi, S. Z. 2016, MNRAS, 461, 2257 [NASA ADS] [CrossRef] [Google Scholar]
 Kovács, G. 2015, A&A, 581, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Krenn, A. F., Fossati, L., Kubyshkina, D., & Lammer, H. 2021, A&A, 650, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kubyshkina, D., Fossati, L., Erkaev, N. V., et al. 2018a, ApJ, 866, L18 [NASA ADS] [CrossRef] [Google Scholar]
 Kubyshkina, D., Fossati, L., Erkaev, N. V., et al. 2018b, A&A, 619, A151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kubyshkina, D., Cubillos, P. E., Fossati, L., et al. 2019a, ApJ, 879, 26 [Google Scholar]
 Kubyshkina, D., Fossati, L., Mustill, A. J., et al. 2019b, A&A, 632, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lammer, H., Selsis, F., Ribas, I., et al. 2003, ApJ, 598, L121 [Google Scholar]
 Lissauer, J. J., Fabrycky, D. C., Ford, E. B., et al. 2011, Nature, 470, 53 [Google Scholar]
 Lissauer, J. J., JontofHutter, D., Rowe, J. F., et al. 2013, ApJ, 770, 131 [Google Scholar]
 Lopez, E. D., & Fortney, J. J. 2014, ApJ, 792, 1 [Google Scholar]
 Loyd, R. O. P., Shkolnik, E. L., Schneider, A. C., et al. 2020, ApJ, 890, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Lozovsky, M., Helled, R., Pascucci, I., et al. 2021, A&A, 652, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 MacDonald, M. G. 2019, MNRAS, 487, 5062 [NASA ADS] [CrossRef] [Google Scholar]
 Mamajek, E. E., & Hillenbrand, L. A. 2008, ApJ, 687, 1264 [Google Scholar]
 Marcy, G. W., Isaacson, H., Howard, A. W., et al. 2014, ApJS, 210, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Masuda, K., Hirano, T., Taruya, A., Nagasawa, M., & Suto, Y. 2013, ApJ, 778, 185 [NASA ADS] [CrossRef] [Google Scholar]
 Maxted, P. F. L., Serenelli, A. M., & Southworth, J. 2015, A&A, 577, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mazeh, T., Holczer, T., & Faigler, S. 2016, A&A, 589, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 McDonald, G. D., Kreidberg, L., & Lopez, E. 2019, ApJ, 876, 22 [Google Scholar]
 Meibom, S., Mathieu, R. D., & Stassun, K. G. 2009, ApJ, 695, 679 [Google Scholar]
 Micela, G., Sciortino, S., Serio, S., et al. 1985, ApJ, 292, 172 [CrossRef] [Google Scholar]
 Mills, S. M., Howard, A. W., Weiss, L. M., et al. 2019, AJ, 157, 145 [NASA ADS] [CrossRef] [Google Scholar]
 Montmerle, T., Ehrenreich, D., Lagrange, A. M., & Matsuyama, I. 2010, in Physics and Astrophysics of Planetary Systems, EAS Publ. Ser. (EDP Sciences), 41, 171 [NASA ADS] [Google Scholar]
 Noyes, R. W., Hartmann, L. W., Baliunas, S. L., Duncan, D. K., & Vaughan, A. H. 1984, ApJ, 279, 763 [Google Scholar]
 Otegi, J. F., Bouchy, F., & Helled, R. 2020, A&A, 634, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Owen, J. E., & Lai, D. 2018, MNRAS, 479, 5012 [Google Scholar]
 Owen, J. E., & Wu, Y. 2017, ApJ, 847, 29 [Google Scholar]
 Owen, J. E., Shaikhislamov, I. F., Lammer, H., Fossati, L., & Khodachenko, M. L. 2020, Space Sci. Rev., 216, 129 [Google Scholar]
 Pallavicini, R., Golub, L., Rosner, R., et al. 1981, ApJ, 248, 279 [NASA ADS] [CrossRef] [Google Scholar]
 Palle, E., Nowak, G., Luque, R., et al. 2019, A&A, 623, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Petigura, E. A., Howard, A. W., Lopez, E. D., et al. 2016, ApJ, 818, 36 [Google Scholar]
 Petigura, E. A., Benneke, B., Batygin, K., et al. 2018, AJ, 156, 89 [Google Scholar]
 Pizzolato, N., Maggio, A., Micela, G., Sciortino, S., & Ventura, P. 2003, A&A, 397, 147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rauer, H., Catala, C., Aerts, C., et al. 2014, Exp. Astron., 38, 249 [Google Scholar]
 Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, J. Astron. Telescopes Instrum. Syst., 1 [Google Scholar]
 Rogers, L. A., Bodenheimer, P., Lissauer, J. J., & Seager, S. 2011, ApJ, 738, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Sandoval, A., Contardo, G., & David, T. J. 2021, ApJ, 911, 117 [CrossRef] [Google Scholar]
 SanzForcada, J., Micela, G., Ribas, I., et al. 2011, A&A, 532, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Skumanich, A. 1972, ApJ, 171, 565 [Google Scholar]
 Stökl, A., Dorfi, E., & Lammer, H. 2015, A&A, 576, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sun, L., Ioannidis, P., Gu, S., et al. 2019, A&A, 624, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Szabó, G. M., & Kiss, L. L. 2011, ApJ, 727, L44 [Google Scholar]
 Tu, L., Johnstone, C. P., Güdel, M., & Lammer, H. 2015, A&A, 577, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vanderburg, A., Becker, J. C., Buchhave, L. A., et al. 2017, AJ, 154, 237 [NASA ADS] [CrossRef] [Google Scholar]
 van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22 [Google Scholar]
 Van Eylen, V., Agentoft, C., Lundkvist, M. S., et al. 2018, MNRAS, 479, 4786 [Google Scholar]
 van Saders, J. L., Ceillier, T., Metcalfe, T. S., et al. 2016, Nature, 529, 181 [Google Scholar]
 Vissapragada, S., JontofHutter, D., Shporer, A., et al. 2020, AJ, 159, 108 [NASA ADS] [CrossRef] [Google Scholar]
 Watson, A. J., Donahue, T. M., & Walker, J. C. G. 1981, Icarus, 48, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Weiss, L. M., Marcy, G. W., Rowe, J. F., et al. 2013, ApJ, 768, 14 [Google Scholar]
 Wright, N. J., Drake, J. J., Mamajek, E. E., & Henry, G. W. 2011, ApJ, 743, 48 [Google Scholar]
Appendix A: Plots of the set priors and obtained posteriors for each considered system
After the detailed discussion commenting the results of the K2285 system in Sect. 4, we present here plots showing the outputs of PASTA referring to all the other considered systems. In general, the posteriors of the stellar parameters agree well with their priors for all analysed systems. In a few cases, namely for KOI94 and Kepler48, the posterior and prior for the q_{SF} and m_{SF} coefficients of Eq. (6) are in tension (last row of Figs. A.5 and A.19). However, the posterior estimates are $({\widehat{q}}_{\mathrm{SF}}\u037e{\widehat{m}}_{\mathrm{SF}})=(3.{71}_{1.58}^{+1.57}\u037e0.820\pm 0.057)$ and $(6.{53}_{1.56}^{+1.61}\u037e0.{926}_{0.057}^{+0.058})$ for KOI94 and Kepler48, respectively, which differ less than 1σ from their respective reference values. We remark that, on the one hand, combining together the results coming from several systems leads to a good agreement between the posterior and prior of the empirically derived conversion coefficients, which confirms the global consistency of PASTA’s framework (see discussion in Sect. 5.1). On the other hand, on a starbystar basis ${\widehat{q}}_{\mathrm{SF}}$ and ${\widehat{m}}_{\mathrm{SF}}$ may occasionally differ from their empirical counterpart, which stresses the importance of giving enough degrees of freedom to the framework.
For all systems but Kepler11, Kepler36, and Kepler48, there is always one planet (two in the case of Kepler411 and KOI94) whose presentday atmospheric and core mass impose a strong constraint on the atmospheric evolution so that the ${f}_{\text{atm}}^{\text{start}}$PDF is well defined. As already seen with K2285, the atmospheric modelling had the positive side effect of constraining the masses of some of the other planets (e.g. M_{b, K − 18} of Kepler18, Fig. A.10, leftmost panel of the second row; M_{e, K − 20} and M_{f, K − 20} of Kepler20, Fig. A.12, two rightmost panels of the second row). We theoretically estimated ${\widehat{M}}_{b,\mathrm{K}18}=7.{9}_{1.1}^{+1.2}\phantom{\rule{0.166667em}{0ex}}{M}_{\oplus}$ (reducing the uncertainties by a factor of ∼ 3), and ${\widehat{M}}_{e,\mathrm{K}20}=0.{650}_{0.061}^{+0.062}\phantom{\rule{0.166667em}{0ex}}{M}_{\oplus}$ and ${\widehat{M}}_{f,\mathrm{K}20}=1.{03}_{0.20}^{+0.22}$ (improving from the upper limits given as priors).
As emphasised in Sect. 5.2, the host stars considered in this work were generally slow rotators when they were young, if compared to their OC counterparts, except for Kepler411 whose median ${\overline{P}}_{\mathrm{rot},150}=4.{0}_{2.7}^{+3.5}$ days is slightly lower than the comparison sample, which has ${\overline{P}}_{\mathrm{rotJ},150}=5.{2}_{4.7}^{+3.2}$. As a consequence of being a moderately fast rotator at 150 Myr, when young Kepler411 was likely in the L_{X} saturation regime, hence the posterior exhibits a peak at high L_{X} values (log L_{X} ∼ 29.5, Fig. A.17, leftmost panel of the second row) that is much more pronounced than the corresponding peaks obtained for the other stars.
Finally, Kepler11 deserves a special, separated comment. The star hosts six planets (from b to g) for which there is no general consensus about their masses in the literature. Therefore, we merged different probability distributions according to the results published by several authors to produce the reference priors for the planetary masses (see Sect. 3). For each planet with a nonflat mass prior (b to f), our posterior theoretical estimates ${\widehat{M}}_{p}$ differ less than 1.5σ from our adopted prior modes (Fig. A.8, second row), mostly as a result of the large prior uncertainties especially on M_{b} and M_{c}. We also note that the ${\widehat{M}}_{p}$ values differ less than 1σ from the corresponding values proposed by Hadden & Lithwick (2014) or Lissauer et al. (2011). The loose constraints we imposed on M_{p} do not enable PASTA to derive sharply peaked ${f}_{\text{atm}}^{\text{start}}$PDFs, except partly for ${f}_{\mathrm{atm},g}^{\mathrm{start}}=0.{043}_{0.023}^{+0.016}$ (see Fig. A.8, last row), which is so distant from its host (${\overline{a}}_{g}=0.4660\pm 0.0040$ AU) that it has basically retained its original atmospheric content (see Kubyshkina et al. 2019b, for a discussion about this). The different evolutionary scenarios that may have characterised the Kepler11 system are also reflected by the lack of constraints on the mass of planet g, whose posterior follows the flat prior, thus adding no information. We remark that all the M_{p} values present in the literature for this system are based on transit time variations, which suffer of degeneracies when used to constrain planetary masses and eccentricities (Hadden & Lithwick 2014). As this is a challenging measurement, it is not surprising to see some kind of tension between PASTA’s outputs and the adopted priors, and a more accurate study of the evolution of the Kepler11 system will need to wait until more precise planetary masses become available.
Fig. A.10. Same as Fig. 2, but for the planetary parameters of the Kepler18 system. 
Fig. A.11. Same as Fig. 1, but for the starrelated properties of Kepler20. 
Fig. A.12. Same as Fig. 2, but for the planetary parameters of the Kepler20 system. 
Fig. A.13. Same as Fig. 1, but for the starrelated properties of Kepler25. 
Fig. A.14. Same as Fig. 2, but for the planetary parameters of the Kepler25 system. 
Fig. A.15. Same as Fig. 1, but for the starrelated properties of Kepler36. 
Fig. A.16. Same as Fig. 2, but for the planetary parameters of the Kepler36 system. 
Fig. A.17. Same as Fig. 1, but for the starrelated properties of Kepler411. 
Fig. A.18. Same as Fig. 2, but for the planetary parameters of the Kepler411 system. 
Fig. A.19. Same as Fig. 1, but for the starrelated properties of Kepler48. 
Fig. A.20. Same as Fig. 2, but for the planetary parameters of the Kepler48 system. 
Fig. A.21. Same as Fig. 1, but for the starrelated properties of WASP47. 
Fig. A.22. Same as Fig. 2, but for the planetary parameters of the WASP47 system. 
All Tables
Saturation values (R_{X, sat}) adopted in our work for different mass ranges, as inferred from McDonald et al. (2019).
Planetary radii R_{p}, masses M_{p}, and semimajor axes a taken from the literature and adopted in this work, followed by the main results obtained with PASTA (median values and 1σ confidence interval), that is the stellar rotation period at 150 Myr (P_{rot,150}), the theoreticallyestimated planetary mass (${\widehat{M}}_{\mathrm{p}}$), and the atmospheric mass fraction at t = t_{disk} = 5 Myr (${f}_{\text{atm}}^{\text{start}}$).
All Figures
Fig. 1. Stellar properties of the K2285 system given by PASTA. First row: P_{rot,150}, t_{⋆}, P_{rot, ⋆}, and M_{⋆}. Second row: L_{X} at 150 Myr, and the x and y exponents adopted by the gyrochronological relations given by Eqs. (1) and (2). Third row: q_{SF} and m_{SF} coefficients of Eq. (6), and the β exponent of Eq. (5). The output posterior PDFs are shown as blue histograms, with the green area representing the 68.3% HPD credible interval. The red curves are the Gaussian or flat (uninformative) imposed priors. When absent (second row), no specific prior was imposed on the parameter; in particular, x and L_{X} are mathematically derived from the PDFs of the jump parameters. The black histogram in the topleftmost panel represents the P_{rot,150} distribution extracted from the sample of Johnstone et al. (2015a) considering stars with a mass differing by less than 0.1 M_{⊙} from the mass of K2285. 

In the text 
Fig. 2. PDFs given by PASTA for the considered planetary parameters in the K2285 system: semimajor axis (first row), mass (second row), and ${f}_{\text{atm}}^{\text{start}}$ (third row). The black curves in each panel of the third row (evident in the second panel, and just barely visible in the others as they are squeezed towards zero) indicate the predicted presentday atmospheric mass fraction. It follows that K2285 b, d, and e have almost entirely lost their atmospheres at some point in time along the system’s evolution. 

In the text 
Fig. 3. Corner plot of the relevant jump parameters returned by PASTA for the K2285 system. 

In the text 
Fig. 4. PDFs (blue histograms) of those jump parameters entering as coefficients in Eqs. (5) and (6). Left and middle panels: intercept q_{SF} and slope m_{SF} of Eq. (6). Right panel: β exponent of Eq. (5). The median and the ±1σ bounds of each PDF are marked as vertical solid and dashed red lines, respectively. The black Gaussians represent the priors inferred from observations for each empirical parameter. 

In the text 
Fig. 5. Top: PDF of the x exponent of our adopted gyrochronological relations expressed in Eqs. (1) and (2) built upon all the stellar targets analysed in this work. The vertical solid red line marks the median of the distribution, while the two vertical dashed red lines marks the σ levels corresponding to the 15.87th and 84.14th percentiles. Bottom: same as the top panel, but for the y exponent of Eq. (2). We note that the top panel is in log scale. 

In the text 
Fig. 6. Top: comparison between the P_{rot,150} PDF obtained by combining those distributions extracted from the analysed sample (blue line) and the P_{rot,150} distribution extracted from Johnstone et al. (2015a), considering stars member of ∼150 Myr old OCs and having the same mass distribution as the planethosting stars analysed in this work. The vertical solid red line marks the median value of the distribution shown by the blue line, while the two vertical dashed red lines indicate the 15.87th and 84.14th percentiles. Bottom: same as the top panel, but for the P_{rot,600} PDF. 

In the text 
Fig. 7. Initial atmospheric mass fractions as a function of planetary semimajor axis (top), planetary mass (middle), and stellar mass (bottom) for the sample of planetary systems considered here. For flat ${f}_{\text{atm}}^{\text{start}}$ PDFs, the uncertainty would span the entire axis and is therefore not displayed for better visualisation; those data points are indicated by empty markers. 

In the text 
Fig. 8. Planetary semimajor axis as a function of planetary mass (top) and stellar mass (middle), and stellar mass as a function of planetary mass (bottom) for the planetary systems considered here. In each panel, the initial atmospheric mass fraction (${f}_{\text{atm}}^{\text{start}}$) is colourcoded and the size of the data points is proportional to the uncertainty (at the σ level) inferred from the ${f}_{\text{atm}}^{\text{start}}$ posterior PDFs. Squares are for flat ${f}_{\text{atm}}^{\text{start}}$ posterior PDFs. 

In the text 
Fig. A.1. Same as Fig. 1, but for the starrelated properties of HD 3167. 

In the text 
Fig. A.2. Same as Fig. 2, but for the planetary parameters of the HD 3167 system. 

In the text 
Fig. A.3. Same as Fig. 1, but for the starrelated properties of K224. 

In the text 
Fig. A.4. Same as Fig. 2, but for the planetary parameters of the K224 system. 

In the text 
Fig. A.5. Same as Fig. 1, but for the starrelated properties of KOI94. 

In the text 
Fig. A.6. Same as Fig. 2, but for the planetary parameters of the KOI94 system. 

In the text 
Fig. A.7. Same as Fig. 1, but for the starrelated properties of Kepler11. 

In the text 
Fig. A.8. Same as Fig. 2, but for the planetary parameters of the Kepler11 system. 

In the text 
Fig. A.9. Same as Fig. 1, but for the starrelated properties of Kepler18. 

In the text 
Fig. A.10. Same as Fig. 2, but for the planetary parameters of the Kepler18 system. 

In the text 
Fig. A.11. Same as Fig. 1, but for the starrelated properties of Kepler20. 

In the text 
Fig. A.12. Same as Fig. 2, but for the planetary parameters of the Kepler20 system. 

In the text 
Fig. A.13. Same as Fig. 1, but for the starrelated properties of Kepler25. 

In the text 
Fig. A.14. Same as Fig. 2, but for the planetary parameters of the Kepler25 system. 

In the text 
Fig. A.15. Same as Fig. 1, but for the starrelated properties of Kepler36. 

In the text 
Fig. A.16. Same as Fig. 2, but for the planetary parameters of the Kepler36 system. 

In the text 
Fig. A.17. Same as Fig. 1, but for the starrelated properties of Kepler411. 

In the text 
Fig. A.18. Same as Fig. 2, but for the planetary parameters of the Kepler411 system. 

In the text 
Fig. A.19. Same as Fig. 1, but for the starrelated properties of Kepler48. 

In the text 
Fig. A.20. Same as Fig. 2, but for the planetary parameters of the Kepler48 system. 

In the text 
Fig. A.21. Same as Fig. 1, but for the starrelated properties of WASP47. 

In the text 
Fig. A.22. Same as Fig. 2, but for the planetary parameters of the WASP47 system. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.