Close-in sub-Neptunes reveal the past rotation history of their host stars: atmospheric evolution of planets in the HD3167 and K2-32 planetary systems

TL;DR

This study uses atmospheric escape modeling and Bayesian analysis to infer the early rotation history of host stars from the atmospheric properties of close-in sub-Neptune planets in the HD3167 and K2-32 systems.

Contribution

It introduces a novel Bayesian framework linking planetary atmospheric evolution to stellar rotation history, applied to specific exoplanet systems.

Findings

HD3167's early stellar flux was 40-130 times solar.

K2-32's early stellar flux was 0.5-4 times solar.

The method constrains stellar rotation periods from planetary atmospheres.

Abstract

Planet atmospheric escape induced by high-energy stellar irradiation is a key phenomenon shaping the structure and evolution of planetary atmospheres. Therefore, the present-day properties of a planetary atmosphere are intimately connected with the amount of stellar flux received by a planet during its lifetime, thus with the evolutionary path of its host star. Using a recently developed analytic approximation based on hydrodynamic simulations for atmospheric escape rates, we track within a Bayesian framework the evolution of a planet as a function of stellar flux evolution history, constrained by the measured planetary radius, with the other system parameters as priors. We find that the ideal objects for this type of study are close-in sub-Neptune-like planets, as they are highly affected by atmospheric escape, and yet retain a significant fraction of their primordial hydrogen-dominated atmospheres. Furthermore, we apply this analysis to the HD3167 and K2-32 planetary systems. For HD3167, we find that the most probable irradiation level at 150 Myr was between 40 and 130 times solar, corresponding to a rotation period of 1.78^{+2.69}_{-1.23} days. For K2-32, we find a surprisingly low irradiation level ranging between half and four times solar at 150 Myr. Finally, we show that for multi-planet systems, our framework enables one to constrain poorly known properties of individual planets.

Tables

Table 1. Table 1: Radii and atmospheric mass fractions following a forward-model evolution for the test planets.

	$R_{\mathrm{pl}}$ [$R_{\oplus }$] / $f_{\mathrm{at}}$ [$\%$]
planet	slow	medium	fast
b	$3.44\pm 0.17$ / $2.3\pm 0.8$	$2.80\pm 0.14$ / ${0.2}_{-0.2}^{+0.8}$	$2.71\pm 0.14$ / ${0.0}_{-0.0}^{+0.9}$
c	$5.26\pm 0.26$ / ${11.7}_{-1.0}^{+1.2}$	$4.12\pm 0.21$ / ${6.9}_{-1.3}^{+1.2}$	$3.22\pm 0.16$ / ${3.1}_{-1.3}^{+1.1}$
d	$6.53\pm 0.33$ / ${20.4}_{-1.0}^{+0.3}$	$4.53\pm 0.23$ / ${11.0}_{-0.8}^{+0.4}$	$3.40\pm 0.17$ / ${5.7}_{-0.7}^{+0.5}$

Table 2. Table 2: Adopted system parameters and data, and obtained results ( x 𝑥 x ) for HD3167 and K2-32. The cases for which R pl subscript 𝑅 pl R_{\rm pl}  reaches the core radius are given in parenthesis.

planet	$M_{\ast }$	age	$P_{\mathrm{rot}}$	$M_{\mathrm{pl}}$	$R_{\mathrm{pl}}$	$\rho$	$d_0$	$\Lambda$	$x$
	[$M_{\odot }$]	[Gyr]	[days]	[$M_{\oplus }$]	[$R_{\oplus }$]	[g/cm3]	[AU]
HD3167b	$0.88\pm 0.02$	$7.8\pm 4.3$	$23.52\pm 2.87$	$5.69\pm 0.44$ (a)	$\left(1.57\pm 0.05\right)$ (a)	${8.00}_{-0.98}^{+1.10}$	$0.01752\pm 0.00063$	$15.1$	–
				$5.02\pm 0.38$ (b)	$\left({1.70}_{-0.15}^{+0.18}\right)$ (b)	${5.60}_{-1.43}^{+2.15}$		$12.3$
					$\left(1.56\pm 0.06\right)$ (c)
					$\left({1.62}_{-0.08}^{+0.17}\right)$ (d)
HD3167c	$0.88\pm 0.02$	$7.8\pm 4.3$	$23.52\pm 2.87$	${8.33}_{-1.85}^{+1.79}$ (a)	${2.74}_{-10}^{+0.11}$ (a)	${2.21}_{-0.53}^{+0.56}$	$0.1806\pm 0.0080$	$39.7$	${0.63}_{-0.08}^{+1.14}$
				${9.80}_{-1.24}^{+1.30}$ (b)	${3.01}_{-0.28}^{+0.42}$ (b)	${1.97}_{-0.59}^{+0.94}$		$42.6$	${1.06}_{-0.26}^{+0.86}$
					${2.85}_{-0.15}^{+0.24}$ (c)
					${2.83}_{-0.14}^{+0.31}$ (d)
HD3167d	$0.88\pm 0.02$	$7.8\pm 4.3$	$23.52\pm 2.87$	$6.90\pm 0.71$ (b)	${\left(1.90\pm 0.07\right)}^{\ast }$	–	$0.07757\pm 0.00027$	–	–
K2-32b	$0.86\pm 0.03$	$7.90\pm 4.50$	${61.10}_{-2.53}^{+3.18}$	$16.50\pm 2.70$ (e)	$5.13\pm 0.28$ (e)	$0.67\pm 0.16$	$0.08036\pm 0.00088$	$29.1$	${0.09}_{-0.07}^{+0.34}$
				$21.10\pm 5.90$ (f)	$5.13\pm 0.28$ (f)	$0.67\pm 0.16$		$37.2$
					$5.62\pm 0.45$ (g)
					${5.17}_{-0.16}^{+0.20}$ (d)
K2-32c	$0.86\pm 0.03$	$7.90\pm 4.50$	${61.10}_{-2.53}^{+3.18}$	(e)	$3.01\pm 0.25$ (e)		$0.1399\pm 0.0015$		–
				(f)	${3.48}_{-0.42}^{+0.97}$ (f)
					$3.32\pm 0.33$ (g)
					${3.12}_{-0.12}^{+0.18}$ (d)
K2-32d	$0.86\pm 0.03$	$7.90\pm 4.50$	${61.10}_{-2.53}^{+3.18}$	$10.30\pm 4.70$ (e)	$3.43\pm 0.35$ (e)	${1.38}_{-0.67}^{+0.92}$	$0.1862\pm 0.0020$	$41.6$	–
				(f)	$3.75\pm 0.40$ (f)
					$3.77\pm 0.41$ (g)
					${3.41}_{-0.14}^{+0.26}$ (d)

Equations

$$P_{\mathrm{rot}}(\tau )=0.407[(B-V{)}_0-0.495{]}^{0.325}{\tau }^{0.566},$$

$$P_{\mathrm{rot}}=\{\begin{array}{ll}{P}_{\mathrm{rot}}^{\mathrm{now}}{\left(\frac{\tau }{T_{\mathrm{age}}}\right)}^{0.566},& \text{τ≥2 Gyr}\\ P_{\mathrm{rot}}^{\mathrm{now}}{\left(\frac{2\mathrm{Gyr}}{T_{\mathrm{age}}[\mathrm{Gyr}]}\right)}^{0.566}{\left(\frac{\tau [\mathrm{Gyr}]}{2\mathrm{Gyr}}\right)}^x,& \text{τ<2 Gyr}\end{array}$$

$$\Lambda =\frac{G{M}_{\mathrm{p}}{m}_{\mathrm{H}}}{k_{\mathrm{b}}{T}_{\mathrm{eq}}{R}_{\mathrm{p}}},$$

Full text

planets and satellites: atmospheres — planets and satellites: physical evolution — planets and satellites: gaseous planets — planets and satellites: individual (HD3167b, HD3167c, HD3167d, K2-32b, K2-32c, K2-32d)

††journal: ApJ

1 Introduction

For late-type stars (i.e., later than mid F-type), rotation rate and high-energy radiation (X-ray and EUV below 912 Å; hereafter XUV) are intimately connected, with faster rotating stars being XUV brighter (e.g., Pallavicini et al., 1981; Pizzolato et al., 2003; Johnstone et al., 2015b). During their main-sequence lifetime, the rotation rate and XUV flux of late-type stars decrease with time. This evolution does not follow a unique path, as stars born with the same mass and metallicity can have widely different initial rotation rates, thus XUV fluxes (e.g., Mamajek & Hillenbrand, 2008; Johnstone et al., 2015b). This non-uniqueness holds up to about 2 Gyr, at which point the evolutionary tracks of both rotation rate and XUV emission converge to one path (e.g., Tu et al., 2015). This implies that, for stars older than about 2 Gyr, it is not possible to infer their past XUV emission from their present stellar properties.

This problem becomes relevant when trying to understand the evolution of planetary atmospheres. The XUV stellar flux is absorbed in the upper atmosphere of planets leading to atmospheric expansion and eventually escape, which has been shown to occur efficiently on close-in Jupiter- and Neptune-mass planets (e.g., Vidal-Madjar et al., 2003; Fossati et al., 2010; Lecavelier des Etangs et al., 2012; Ehrenreich et al., 2015). Furthermore, atmospheric escape plays a significant role in shaping the observed exoplanet population and atmospheric properties (e.g., Lopez & Fortney, 2013, 2014; Jin et al., 2014; Tian & Ida, 2015; Owen & Wu, 2017; Owen & Lai, 2018; Jin & Mordasini, 2018). Since in most cases atmospheric mass-loss rates are directly proportional to their XUV irradiation, escape is most effective when planets are young, thus the importance of knowing the past evolution of the stellar XUV radiation.

Luckily, for stars hosting close-in planets with hydrogen-dominated atmospheres, we can infer the stellar XUV irradiation tracks by modeling the planets’ atmospheric evolution. The amount of XUV flux received by a planet during its lifetime determines the amount of hydrogen remaining in its envelope: since the radius of planets with masses below 0.35 Jupiter masses strongly correlates with its hydrogen envelope content (e.g., Lopez & Fortney, 2014; Hatzes & Rauer, 2015), the currently observed planetary radius constrains the past evolution of the stellar XUV emission.

Here, we present a framework, based on mass-loss rates extracted from hydrodynamic simulations (Kubyshkina et al., 2018b), enabling the modeling of the atmospheric evolution of super-Earths and sub-Neptunes that we apply to the HD3167 and K2-32 planetary systems. Section 2 summarises the algorithms and tools employed to compute the planetary atmospheric evolutionary tracks. Section 3 shows the results of tests conducted on mock systems. Section 4 presents the results obtained for the HD3167 and K2-32 planetary systems, while Section 5 gives our discussion and conclusions.

2 Atmospheric Evolution Model

To model the evolution of planetary atmospheres, we employ the scheme described by Kubyshkina et al. (2018b). Starting from a set of initial system parameters, as a planet loses mass from its hydrogen envelope, at each time-step the code extracts the atmospheric mass-loss rate from an analytical approximation based on a grid of one dimensional hydrodynamic models, uses this value to update the atmospheric mass fractions, and estimates the planetary radius, which we then compare to the observed radius once the age of the system has been reached. Here, we summarise the key steps and describe the applied improvements.

There are three key ingredients: a model of the stellar flux evolutionary track, a model relating planetary parameters and atmospheric mass, and a model computing atmospheric escape rates.

For late-type stars, the stellar XUV flux out of the saturation regime depends on the stellar mass and rotation period, where the latter is time-dependent. Mamajek & Hillenbrand (2008) suggested the following analytical approximation for the evolution of the stellar rotation period (in days)

[TABLE]

where $\tau$ is the stellar age (in Myr), and $(B-V)_{0}$ is the reddening-free stellar color. Equation (1) represents the average approximation based on a large set of late-type dwarfs. However, the evolutionary tracks for stellar rotation rates are non-unique (Johnstone et al., 2015b; Tu et al., 2015). To account for the different rotation-rate histories of different stars (i.e., from slow to fast rotator), we model the rotation period as a power law in $\tau$, normalised such that the rotation period at the present age ($T_{\rm age}$) is consistent with the measured stellar rotation period ($P_{\rm rot}^{\rm now}$), obtaining

[TABLE]

where the exponent $x$ is a positive value, typically ranging between 0 and $\sim$2, controlling the stellar rotation period at ages younger than 2 Gyr. Figure 1 shows examples of the parametric $P_{\rm rot}$ curves for a range of $x$ values equal to 0, 1, and 2.

We then derive the stellar X-ray and EUV luminosity from the rotation period following Pizzolato et al. (2003) and Wright et al. (2011) that relate rotation rates and stellar masses to X-ray luminosity, including saturation effects for the fastest rotators, and Sanz-Forcada et al. (2011) relating X-ray and EUV fluxes. To account for variations of the stellar bolometric luminosity with time, we use the MESA Isochrones and Stellar Tracks (MIST, Paxton et al., 2018) model grid.

To estimate the planetary mass-loss rate at a given time, we use the analytic formulas provided by Kubyshkina et al. (2018c), which are based on a grid of one-dimensional hydrodynamic upper-atmosphere models (Kubyshkina et al., 2018b). The analytic approximation requires as input stellar XUV flux, planetary mass ($M_{\rm pl}$), radius ($R_{\rm pl}$), equilibrium temperature ($T_{\rm eq}$), and orbital separation ($d_{0}$).

To estimate the atmospheric mass of a planet $M_{\rm atm}$ as a function of planetary mass, radius, and $T_{\rm eq}$, we pre-compute $M_{\rm atm}$ following Johnstone et al. (2015) for the range of parameters used in Kubyshkina et al. (2018c), among which we further interpolate during an evolution run. For a given core mass, we determine the core radius $R_{\rm core}$ assuming an Earth-like density, as for planets with hydrogen-dominated envelopes the core composition has little influence on planetary size (e.g., Lopez & Fortney, 2014; Petigura et al., 2016).

We combine these ingredients to compute the planetary atmospheric evolutionary tracks. For each run, we assume that the planetary orbital separation and stellar mass remain constant, and ignore the contribution of gravitational contraction and radioactive decay on the planetary $T_{\rm eq}$ during the first phases of evolution. We further assume that every planet accreted a hydrogen-dominated atmosphere from the proto-planetary nebula. We begin our simulations at 5 Myr, which is the typical lifetime of protoplanetary disks (Mamajek, 2009). We also consider an initial planetary radius corresponding to a value of the restricted Jeans escape parameter $\Lambda$ (computed for the planetary core mass and $T_{\rm eq}$ at the beginning of the simulation) of 5, where (Jeans, 1925; Fossati et al., 2017)

[TABLE]

with $G$ gravitational constant, $m_{\rm H}$ mass of an hydrogen atom, and $k_{\rm b}$ Boltzmann constant. We come back to this point in Section 3.

At each step of the evolution, first we compute the mass-loss rate based on the stellar flux and system parameters, which we use to update the atmospheric mass fraction, and the planetary radius. We adjust the time step, such that the change in atmospheric mass is less than 5% of $M_{\rm atm}$.

Finally, we use a Bayesian approach to constrain the evolutionary track of the stellar XUV luminosity fitting the currently observed planetary radius. To this end, we combine the planetary evolution model with the open-source Markov-chain Monte Carlo (MCMC) algorithm of Cubillos et al. (2017), to compute the posterior distribution for the stellar rotation rate and the considered system parameters.

For each MCMC run, we let the planetary mass, age of the system, present-time rotation period, orbital separation, and stellar mass as free parameters, with Gaussian-like priors according to the measured values and uncertainties. Within our scheme, we adopt a stellar radius from the MESA evolutionary tracks, which depends on the stellar mass. Similarly, we compute $T_{\rm eq}$ derived from the adopted stellar radius and $d_{0}$, assuming zero Bond albedo and full energy redistribution. The stellar radius may also be implicitly included in the computation in case only an upper limit for $P_{\rm rot}^{\rm now}$ has been given. We parameterise the stellar rotation evolution using $x$ of Equation 2, adopting a uniform prior ranging between 0 and 2. However, we provide the results in terms of the stellar rotation period at the age of 150 Myr, because this is more intuitive than $x$ and it allows one to compare the results with the distribution of stellar rotation periods measured in young open clusters.

3 Validation

3.1 Initial conditions

Kubyshkina et al. (2018b) showed that for planets less massive than 6 $M_{\oplus}$ and orbital separations smaller than 0.1 AU, the initial radius has no impact on the final radius at the end of the evolution. This is because, for a given $M_{\rm pl}$ and $d_{0}$, a larger (smaller) radius, hence larger (smaller) $M_{\rm atm}$, leads to a higher (lower) mass-loss rate. Here, we look more thoroughly for the $M_{\rm pl}$ and $d_{0}$ values above which the atmospheric evolution drawn by our framework starts to become sensitive to the initial radius.

The main parameters controlling atmospheric escape at a given XUV irradiation and $R_{\rm pl}$ are $M_{\rm pl}$ and $d_{0}$. Therefore, we consider systems composed of a 0.8 $M_{\odot}$ medium rotator ($x=1$) and planets with masses ranging between 1 and 20 $M_{\oplus}$ and orbital separations ranging between 0.02 and 0.50 AU. For each $M_{\rm pl}$–$d_{0}$ pair, we computed a series of evolutionary tracks over a range of initial planetary radii, looking for the time at which the evolutionary tracks merge ($T_{\rm merge}$), namely the atmospheric mass fractions become the same within 2%, and the minimum initial radius for which this happens.

We set the upper limit of the tested initial radii at 10 $R_{\oplus}$, the maximum $R_{\rm pl}$ available in the grid of upper atmosphere models (Kubyshkina et al., 2018b). We set the lower limit of the tested initial radii based on the minimum expected atmospheric mass fraction $f_{\rm at}=M_{\rm atm}/M_{\rm p}$. Since within the protoplanetary disc the mass of the initially accreted envelope grows with increasing planetary mass, we expect $f_{\rm at}>10^{-4}$ for $M_{\rm pl}$ $\leq$ 5 $M_{\oplus}$ and $f_{\rm at}=10^{-2}$ for higher masses (Stökl et al., 2015).

Figure 2 shows that $T_{\rm merge}$ increases with $M_{\rm pl}$ and $d_{0}$, and that it remains always below $\sim$2 Gyr. Figure 2 further shows the minimum atmospheric mass fraction $f_{\rm at,min}$ required for the evolutionary tracks to merge as a function of $M_{\rm pl}$ and $d_{0}$. Both $T_{\rm merge}$ and $f_{\rm at,min}$ increase with increasing $M_{\rm pl}$ and $d_{0}$. For systems with $M_{\rm pl}$ above 5 $M_{\oplus}$ and $d_{0}$ larger than 0.2 AU, $f_{\rm at,min}$ lies above the atmospheric mass fraction that can be expected to be originally accreted by each planet before the dispersal of the protoplanetary disc based on Stökl et al. (2015). Therefore, the results of the atmospheric evolution drawn by our framework start to become sensitive to the initial radius for planets more massive than 5 $M_{\oplus}$ and lying at orbital separations larger than 0.2 AU.

When applied to real cases, the results presented in Figure 2 should be considered as indicative. This is because the mass and rotation rate, thus XUV emission, of a given host star are possibly different from the values considered here. However, we can conclude that for sub-Neptune-like planets within 0.2 AU the initial conditions are in general not important, while for heavier and/or more distant planets this is not the case. We, however, evaluate in each case the impact of the choice of the initial conditions on the results.

3.2 Reproducibility of mock systems

To test whether and in which cases the retrieval approach can constrain the stellar XUV luminosity evolutionary tracks, we carry out an injection-retrieval test. To do so, we create a 5 Gyr old system, with a stellar mass ($M_{*}$) of 0.8 $M_{\odot}$ and $P_{\rm rot}^{\rm now}=30$ days. We simulate three planets (labeled “b”, “c”, and “d”) with masses of 20, 15, and 8 $M_{\oplus}$ at orbital separations of 0.05, 0.1, and 0.2 AU, respectively. These parameters make sure that the initial planetary radius, set as the value for which $\Lambda=5$, has no influence on the results. The core radii of the three planets are 2.71, 2.47, and 2.00 $R_{\oplus}$, respectively.

3.2.1 Typical uncertainties

To estimate the system-parameter uncertainties, we look at the typical measured uncertainties listed in the NASA Exoplanet Archive111http://exoplanetarchive.ipac.caltech.edu. For planets similar to those considered in this test (i.e., $M_{\rm pl}$ below 25 $M_{\oplus}$, excluding upper limits, $T_{\rm eq}$ between 300 and 2000 K, stellar mass between 0.4 and 1.3 $M_{\odot}$), we adopt stellar mass uncertainties of 4% (signal-to-noise ratio, S/N, of 25), orbital separation uncertainties of 1.5% (S/N = 66.7), planetary mass uncertainties of 20% (S/N = 5), and stellar age uncertainties of 30% (S/N = 3.3). Note that the estimated age uncertainties vary significantly; however, this uncertainty has the largest influence only for systems much younger than the focus of this work (see below).

Finally, for the measured rotation period it is not possible to reliably associate a typical uncertainty because the number of stars for which $P_{\rm rot}^{\rm now}$ has been measured is small, and in most cases $P_{\rm rot}^{\rm now}$ either corresponds just to an upper limit or has been estimated from approximations, such as that of Mamajek & Hillenbrand (2008). On the basis of the few stars for which $P_{\rm rot}^{\rm now}$ has been measured, we adopt an uncertainty of 15% (S/N $\approx$ 6.7).

3.2.2 Results

We simulate three evolution scenarios with $x$ = 0.1, 0.566, and 1.5, corresponding to an initially “slow”, “medium”, and “fast” rotator or alternatively rotational periods at 150 Myr ($P_{\rm rot}^{150}$) of 13.8, 4.1, and 0.4 days, respectively. The planetary radii and atmospheric mass fractions (accounting for the uncertainties on the parameters) obtained for these scenarios and for planets “b”, “c”, and “d” following a forward model are listed in Table 1. When orbiting a medium or a fast rotator, planet “b” loses (almost) completely its atmosphere as, within the uncertainties, its atmospheric mass fraction remains below 1.0%.

For all measured/derived system parameters we set a gaussian-like prior defined by the injected values and their uncertainties. For the parameter $x$, which defines the evolution of the stellar rotation period, we set a uniform prior between 0 and 2.

As a typical example of the results obtained for the test systems, Figure 3 shows the MCMC posterior distributions for the case of the “c” planet orbiting the medium rotator. The results obtained for the other tests are shown in the Appendix (Figures 10–17).

In Figure 3, the top row shows the posterior distribution obtained when all system parameter were set free. To identify whether there are dominant parameters driving the results and which ones they are, we repeated the analysis fixing one parameter at a time. For completeness, we run a test fixing all parameters except for the stellar rotation period at 150 Myr, recovering the injected $P_{\rm rot}$ distribution.

For each modelled case, the posterior distribution for the stellar rotation period matches the injected value within the 68% highest posterior density (HPD) credible interval. Our tests indicate that planetary mass is the main parameter influencing the posterior distribution of the rotation period, with the other parameters playing a minor role.

Our tests also show that the posterior distribution of all parameters, except for the planetary mass, are consistent with the input distribution. For the planetary mass, we find that the posterior distributions are consistently shifted towards higher values than their prior distributions; however, the results are consistent within 1$\sigma$ with the prior. The pair-wise posterior distributions (Figure 4) reveal a well-defined degeneracy between $P_{\rm rot}^{150}(M_{\rm pl})$ and $M_{\rm pl}$, which occurs because one can arrive at the same planetary radius through a lower (higher) planetary mass and a slower (faster) rotator, thus lower (higher) past XUV irradiation. This explains the shift between prior and posterior distributions. This effect is more evident for closer-in planets and for slower rotators. In the case of the slow-rotator simulation, we also find a slight shift of the posterior of the stellar rotation period compared to the injected value, but still lying well within 1$\sigma$. This indicates that the algorithm prefers to reach the final planetary radius with a slightly heavier planet and a slightly faster rotator.

For the planets orbiting the fast rotator, the posteriors for planetary mass are narrower than the prior. This happens also for the test planet “b”, which is the closer-in of the three, orbiting the medium rotator (Figure 13), for which the final atmospheric mass fraction is below 1%. This happens because the thin atmosphere remaining at the end of the simulation can be obtained only through a narrow range of planetary masses. This indicates that our model imposes constraints on the plausible physical parameters of the system, leading to constraints on the planetary masses.

The left and middle panels of Figure 4 show the pair-wise posterior distributions of $P_{\rm rot}^{150}$ vs. $M_{\rm pl}$ and $T_{\rm age}$ vs. $M_{\rm pl}$, for the medium-rotator case when all parameters are free. The correlation is rather narrow, meaning that, despite the whole range of possible stellar rotation rates is wide, at each given planetary mass the range of possible stellar rotation rates is narrow. For the case of the slow rotator, we also find a correlation between $M_{\rm pl}$ and $T_{\rm age}$. This occurs because the atmospheric mass fraction changes more slowly over time for the slow rotator. A correlation between $P_{\rm rot}^{150}$ and $T_{\rm age}$ is the direct consequence of the other two correlations (right panel of Figure 4). The other system parameters do not show significant correlations.

For some of the test planets, such as planets “b” and “c” orbiting the fast rotator, we obtained broad posterior distributions, thus weak constraints on the stellar rotation rate history. The case of planet “b”, the closer-in and most massive planet, orbiting the fast or medium rotator lead to the (almost) complete escape of the atmosphere. Therefore, even though formally the injected value of the rotation period lies within the 68% HPD credible interval, we can just obtain an upper limit for the stellar rotation period. Our tests indicate that for planets with a thin hydrogen envelope, a meaningful use of our scheme requires a high accuracy of the system parameters: the thinner the hydrogen atmosphere, the higher the required accuracy on the system parameters.

To make sure that our results do not depend on the adopted masses and positions of the test planets, we carry out additional runs considering test planets orbiting at distances from the host star ranging between 0.1 and 0.25 AU and having masses ranging between 5 and 20 $M_{\oplus}$. The results of these tests, not shown here, are consistent with those presented above, confirming that the accuracy with which we are able to constrain the stellar rotation period depends mainly on the accuracy of the others parameters (predefined outside of our framework) and the uncertainty of the data.

3.2.3 Influence of the uncertainties

From the test runs, illustrated in Figures 3 and 10–17, we expect that the main uncertainties affecting our ability to constrain the stellar rotation history are those on $M_{\rm pl}$ and $R_{\rm pl}$. Therefore, we performed additional tests to examine how the results change by changing the uncertainties on these values. For the planetary radius, we varied the uncertainty between 2 and 50%, while for the planetary mass, we varied the uncertainty between 5 and 50%. In addition, because of the large spread with which ages of planetary systems are derived, we present the results of additional runs considering age uncertainties of 50 and 100%.

Figure 5 shows the results of these tests. The top-left panel is the reference distribution, which coincides with the top-left panel of Figure 3. The first, second, and third columns present the results obtained varying the uncertainties on the system’s age, planetary radius, and planetary mass, respectively. There is little variation in the resulting posterior distributions for $P_{\rm rot}^{150}$ when varying the system’s age uncertainty, because the main changes in atmospheric mass fraction occur during the early stages of evolution—before 2 Gyr. Therefore, the uncertainty on the system’s age plays a role only for very young systems. As expected, increasing the uncertainty on the planetary radius and mass broadens the posterior distribution of $P_{\rm rot}^{150}$, becoming nearly flat for a radius and/or mass uncertainties of 50% (not shown). Improving the planetary-mass uncertainty has the largest impact on the expected $P_{\rm rot}^{150}$ posterior distribution. For the real planets examined in this study, the uncertainties on planetary mass are 13–16%.

3.2.4 Simultaneous modelling of multiple planets in the same

system

Here we explore what can be obtained by analysing simultaneously two or more planets orbiting the same star. Figure 6 shows the posterior distributions for $P_{\rm rot}^{150}$ obtained by simultaneously modelling the “b”, “c”, and “d” test planets orbiting a slow, medium, and fast rotator. For comparison, Figure 6 further shows the distributions obtained analysing the planets separately. As one can expect, considering simultaneously multiple planets narrows, and thus improves, the results. This can be best seen in the middle panel (i.e., medium rotator), where the innermost planet constrains the upper boundary of $P_{\rm rot}^{150}$, while the outer planets constrain the lower boundary of $P_{\rm rot}^{150}$, and in the right panel, for which the final distribution is significantly narrower.

The mock system we used here has been constructed in such a way that the uncertainties on the planetary parameters are homogeneous among the three planets. However, this is unlikely to be the case for real systems. Therefore, we expect that a meaningful simultaneous analysis of multi-planet systems requires that the masses and radii of the planets considered in the analysis have sufficiently low uncertainties.

The results obtained on multi-planet systems are not different from those of single-planet systems. This is not surprising given that those results correspond to the superposition of what is separately obtained for each single planet in the system. This indicates that for multi-planet systems, having one planet with a good mass estimate is enough to constrain the stellar rotation history, even if modeling all planets in a system simultaneously.

3.2.5 Constraining poorly known planetary masses in multi-planet systems

Above, we have shown that for systems with well estimated planetary masses a simultaneous analysis of all planets may better constrain the stellar rotation history. Here we consider the case in which one or more planets in a system have a poorly constrained or unknown mass. To consider this case, we slightly modified the MCMC algorithm as follows. We turned the poorly constrained or unknown planetary mass into the retrieval parameter and set the history of the stellar rotation period (more precisely $x$) as prior, where its distribution has been derived employing the planet with the well constrained mass.

Figure 7 shows as an example the results obtained to recover the mass of the test planet “c” lying within the mock system with a medium rotator, employing the value of the exponent $x$ derived from the analysis of planet “d”, and otherwise, the results obtained for the mass of the planet “d”, assuming the exponent $x$ resolved for the planet “c”. For both test planets, we considered a uniform mass prior on the planetary mass ranging between 1 and 39 Earth masses. Figure 7 shows that the true masses lie within the obtained 68$\%$ HPD credible intervals.

4 Application to Observed Planetary Systems

We now apply our atmospheric evolution framework to the HD3167 and K2-32 planetary systems aiming at estimating the rotation rate evolution of their host stars. Table 2 lists the adopted system parameters and results. Although in Section 3.2.4 we conclude that it is possible to derive the stellar rotation history of the host star with a simultaneous analysis of all planets in the system, here we consider just one planet in each system. This is because for each of the HD3167 and K2-32 planetary systems just one planet holds an hydrogen atmosphere and/or has a mass which has been measured to better than 30$\%$. For the K2-32 system, we use the $P_{rot}$ distribution obtained using the planet with a well-measured mass to constrain the rotation history of the host star and the mass of the other planets in the system (see Section 3.2.5).

4.1 HD 3167

The HD3167 planetary system is composed of three sub-Neptune-mass planets, where the innermost (HD3167b) and outermost (HD3167c) planets transit their host star (Vanderburg et al., 2016; Gandolfi et al., 2017; Christiansen et al., 2017). All planets in the system have masses derived from radial-velocity measurements.

HD3167b has a measured average density consistent with a rocky composition, suggesting that the planet has lost its primordial hydrogen-dominated atmosphere (Gandolfi et al., 2017). We indeed find that the hydrogen envelope of HD3167b escapes completely within the first 0.01 Myr for any reasonable $x$ value (i.e., $0<x<2$). We cannot use this planet for our analysis as it would return a flat distribution of the stellar rotation period.

HD3167c has an average density of $\sim$2 g cm*-3*, suggesting that the planet still retains part of its primary hydrogen-dominated atmosphere. The evolutionary tracks do not depend on the initial conditions, except for a slow rotator and a planetary mass above 10.9 $M_{\oplus}$, for which the atmospheric evolutionary tracks merge at $f_{\rm at}$ $\approx$ 7$\times$10*-2* and at an age of $\approx$1 Gyr. Figure 8 shows the posterior distributions obtained adopting the system parameters given by Christiansen et al. (2017), indicating that the host star was likely to be a fast to moderate rotator. In particular, we obtain that at 150 Myr the star emitted 38–128 times the XUV flux of the current Sun, where the upper boundary corresponds to the saturation limit. The pair-wise posterior distributions of the parameters involved in the modeling are shown in Appendix (Figure 18). We run the simulations again adopting the system parameters of Gandolfi et al. (2017) obtaining consistent results, except for a slightly broader $P_{\rm rot}^{150}$ posterior distribution towards larger values (i.e., slower rotator) due to the lower planetary mass. To demonstrate the constraining power of our technique, Figure 8 shows the results obtained from the analysis of HD3167c in comparison with the distribution of stellar rotation periods at an age of 150 Myr measured for 0.8 and 0.9 solar mass stars members of open clusters (Johnstone et al., 2015b).

HD3167d is not transiting and therefore does not have a measured radius. Nevertheless, the stellar rotation rate evolutionary tracks based on the results obtained from HD3167c lead us to constrain the possible atmospheric evolution of HD3167d suggesting that the planet should have completely lost its primary, hydrogen-dominated atmosphere within 0.02–1.48 Gyr.

4.2 K2-32

The K2-32 planetary system is composed of three transiting sub-Neptune-mass planets. Radial-velocity measurements placed loose constraints on the planetary masses (Dai et al., 2016; Petigura et al., 2017): K2-32c has an upper-mass limit, whereas for K2-32d the mass is measured with a 50% uncertainty. The masses of these planets are too poorly constrained to be useful for our analysis. Therefore, we focused our analysis on K2-32b, which has a low average density of 0.67$\pm$0.16 g cm*-3*, strongly indicative of the presence of an extended hydrogen envelope. For the analysis, we consider the system parameters given by Petigura et al. (2017).

Our results (Figure 9) strongly favor the slow rotator scenario with the posterior distribution indicating $P_{\rm rot}^{150}$ values lying between 10 and 26 days, which coincides with the slow rotation edge of the distribution given by Johnstone et al. (2015b), and corresponds to XUV fluxes at 150 Myr between 0.5 and 4.0 times those of the present Sun. This means that, within our framework, the existence of such a close-in, low-density planet is possible only if the star has always been a very slow rotator. The pair-wise posterior distributions of the parameters involved in the modeling are shown in the Appendix (Figure 19).

Given the weak mass constraints for K2-32c and K2-32d, we used the results on the stellar rotation rate history obtained from the analysis of K2-32b to place further constraints on the estimated planetary masses, as described in Section 3.2.5. The low measured average density of K2-32d (1.38${}^{+0.92}_{-0.67}$ g cm*-3*) suggests that the planet has a significant hydrogen-dominated envelope. Assuming the range of possible stellar rotation rates found from the analysis of K2-32b, we expect the mass of K2-32d to range between 5.60 and 6.87 $M_{\oplus}$.

For K2-32c, by adopting a mass equal to the upper mass limit and the results on the stellar rotation rate derived from K2-32b, we obtain a planetary radius significantly larger than the measured one, with the planet retaining nearly 50% of its initial envelope. To fit the observed $R_{\rm pl}$, we require planetary masses ranging between 4.55 and 7.29 $M_{\oplus}$. By adopting instead the planetary radius given by Mayo et al. (2018), the range of possible planetary masses shifts slightly towards larger values, namely 4.76–7.52 $M_{\oplus}$.

We further test our results by carrying out an additional MCMC run including all three planets. With this, we aim at confirming that simultaneously considering all three planets does not affect the results on the stellar rotation history, and test if the posterior distributions for the masses of K2-32c and K2-32d are consistent with the above estimates. The results of this test are shown in Figure 20. As expected, the $P_{\rm rot}^{150}$ distribution is identical to that obtained considering only K2-32b. Furthermore, the 68$\%$ HPD credible intervals for the masses of planets c and d lie between 6.16 and 7.68 $M_{\oplus}$ and between 5.36 and 6.73 $M_{\oplus}$, respectively. Except for a slight difference in the lower mass boundary for planet c, these values are comparable to those obtained using the algorithm described in Section 3.2.5.

5 Discussions and conclusions

The currently observed atmospheric properties of exoplanets strongly depend on their mass-loss history, which can reveal the otherwise unknown evolution history of the stellar high-energy radiation. We developed a framework to compute planetary radii as a function of time and coupled it to a Bayesian sampler to infer the rotation rate history of the host star, which is constrained by the measured basic system parameters. We applied this framework to the HD3167 and K2-32 systems.

Despite having very similar masses, we find that the K-type stars HD3167 and K2-32 had a widely different initial rotation-rate, thus high-energy emission. Most likely, HD3167 evolved as a medium (or fast) rotator, while K2-32 evolved as an extremely slow rotator. We further showed that our evolution model can lay further constraints for planetary systems in which some of the planetary parameters are poorly measured. For the K2-32 system, the stellar rotation rate constrained by K2-32b, leads to masses ranging between 6.16 and 7.68 $M_{\oplus}$ for K2-32c, and masses ranging between 5.36 and 6.73 $M_{\oplus}$ for K2-32d.

We found that in comparison to single-planet systems, multi-planet systems have a major advantage. Multi-planet systems containing more than one planet with well measured mass enable to more precisely constraint the stellar rotation history. Multi-planet systems containing one planet with well-measured mass and other planets with poorly measured masses our technique enables one to use the former to derive the stellar rotation history, which can then be used to constrain the mass of the latter planets.

We have shown that atmospheric evolution modeling of close-in planets has the power to uncover the past rotation history of late-type stars. It enables one also to constrain poorly known parameters (i.e., mass and/or radius) of some of the planets in a system. These results can be extended to other systems provided that they host a close-in planet with a mass close to that of Neptune or smaller and with a hydrogen-dominated envelope. A systematic analysis of a large sample of systems will also help to set constraints on the different model assumptions adopted by this framework. Furthermore, by comparing the output stellar rotation period distribution at a given age with what is measured in clusters, or by comparing the results obtained from two planets in a system, it would be possible to identify planets that evolved differently from what expected. These planets might have been affected, for example, by strong impacts, or experienced an early dispersal of the protoplanetary nebula, or underwent strong hydrogen release from the solidifying rocky core (e.g., Chachan & Stevenson, 2018; Bonomo et al., 2019). Firmly identifying these “odd” systems would be the first step towards studying them to better understand the evolution of planetary systems.

Our results show also the importance of precisely constraining system parameters, particularly the system’s age and planetary mass and radius. The ESA PLATO mission (Rauer et al., 2014) will be in the best position to simultaneously discover transiting planetary systems and provide system parameters with the necessary accuracy to pin down the evolutionary track of the planetary atmosphere and stellar rotation rate.

We acknowledge the FFG project P853993, the FWF/NFN projects S11607-N16 and S11604-N16, and the FWF projects P27256-N27 and P30949-N36. N.V.E. acknowledges support by the Russian Science Foundation grant No 18-12-00080. We thank the anonymous referee for insightful comments.