© ESO 2021

1 Introduction

The great variety of planetary systems discovered in the last two decades, with populations of gaseous, rocky planets orbiting different type of stars (e.g. Perryman 2011; Howard et al. 2012), has made it necessary to address several open questions about the origin and evolution of these systems as a whole. In this context, a natural synergy has been established between the asteroseismology and exoplanetology communities (e.g. Christensen-Dalsgaard et al. 2010; Batalha et al. 2011; Huber et al. 2013a,b; Huber 2018; Campante et al. 2018), providing a highly accurate characterisation of the host star, which leads to a more accurate determination of the planetary parameters (e.g. Reese et al. 2012; Buldgen et al. 2015). Developing more refined and sophisticated stellar models, in particular regarding the rotational properties and the exchange of angular momentum (AM) between the host star and the planetary orbits, is crucial to better constrain the evolutionary history of the systems.

In this framework, we aim to revisit the Kepler-444 system. We have provided a detailed seismic characterisation in Buldgen et al. (2019, Paper I hereafter), leading to new determinations for the mass, radius, and age of the host star (0.754 ± 0.030 M_⊙, 0.753 ± 0.010 R_⊙,11.0 ± 0.8 Gyr), and hence new estimates for the mass of the Kepler-444-d and e planets. These values were found to be in good agreement with previous determinations obtained by Campante et al. (2015). Here we would like to take a step further, studying the orbital evolution of Kepler-444-d and Kepler-444-e and the eventual impact of atmospheric evaporation.

Kepler-444 is one of the oldest planetary system known thus far. This is a rather peculiar system, composed of five sub-Earth-sized planets orbiting the primary star Kepler-444-A, within 0.1 AU and P_orb < 10 days, along with two spatially unresolved M dwarfs (Kepler-444-BC) at a projected distance of 66 AU (Lillo-Box et al. 2014; Dupuy et al. 2016). The highly eccentric stellar orbit (e = 0.864 ± 0.023) brings the M dwarf binary companion very close to the planetary system every ~460 yr (Campante et al. 2015; Mills & Fabrycky 2017), within $5_{-1.0}^{+0.9}$ AU. This likely primordial orbital configuration is expected to have strongly impacted the protoplanetary disc, truncating it at ≈2 AU and severely depleting it of the solid material from which the total ≈1.5 M_⊕ planetary mass would have formed (Dupuy et al. 2016). This particular configuration makes the early history of the system puzzling since the formation and migration of planets in such a truncated disc are not well understood (Papaloizou 2016; Mills & Fabrycky 2017). Some studies suggest that such tightly packed planetary systems would have formed locally; however, the estimated densities for the two outermost planets Kepler-444-d and e, being in between a Mercury-like rock/iron ratio and water world composition, would accredit the hypothesis of a formation beyond the snowline (Mills & Fabrycky 2017; Buldgen et al. 2019).

Bourrier et al. (2017) detected significant variations in the Kepler-444-A Ly-α signal at different observational epochs. They observed flux variations of the order of 20 and 40% during the transiting and non-transiting periods, respectively. These variations could be related to stellar activity, explaining the signal detection in and out of transit. However, Bourrier et al. (2017) argue that this scenario is unlikely given the quietness of this old main-sequence star. Instead, they propose that the ~20% variation could be due to the presence of escaping hydrogen-rich atmospheres from the two outermost planets (Kepler-444-e and Kepler-444-f) and the ~40% variation could result from the presence of an undetected Kepler-444-g planet as of yet at a larger orbital distance. If the planetary origin of the signal is confirmed, Bourrier et al. (2017) propose that Kepler-444-e and f could be examples of ocean planets. They would have formed beyond the snowline and migrated to the current orbital distance, retaining a substantial amount of water ice to replenish the escaping atmosphere along the ~11 Gyr of evolution (Bourrier et al. 2017).

In this study, we use the comprehensive asteroseismic characterisation of Kepler-444-A performed in Paper I to provide a detailed modelling of this planetary system, which accounts for the rotational evolution of the host star together with the impact of tides and X-ray and extreme ultraviolet (EUV) fluxes on the planets. We use the Geneva stellar evolution code (Eggenberger et al. 2008) that follows the rotational evolution of the star by computing the internal angular momentum (AM) transport through shear instability, meridional circulation, and magnetic instabilities (Eggenberger et al. 2010a,b, 2019). These stellar models¹ are then coupled to the orbital evolution code (Privitera et al. 2016a,b; Meynet et al. 2017; Rao et al. 2018) to study the impact of dynamical and equilibrium tides on the architecture of the system. We also study the impact of atmospheric evaporation on Kepler-444-e, using the prescription for the hydrodynamical escape in the energy-limited regime (Lecavelier Des Etangs 2007; Erkaev et al. 2007). We investigate the impact of different rotational histories, XUV-luminosity prescriptions, and heating efficiencies on the estimation of the initial planetary mass.

2 Method and implemented physics

In order to characterise the evolution of a planetary system, providing robust and precise stellar parameters for the host star through asteroseismic modelling is crucial. Accurate parameters of the host star (in particular, the mass, radius, and age) are indeed needed to better understand the global evolution of the planetary system. In the present study, the stellar parameters determined in Paper I were used as key inputs and constraints to compute rotating stellar models of Kepler-444. The computation of the evolution of the host star was then coupled to the computation of the evolution of the orbits of the planets to obtain a global view of the evolution of the system. In what follows, we first briefly recall the physics implemented in the stellar and orbital evolution codes.

2.1 Geneva stellar evolution code

Stellar models are computed by using a version of the Geneva stellar evolution code (GENEC) specifically optimised for the evolution of low- and intermediate-mass stars. The physics implemented in the code is described in detail in Eggenberger et al. (2008); here, we only briefly recall the fundamental properties that are of interest for the present study. Rotational effects are included with the hypothesis of shellular rotation (Zahn 1992). The internal AM transport was followed simultaneously with the evolution of the star by accounting for meridional currents, transport by the shear instability, and transport by magnetic fields. The transport by magnetic fields was taken into account in the context of the Tayler-Spruit dynamo (Spruit 2002). While correctly accounting for the approximately uniform rotation in the radiative interior of the Sun, we recall that the exact description of this process as well as its viability under realistic stellar conditions is still a matter of debate (Braithwaite 2006; Zahn et al. 2007; Braithwaite & Spruit 2017; Fuller et al. 2019). This of course illustrates that the general modelling of the impact of magnetic fields on internal stellar properties remains an open question. The equation describing the transport of AM is then given by the following: $$\rho \frac{\text{d}}{\text{d}t}{\left(r^2\Omega \right)}_{M_r}=\frac{1}{5r^2}\frac{\partial }{\partial r}\left(\rho r^4\Omega U(r)\right)+\frac{1}{r^2}\frac{\partial }{\partial r}\left(\rho \left(D_{\text{shear}}+{\nu }_{\text{TS}}\right)r^4\frac{\partial \Omega }{\partial r}\right),$$(1)

with r being the radius, and ρ(r) and Ω(r) being the mean density and angular rotation on an isobar, respectively. The velocity U(r) corresponds to the radial component of the meridional circulation velocity in the radial direction and D_shear is the diffusion coefficient for AM transport by the shear instability (Eggenberger et al. 2010c). The viscosity term ν_TS was included to account for the AM transport by magnetic fields, which only operates when the radial differential rotation is larger than a given threshold (Spruit 2002). Atomic diffusion was included and diffusion coefficients were computed following the prescription by Paquette et al. (1986). A solar calibrated value of the mixing-length parameter was used for convection. In addition to internal AM transports, the braking of the stellar surface rotation due to magnetised winds was taken into account by using the braking law of Matt et al. (2015, 2019), which describes the global angular momentum loss of Sun-like stars, defined in this case as stars with masses ranging between about 0.5 and 1.3 M_⊙, with outer convective envelopes that are magnetically active (Matt et al. 2015, 2019). The equations for the torque read as follows: $$\frac{\text{d}J}{\text{d}t}=\{\begin{array}{ll}-T_{\odot }{\left(\frac{R}{R_{\odot }}\right)}^{3.1}{\left(\frac{M}{M_{\odot }}\right)}^{0.5}{\left(\frac{{\tau }_{\text{cz}}}{{\tau }_{\text{cz}\odot }}\right)}^{\text{p}}{\left(\frac{\Omega }{{\Omega }_{\odot }}\right)}^{\text{p+1}},\hfill & if \left(\text{Ro}>{\text{Ro}}_{\odot }/\chi \right),\hfill \\ -T_{\odot }{\left(\frac{R}{R_{\odot }}\right)}^{3.1}{\left(\frac{M}{M_{\odot }}\right)}^{0.5}{\chi }^p\left(\frac{\Omega }{{\Omega }_{\odot }}\right),\hfill & \text{if} \left(\text{Ro}\le {\text{Ro}}_{\odot }/\chi \right).\hfill \end{array}$$(2)

We note that R_⊙ and M_⊙ are the radius and the mass of the Sun, R and M are the radius and the mass of the stellar model, τ_cz is the convective turnover timescale, and Ro is the Rossby number, defined as the ratio between the stellar rotational period and the convective turnover timescale (Ro = P_⋆∕τ_cz). The quantity χ ≡ Ro_⊙∕Ro_sat indicates the critical rotation rate for stars with a given ${\tau }_{\text{cz}}/{\tau }_{{\text{cz}}_{\odot }}$, defining thetransition from the saturated to unsaturated regime. In this work, χ is considered equal to 10 or 13 depending on the relation adopted between the Rossby number and the predicted X-ray luminosity (see Sect. 2.2). The exponent p is considered equal to 2.3 and the constantT_⊙ was calibrated in order to reproduce the solar surface rotation rate (Eggenberger et al. 2019). We recall that, once this calibration is performed, this braking law is expected to correctly predict the global angular momentum loss of stars with different masses (Matt et al. 2015). This is due to the fact that the effects of a change in mass on the physical properties of convective envelopes are consistently taken into account in this expression through the computation of the convective turnover timescales. In addition to the internal transport of AM and to the torque due to magnetised winds, the exchanges of AM between the star and the planetary orbits due to equilibrium and dynamical tides are taken into account as described in Sect. 2.2.

Solar models computed with the physics described above have been proven to be in very good agreement with helioseismic data in reproducing the approximately flat rotation profile in the solar radiative zone, thanks to the contribution of magnetic instabilities in the transport of AM (Eggenberger et al. 2005, 2019). Moreover, the same models are also able to correctly reproduce the distribution of surface rotation rates observed for solar-like stars in open clusters (Eggenberger et al. 2019). The capability of these models to reproduce surface and internal rotational properties of solar-like stars well is of key importance for modelling the rotational evolution of a planet host star and to study its impact on the global evolution of an exoplanetary system.

2.2 Orbital evolution code

The evolution of exoplanetary systems was computed using the orbital evolution code of Privitera et al. (2016a,b), Meynet et al. (2017) and Rao et al. (2018), considering the case of a planet on a circular, coplanar orbit around its host star. A full model of the host star was computed as explained above, giving as initial inputs the stellar parameters obtained from the asteroseismic characterisation of the system, going from the pre-main sequence (PMS) phase to the end of the main sequence (MS) or up to the red giant branch (RGB) tip, depending on the physical case we are interested in. Stellar models were then coupled with the orbital evolution code in which we gave the initial orbital distance, the mass of the planet, and the initial angular velocity of the host star as input. In our simulations, we did not account for orbital migration mechanisms occurring across the protoplanetary discs; we did indeed start our computations once the protoplanetary disc had completely disappeared and followed the eventual, subsequent migration of the planets related to the impact of tidal dissipation in the host star convective envelope. We recall that in our study, we do not consider possible mutual interactions among planets. The physics of the code is described in Rao et al. (2018). In the following, we briefly recall the fundamental equations. The total change of the orbital distance is as follows: $$\left(\stackrel{.}{a}/a\right)=-\frac{{\dot{M}}_{\star }+{\dot{M}}_{\text{pl}}}{M_{\star }+M_{\text{pl}}}-\frac{2}{M_{\text{pl}}{v}_{\text{pl}}}\left[F_{\text{fri}}+F_{\text{gra}}\right]+{\left(\dot{a}/a\right)}_{\text{t}},$$(3)

where Ṁ_⋆ = −Ṁ_loss is the stellar mass loss rate, M_pl and Ṁ_pl are the planetary mass and the rate of change of the planetary mass, and v_pl is the planetary orbital velocity. We note that F_fri and F_gra are the frictional and gravitational drag forces, respectively, whose expressions are the same as in Villaver & Livio (2009), Mustill & Villaver (2012) and Villaver et al. (2014). The effect of dynamical and equilibrium tides is included in the term ${\left(\dot{a}/a\right)}_t$. The expression for the equilibrium tides (Zahn 1966, 1977, 1989; Alexander et al. 1976; Livio & Soker 1984; Villaver & Livio 2009; Mustill & Villaver 2012; Villaver et al. 2014) is the same as in Privitera et al. (2016b). The effect of equilibrium tides is only taken into account when an external convective envelope is present. The expression for equilibrium tides is as follows: $${\left(\dot{a}/a\right)}_{\text{eq}}=\frac{f}{{\tau }_{\text{cz}}}\frac{M_{\text{env}}}{M_{\star }}q(1+q){\left(\frac{R_{\text{star}}}{a}\right)}^8\left(\frac{{\Omega }_{\star }}{{\omega }_{\text{pl}}}-1\right),$$(4)

where f is a numerical factor obtained from integrating the viscous dissipation of the tidal energy across the convective zone (Villaver & Livio 2009), M_env is the mass of the convective envelope, q is the ratio between the mass of the planet and the one of the star (q = M_pl∕M_⋆), Ω_⋆ is the angular velocity at the stellar surface, and ω_pl = 2π∕P_orb is the orbital angular velocity of the planet; τ_cz is the convective turnover timescale already introduced in the preceding section.

In addition to equilibrium tides, we also accounted for the impact of dynamical tides, namely for the frequency-average tidal dissipation of inertial waves excited in the convective envelope of the rotating star, whose restoring force is the Coriolis force (Bolmont & Mathis 2016; Gallet et al. 2017; Bolmont et al. 2017; Benbakoura et al. 2019). Their effect is accounted for when the condition ω_pl < 2 Ω_⋆ is satisfied since a planetary companion in a circular-coplanar orbit around a uniformly rotating host star is able to excite inertial waves when the orbital frequency is lower than twice the stellar angular rotation (Ogilvie & Lin 2007). We use the following expression for dynamical tides given by Ogilvie (2013) and Mathis (2015): $${\left(\stackrel{.}{a}/a\right)}_{\text{dyn}}=\left(\frac{9}{2Q_{\text{d}}^{\prime }}\right)q{\omega }_{\text{pl}}{\left(\frac{R_{\star }}{a}\right)}^5\frac{({\Omega }_{\star }-{\omega }_{\text{pl}})}{\mid {\Omega }_{\star }-{\omega }_{\text{pl}}\mid },$$(5)

with $Q_{\text{d}}^{\prime }=3/(2D_{\omega })$ and $D_{\omega }=D_{0\omega }{D}_{1\omega }{D}_{2\omega }^{-2}$. The ‘D’ terms are as follows: $$\{\begin{array}{l}{D}_{0\omega }=\frac{100\pi }{63}{ϵ}^2\frac{{\alpha }^5}{1-{\alpha }^5}{(1-\gamma )}^2,\hfill \\ D_{1\omega }={(1-\alpha )}^4{\left(1+2\alpha +3{\alpha }^2+\frac{3}{2}{\alpha }^3\right)}^2,\hfill \\ D_{2\omega }=1+\frac{3}{2}\gamma +\frac{5}{2\gamma }\left(1+\frac{\gamma }{2}-\frac{3{\gamma }^2}{2}\right){\alpha }^3-\frac{9}{4}\left(1-\gamma \right){\alpha }^5,\hfill \end{array}$$(6)

where α = R_c∕R_⋆, β = M_c∕M_⋆, $\gamma =\frac{{\alpha }^3(1-\beta )}{\beta (1-{\alpha }^3)}$, and $ϵ=\frac{{\Omega }_{\star }}{\sqrt{\frac{{\text{GM}}_{\star }}{{\text{R}}_{\star }^3}}}$. We note that M_c and R_c are the mass and the radius of the radiative core, considered as the region of the star below the base of the convective envelope. The term D_ω was computed in Ogilvie (2013) as the frequency-averaged tidal dissipation. This dissipation is referred to as the kinetic energy of the wave-like solution, which was obtained from the dynamical equations system for the stellar envelope (Ogilvie 2013; Mathis 2015).

Tides become inefficient when the planetary distance is equal to the corotation radius, namely when the orbital frequency of the planet is equal to the stellar surface angular rotation (ω_pl = Ω_⋆). The effect that tides have is to widen the orbit when the planet is beyond the corotation radius and to shrink it when the planet is inside.

We also consider the evolution and impact of the stellar XUV flux on the exoplanetary atmosphere. Following the work by Tu et al. (2015), we first used the formula for the computation of the XUV flux of Wright et al. (2011). This prescription gives the ratio of the X-ray luminosity with respect to the bolometric luminosity of the star (R_X = L_X∕L_⋆) in saturated and unsaturated regimes as a function of the Rossby number. The value of the Rossby number establishes the transition from one regime to the other (see Sect. 2.1): $$\{\begin{array}{ll}{R}_{\text{X}}={10}^{-3.13}\times {\left(\frac{\text{Ro}}{{\text{Ro}}_{\text{sat}}}\right)}^{-\beta },\hfill & \text{if}\text{Ro}>{\text{Ro}}_{\odot }/\chi \hfill \\ R_X={10}^{-3.13},\hfill & \text{if}\text{Ro}\le {\text{Ro}}_{\odot }/\chi .\hfill \end{array}$$(7)

A value of χ = 10 was used with this expression for the transition from a saturated to unsaturated regime. In the saturated regime, a constant value R_X = 10^−3.13 was adopted as given by Wright et al. (2011). In the unsaturated regime, a solar-calibrated value was used for the exponent β with a value of −3.25. We note that this value slightly differs from the value of −2.7 reported by Wright et al. (2011), which reflects differences in the computation of the Rossby numbers. In the present case, the Rossby numbers were directly computed from the properties of our rotating models. As in Tu et al. (2015), we also used the prescription of Sanz-Forcada et al. (2011) to obtain the corresponding EUV luminosity: $${\log }_{10}(L_{\text{EUV}})=4.8+0.86{\log }_{10}(L_{\text{X}}).$$(8)

The expressions given by Eqs. (7) and (8) were used by default in the present study to compute the X-ray and EUV luminosities from the structural and rotational properties of our models. However, more recent studies suggest different prescriptions for these relations, in particular regarding the hypothesis of a constant value of the R_X ratio in the saturated regime (e.g. Reiners et al. 2014; Johnstone et al. 2021). Following the work by Johnstone et al. (2021), we then also included a more recent prescription defining the evolution of the X-ray luminosity as a function of the Rossby number: $$\{\begin{array}{ll}{R}_{\text{X}}={\text{C}}_1 {\text{Ro}}^{{\beta }_1},\hfill & \text{if}\text{Ro}>{\text{Ro}}_{\odot }/\chi \hfill \\ R_{\text{X}}={\text{C}}_2 {\text{Ro}}^{{\beta }_2},\hfill & \text{if}\text{Ro}\le {\text{Ro}}_{\odot }/\chi .\hfill \end{array}$$(9)

The main difference between this expression and Eq. (7) is related to the fact that a dependence on the Rossby number is also taken into account in the saturated regime. As for the case of the expression of Wright et al. (2011), the values of the different coefficients entering Eq. (9) cannot be directly taken from the work of Johnstone et al. (2021) due to differences in the way Rossby numbers are computed. Based on the work of Johnstone et al. (2021), we then determined the constants and exponents of the broken power law of Eq. (9) by doing a linear fit of the R_X −Ro distribution of stars taken from the catalogue of Wright et al. (2011). A value of χ = 13 is then obtained for the transition between the saturated and unsaturated regime with the following values for the coefficients: C₁ = 6.27 × 10⁻⁶, C₂ = 3.71 × 10⁻⁴, β₁ = −2.275, and β₂ = −0.177. The non-zero value of β₂ confirms that a small dependence on the rotation rate is favoured over a constant value of R_X in the saturated regime, which is in good agreement with the results of Johnstone et al. (2021). Moreover, the value determined for β₁ indicates that a steeper slope is obtained in the unsaturated regime when the relation is calibrated to reproduce the solar values as done previously in Eq. (7). As a consequence, using Eq. (9) instead of Eq. (7) leads to higher values for the X-ray luminosity for old main-sequence solar-type stars (see discussion in Sect. 4.3 for the impact of these different prescriptions in the specific case of Kepler-444).

In addition to suggesting a new R_X – Ro relation, Johnstone et al. (2021) also reconsidered the relation between the EUV and X-ray fluxes. Consequently, when Eq. (9) is adopted for the computation of the X-ray luminosity, the corresponding equations for the EUV fluxes for the hard (10–36 nm) and soft (36–92 nm) components, respectively, are also the same as in Johnstone et al. (2021) and read as: $$\begin{array}{ll}{\log }_{10}(F_{\text{EUV},\text{1}})\hfill & =2.04+0.681{\log }_{10}(F_{\text{X}}),\hfill \\ {\log }_{10}(F_{\text{EUV},\text{2}})\hfill & =-0.341+0.920{\log }_{10}(F_{\text{EUV},\text{1}}).\hfill \end{array}$$(10)

Fig. 1 Evolutionary tracks in the HR diagram obtained with the GENEC (in blue) and CLES (in salmon) evolution codes together with the location of the optimal asteroseismic solution (red dot). The tracks show the evolution from the PMS to part of the MS phase. The GENEC track in particular corresponds to the moderate rotation case.

3 Rotational and orbital evolution

3.1 Rotational history of the system

Rotating models of Kepler-444 have been computed with the GENEC code using the input physics described in the previous section. The initial parameters of these models are taken from Paper I in order to reproduce the optimal asteroseismic solution. The rotational history of Kepler-444-A being unknown, we computed models representative of an initially slow, moderate, and fast rotating star by using different values for the initial rotation velocity and disc lifetimes during the PMS. The choice of these values was determined from surface rotation rates observed for solar-type stars in open clusters of various ages as was done in Eggenberger et al. (2019). Models corresponding to slow, medium and fast rotators were then computed with an initial rotation rate Ω_ini of 3.2, 5, and 18 Ω_⊙, respectively. Figure 1 shows the evolution in the HR diagram of the moderate rotating case (blue line) together with the evolutionary track corresponding to the non-rotating model computed with the CLES code in Paper I (salmon line). The location in the HR diagram of the optimal asteroseismic solution is denoted by a red dot.

In the top panel of Fig. 2, we show the evolution of the surface rotation as a function of the age for the slow, moderate, and fast rotating cases. After the disc-locking phase, during which the surface angular velocity remains constant, the star experiences a spin up due to the rapid contraction on the PMS phase. The surface rotation reaches a peak at approximately 48 Myr, then it starts to decrease due to the braking of the surface induced by magnetised stellar winds. Independently from the initial value of Ω_ini, all the surface rotation tracks then overlap after ~10³ Myr, converging and reproducing the surface rotation of the star well at the age of the system: Ω_K444 = 1.47 × 10⁻⁶ Hz (P_rot = 49.40 ± 6.04 d, Mazeh et al. 2015). We thus find that rotating models able to correctly reproduce the rotational properties of the Sun and the evolution of surface rotation rates of solar-type stars in open clusters predict a surface rotation rate for Kepler-444-A that is in perfect agreement with the observed value. It is important to note that this result has been obtained without any adjustment to the values of the parameters describing the efficiency of surface magnetic braking or internal AM transport. This implies that the observed rotation rate of Kepler-444-A is perfectly in line with the rotational properties expected for such an old K0-type star. We thus find that Kepler-444-A does not exhibit a peculiarly fast rotation compared to similar old stars, which suggests that there is no reason for this star to be especially active. This is an important result regarding the physical nature of the Ly-α variations observed for this star.

The degeneracy in the initial rotation of the star due to the braking of the surface by magnetised winds makes tracing back the rotational history of Kepler-444-A difficult which, therefore, can only be better constrained by looking at the indirect effects of rotation on the planetary system evolution. For example, the impact of dynamical and/or equilibrium tides on the architecture of the system is strongly related to the rotational properties of the host star, as well as the emission of high-energy luminosity and its effect on atmospheric escape. In the bottom panel of Fig. 2, the evolution of the EUV (solid line) and X-ray luminosity (dotted line) is shown, relative to the three initial surface rotations considered. These X-ray and EUV luminosities were computed according to Eqs. (7) and (8) (unless explicitly mentioned, the X-ray and EUV luminosities considered in the present study correspond to the ones predicted by these equations). During the early stages of the evolution, all rotators are in the saturated regime with emitted luminosities following the same trend. The higher the Ω_ini value is, the longer the duration of the saturated regime. Therefore, the slow rotating case is the first one to switch to unsaturated conditions at ~ 125 Myr, followed by the medium and the fast ones. In the case of a fast initial rotation, the planets of the system would be exposed to high-energy irradiation for a rather prolonged period due to the longer duration of the saturated regime.

Fig. 2 Top panel: surface angular velocity evolution for Kepler-444 in the case of a fast rotator (Ωini = 18 × Ω⊙, solid blue line), a medium rotator (Ωini = 5 × Ω⊙, solid magenta line), and a slow rotator (Ωini = 3.2 × Ω⊙, solid black line). The red dot represents the surface angular velocity of Kepler-444-A at the age of the system (11 Gyr). The red-dashed line shows the critical velocity limit. The vertical black-dotted line indicates the age of the system, while the horizontal red-dotted line shows the value of the surface rotation at the age of the system. Bottom panel: evolution of the EUV (solid lines) and X-ray luminosities (dashed lines) computed according to Eqs. (7) and (8) relative to the different rotational histories considered. The vertical black-dotted line indicates the age of the system.

3.2 Orbital evolution: Kepler-444-d and Kepler-444-e

We studied the impact of dynamical tides on the total change of the orbital distance for the two planets for which an estimate of the mass is available, namely Kepler-444-d and Kepler-444-e. We first chose the current measured values of the orbital distances as initial input for our computations (see Table 1) and considered the masses to remain constant along the entire evolution. For the planetary masses, we took the value at the errorbar upper limit (M_d = 0.1016 and M_e = 0.0921 M_⊕), while for the host star we considered the case of a fast rotator (Ω_ini = 18 × Ω_⊙) in order to maximise the impact of dynamical tides. As a result of the simulation, we did not find any significant change in the orbital evolution of the system related to dynamical tides (see Fig. 3). This result is due to the low mass ratio between planets and the host star (q = m_pl∕M_⋆), which substantially weakens the impact of dynamical tides. Moreover, the initial orbital distances of such low mass planets are too large for dynamical tides to have an appreciable impact on the evolution of the architecture of the system.

To investigate these two points in more detail, additional computations were first performed for the evolution of a planet with the same initial distance as Kepler-444-d (i.e. 0.06 AU), but for the following different values of the planetary mass: 10, 100, 317, and 1000 M_⊕. The results of these computations are shown in Fig. 4. Even in the case of a planet of 10 M_⊕, no change in the orbital distance is seen. When the mass is increased to 100 M_⊕, a slight increase in the orbital distance is obtained as a result of the more efficient dynamical tides. A significant impact of dynamical tides on the evolution of the semi-major axis is only seen for more massive planets and, in particular, for the 1000 M_⊕ planet. This explains the evolution at a constant orbital distance obtained for the light planets of the Kepler-444 system seen in Fig. 3.

In the same way, the impact of the initial orbital distance on the efficiency of tides is studied by computing models with a mass fixed to the maximal value allowed for Kepler-444-d (0.1016 M_⊕). Figure 5 shows the orbital evolution corresponding to different values of the initial semi-major axis. We find that only for values of the initial orbital distance lower than about 0.012 AU, a significant change is observed for the orbital evolution. This is again a much lower value than the current orbital distance of the planets in the Kepler-444 system. This illustrates the fact that dynamical tides can only have an impact on very close-in and massive planets, which is not the case for the Kepler-444 system.

Fig. 3 Top panel: orbital evolution of Kepler-444-d and e (solid red and green line, respectively), starting with the current value of the orbital distance ($a_{{\text{ini}}_{\text{d}}}=0.06$, $a_{{\text{ini}}_{\text{e}}}=0.0696$), at a fixed mass (Md = 0.1016 M⊕, Me = 0.0921 M⊕) and host star initial angular rotation Ωin = 18 × Ω⊙. The vertical red-dotted line corresponds to the age of the system (age = 11 Gyr). The blue-dashed lines indicate the evolution of the corotation radii, the blue-dotted ones indicate the evolution of the minimum initial orbital distances below which dynamical tides are no longer active. The magenta area represents the extension of the stellar convective envelope, while the cyan area represents the stellar radiative core. Bottom panel: orbital angular velocity of the planets and the double of the host star angular velocity (blue-dashed line) as a function of time. Dynamical tides are active as long as the orbital frequency of the planets is lower than double the host star surface angular rotation.

Fig. 4 Impact of a change in the planetary mass on the orbital evolution for a planet located at the same distance as Kepler-444-d (initial value fixed to $a_{{\text{ini}}_{\text{d}}}=0.06$ AU). The four red lines span a range of masses from 10 to 1000 M⊕, with larger masses corresponding to more significant increases in the distances. The dashed lines, dotted lines, and shaded regions are the same as in Fig. 3 (top panel).

Fig. 5 Impact of a change in the initial orbital distance for a planet with the maximal value of the mass allowed for Kepler-444-d (mass of the planet fixed to 0.1016 M⊕). The dashed lines, dotted lines, and shaded regions are the same as in Fig. 3 (top panel).

3.3 Orbital evolution: Kepler-444-b

We finally investigated the specific case of the innermost planet of the system (Kepler-444-b). As discussed above, the impact of dynamical tides can be more significant due to the closer orbit of Kepler-444-d to the host star. Although there is no mass measurement for Kepler-444-b, it is possible to have an estimate of this quantity using the radius-mass scaling law obtained by Valencia et al. (2006) for super-Mercury planets ($R_{\text{pl}}\propto M_{\text{pl}}^{0.3}$). The mass of Kepler-444-b computed by means of this relation is about m_b = 0.048 M_⊕. As shown in Fig. 5, no change in the orbit is expected for such a light planet at the initial orbital distance of Kepler-444-b (about 0.04 AU). This is confirmed by the specific computation of the orbital evolution of this planet using the method described in the previous section (as in the case of Kepler-444-d and Kepler-444-e), according to which the impact of dynamical tides on the planet orbit is also negligible in the case of Kepler-444-b, even considering an initial fast rotation for the host star.

4 Impact of evaporation

In the previous sections, we studied the orbital evolution of Kepler-444-d, Kepler-444-e, and Kepler-444-b using the current measurements for the orbital distances and planetary masses as initial inputs and testing the impact of different initial rotations of the host star. We did not find any appreciable changes in the planetary orbits. This result is mainly due to the low values of the planetary masses, or more precisely to the low values of the mass ratio q = M_pl∕M_⋆ and the relatively large values for the initial orbital distances. In the preceding studies, we however made the simplifying assumption of a constant planetary mass along the whole evolution.

If we want to have a more comprehensive study of the system, we should include the impact of atmospheric evaporation. This process maybe extremely relevant in changing the fate of the system, especially in the present case of low-mass close-in planets (a ≤ 0.1AU) that are strongly exposed to high-energy stellar radiation, eventually undergoing atmospheric hydrodynamical escape (e.g. Ehrenreich et al. 2015; Jin & Mordasini 2018). The atmospheric escape may notably influence the evolution of the planetary radius and mass, with distinct efficiencies depending on multiple parameters related to both the planet properties and the rotational evolution of the host star. According to Jin & Mordasini (2018), low-mass low-density planets are most likely to entirely lose their envelopes on relatively short timescales (2 Myr) due to their low gravitational binding energy. In case the planets formed outside of the ice line, they could have retained a significant water ice content while migrating towards the host star. This ice, which is heated up by the stellar radiation, can replenish the thin planetary atmosphere with H during the long evolution of the system. This would allow one to explain why, despite its old age, there would still be a phenomenon of hydrogen escape from these low-mass planets.

Using the Space Telescope Imaging Spectrograph of the Hubble Space Telescope, Bourrier et al. (2017) observed variations in the HI Lyman-α line (1215.6701 Å) of Kepler-444-A. They detected significant variations at different observational epochs. These variations are up to the ~ 20% during the transit of Kepler-444-e and Kepler-444-f. When no transit occurs, the variation is even stronger, up to the ~ 40% (Bourrier et al. 2017). The nature of these variations, which were observed during in and out of transit epochs of the two planets, could be attributed to stellar variability. However, their amplitude over short timescales (~ 2−3 h) seems to be too intense for such an old main sequence K0-type star (Bourrier et al. 2017). Another interesting mechanism eligible to explain these variations is the presence of neutral H-rich exospheres trailing the Kepler-444 planets. According to this hypothesis, the authors attribute the 20% variation to the transit of Kepler-444-e and Kepler-444-f, while the 40% variation is attributed to the partial transit of an undetected outermost planet as of yet, Kepler-444-g. They performed simulations with the EVaporating Exoplanets (EVE) code to interpret the Ly-α variations in terms of exospheric absorption, finding that these observations could be explained by neutral hydrogen mass loss rates of the order of 10⁷−10⁸ g s⁻¹ for Kepler-444-d, e, and f. More observations of the systems would be necessary to eventually corroborate or definitely discard the hypothesis of the presence of H-rich exospheres and of a sixth planet Kepler-444-g. Eventual detections of hydrogen exospheres would give some hints as to the possibility of these planets being formed with a Ganymede-like composition and evolving while retaining a substantial amount of water-ice.

As discussed in Sect. 3.1, our detailed modelling of the rotational evolution of Kepler-444-A does provide support to the hypothesis for which the HI-Lyα variations are not related to stellar activity, showing that this star is not rotating particularly fast and that its surface angular velocity is perfectly in line with what is expected for the rotation rate of similar old stars. In the following, we aim to investigate whether the observed signal can be interpreted as due to the absorption by neutral H-rich exospheres arising from the dissociation of water molecules from oceans at the surface of the planet, under the assumption of being formed with a Ganymede-like composition (Bourrier et al. 2017).

In the next section, we thus study the impact of atmospheric evaporation on the evolution of the Kepler-444-e planet and compare ourestimate for the mass loss rate and total mass loss to the ones reported in Bourrier et al. (2017). We focus our study specifically on this planet because we can benefit from the availability of the planetary mass measurement and mass loss rate estimations as constraints to our computations.

4.1 Kepler-444-e evaporation

We computed a model for the orbital evolution of Kepler-444-e, this time accounting for atmospheric evaporation. We used the prescription for planetary atmospheres mass loss in the energy-limited regime configuration (Lecavelier Des Etangs 2007; Erkaev et al. 2007)as in the following: $${\stackrel{.}{M}}_{\text{tot}}=\eta {\left(\frac{R_{\text{XUV}}}{R_{\text{p}}}\right)}^2\frac{3}{4}\frac{F_{\text{XUV}}(\text{sma})}{G{\varrho }_{\text{pl}}{K}_{\text{tide}}},$$(11)

where the factor R_XUV∕R_p accounts for the planet cross-section increase due to EUV radiation, K_tide accounts for the enhancement in the mass loss rate of close-in planets, due to the proximity of the Roche-lobe surface to the planetary surface (Erkaev et al. 2007), F_XUV is the high-energy flux emitted by the host star, ϱ_pl is the average density of the planet, and G is the gravitational constant. The factor η represents the heating efficiency, namely the percentage of stellar energy used to heat up the planetary atmospheric layers. In Bourrier et al. (2017), the estimate for the current atmospheric mass loss rate is ${\stackrel{.}{M}}_{\text{e}}^{\text{tot}}=\eta 1.7\times {10}^9$ g s⁻¹, which is calculated by taking into account their derived values for the current X-ray and EUV stellar fluxes (F_X(5−100 Å) ~ 0.1, $F_{\text{EUV}}(100-912 \AA )={1.3}_{-1.1}^{+2.6}$, in units of (erg cm⁻² s⁻¹), given at 1 AU) and setting the factors R_XUV∕R_p and K_tide to unity. They assumed that these fluxes remain constant along the whole evolution. They made the hypothesis that the planet formed beyond the ice line, with a Ganymede-like composition, for which the mean density is 1.94 g cm⁻³ and half of the planetary mass is made of water ice. Using a conservative heating efficiency (η = 0.1), they estimated Kepler-444-e to have lost approximately 27% of its initial water content after ~11 Gyr.

In the following, we present our results concerning the estimation of the Kepler-444-e atmospheric mass loss. This time, they were computed by considering the evolution of the stellar rotation and the related emitted XUV flux.

4.2 Impact of the rotational history

The hypothesis of a constant XUV flux along the evolution of the system, as described in the previous section, does not provide a realistic scenario characterising the high-energy emission of the host star since the XUV luminosity may vary over a range of several orders of magnitude along its evolution, as shown in Fig. 2. In this context, it appears fundamental to consider the rotational history of the star and the related evolution of XUV luminosity in order to have a more realistic estimation of the planetary mass loss (see e.g. Johnstone et al. 2015).

As a starting point, we used the estimation of the initial planetary mass given in Bourrier et al. (2017), namely $M_{{\text{pl}}_{\text{in}}}=0.04344 M_{\oplus }$², computed inthe hypothesis of a constant XUV flux. Considering this value as the initial mass for Kepler-444-e, we tested whether by includingan evolving XUV flux (linked to the rotational history of the host star) the planet would have survived to complete evaporation³ and if not, what would have been the timescale needed to remove the totality of its mass. We thus took into account the evolution of the stellar XUV flux as predicted by our modelling of the evolution of the system using Eqs. (7) and (8). In all these computations, the heating efficiency (η) was fixed to 0.1. Since the rotational history of the host star is unknown, we decided to first consider the case of a slow rotator (Ω_in = 3.2 Ω_⊙) in order to test the impact of the related high-energy flux emission. As a result, we determined that in this case the planet would be completely evaporated after ~50 Myr. This result shows how the inclusion of an evolving XUV leads to much stronger lifetime mass losses, even in the case of a slow rotator⁴. According to such a short timescale, it looks unlikely that the planet could have retained any atmosphere after ~ 11 Gyr of evolution for an initial mass of 0.04344 M_⊕.

We may suppose that Kepler-444-e formed with an initial mass larger than 0.04344 M_⊕. This mass indeed was estimated integrating backwards in time the mass loss rates derived in the assumption of a constant XUV flux. Such an assumption entails an underestimation of the effective XUV flux and the related mass loss rates, especially at the early stages of the evolution.

This time, we integrated backwards in time the mass loss rates obtained in the case of a slow rotator and looked for the range of initial masses that reproduce the mass of the planet at the age of the system (~ 11 Gyr). We found that this range is $\left(0.255~0.34\right) M_{\oplus }$ (see Fig. 7). At this point, we tested whether this range of initial masses was compatible with the assumption that Kepler-444-e formed with an initial Ganymede-like composition, namely as a differentiated planet with a 50% water-ice content. In a recent study, Michel et al. (2020) found that the water mass fraction of solid planets that formed beyond the ice line, orbiting a thick-disc population’s host star, represents ~ 70% of the total planetary bulk composition. We may consider this percentage as representative of the maximum fraction of the total mass loss allowed for Kepler-444-e, with respect to the initial mass.

We thus compared the range of initial masses found in the case of a slow rotator with the present mass of Kepler-444-e at the extremes of the errorbar (namely at 0.015 and 0.0921 M_⊕), and we found that it would correspond to mass losses of 94% and 73% over the lifetime of the planet, with respect to the initial mass. The value that reproduces a final mass of 0.0336 M_⊕ within this range of initial masses is $M_{{\text{pl}}_{\text{in}}}=0.28 M_{\oplus }$, which leads to a percentage of total mass loss of 88%.

This result shows that in the case of a slow rotator, we predict initial masses that are not compatible with the assumption of an initial Ganymede-like composition. Indeed, we found total mass losses that are globally above the maximum percentage ($>70\% M_{{\text{in}}_{\text{pl}}}$) allowed for such a planet. We would expect to obtain even stronger mass losses and consequently to predict larger initial masses in the case of moderate (Ω_in = 5 Ω_⊙) and fast (Ω_in = 18 Ω_⊙) rotators, given that a star remains in the saturation regime of emission longer for higher values of Ω_in. According to our results, we disfavour the hypothesis according to which Kepler-444-e has retained a fraction of its initial water-ice content after ~ 11 Gyr, due to the strong mass losses experienced by the planet, even in the case of a slow-rotating host star.

In this context, we are interested in exploring how much the lifetime mass loss of Kepler-444-e may change when considering the host star has an initial rotation rate smaller than 3.2 Ω_⊙. We may assume that the host star had an initial surface rotation rate equal to the one of the present Sun (Ω_in = Ω_⊙). Such a star is a super-slow rotator. Statistically, super-slow rotators appear to be much rarer than slow, medium, and fast rotators looking at the observed distribution of surface rotations for stars in star-forming regions and young open clusters. Based on the study by Gallet & Bouvier (2015), we indeed note that no young stars (i.e. with an age less than 4 Myr) in the mass range of 0.7−0.9 M_⊙ are observed with a surface angular velocity lower than the initial value of 1 Ω_⊙ adopted for the super-slow rotating case, and that the mean value for the 25th rotational percentile is found to be about 3 Ω_⊙ for these stars. Nevertheless, it might still be interesting to estimate the mass loss rates of Kepler-444-e, assuming this peculiar stellar rotational history. In Fig. 6, we show the evolution of the X-ray and EUV luminosity for a super-slow rotator. We notice that for such a low value for the initial surface rotation rate, the emission of high energy occurs in non-saturated conditions for the whole evolution of the system. In this case, we predict a range of initial masses of (0.147~0.22) M_⊕ to which total mass losses of 89% and 58% of $M_{{\text{pl}}_{\text{in}}}$ correspond when compared with the final masses of 0.015 and 0.0921 M_⊕, respectively.We see that in this case, we predict an initial mass of the planet (0.22 M_⊕) able to reproduce the present mass of Kepler-444-e at the upper limit of the errorbar (0.0921 M_⊕) that is compatible with an initial Ganymede-like composition, namely we predict a percentage of the total mass loss below 70%. There is actually a restricted range of initial masses that provides the same kind of result for a corresponding range of final masses: For $M_{{\text{in}}_{\text{pl}}}\in \left[0.18,0.22\right] M_{\oplus }$, we reproduce $M_{\text{fin}}\in \left[0.053,0.0921\right] M_{\oplus }$ and total mass losses below 70%.

In Fig. 7, we represent our predictions and evolution of the planetary mass for different rotational histories of the host star. In particular, we show the results obtained for the super-slow, slow, and medium rotators (yellow, cyan, and green shaded areas, respectively). For the super-slow rotator, we also indicate the range of initial masses that would be compatible with the initial composition assumed for Kepler-444-e (black hatched area).

Fig. 6 Evolution of the EUV and X-ray luminosity (solid and dashed line, respectively) computed using Eqs. (7) and (8) for the super-slow rotating case (Ωin = Ω⊙). The verticalline indicates the age of the system.

Fig. 7 Evolution of the planetary mass estimated by considering different rotational histories for η = 0.1 and RXUV∕Rpl = 1. The yellow shaded area shows the evolution for Ωin = 1 Ω⊙, the cyan shaded area for Ωin = 3.2 Ω⊙, and the green shaded area for Ωin = 5 Ω⊙. The black hatched area shows the evolution for Ωin = Ω⊙, which is the range of initial masses compatible with an initial Ganymede-like composition.

4.3 Impact of XUV-luminosity prescriptions

As mentioned in previous sections, the rotational history of a given star is unknown. Nevertheless, for stars hosting planetary systems, it may be possible to infer some indications about preferential rotational histories by studying their direct impact on the evolution of the planets (Kubyshkina et al. 2019). In the previous section, we found that in the case of a slow rotator (Ω_in = 3.2 Ω_⊙), we predict initial masses for Kepler-444-e that are not compatible with an initial Ganymede-like composition. We thus considered the particular case of a super-slow rotator, finding a restricted range of initial masses ($0.18\lesssim M_{{\text{in}}_{\text{pl}}}/M_{\oplus }\lesssim 0.22$) with which our results are compatible, with total mass losses below the 70% threshold of the initial mass. According to this result, we might establish that the super-slow rotator case would constitute a preferential pathway for the rotational history of the host star.

However, there is another important physical ingredient intervening in the evaluation of the atmospheric mass loss that may change the estimation of the initial planetary mass significantly, namely the XUV-luminosity prescription. It is possible that the estimation of the planetary mass loss might vary considerably, depending on the XUV-luminosity prescription used, even if the same rotational history is considered for the host star. So far, we have used Eq. (7) for the computation of the X-ray luminosity and Eq. (8) for the derivation of the corresponding EUV luminosity (Wright et al. 2011; Sanz-Forcada et al. 2011). It is worth investigating how the results obtained for the evaporation of Kepler-444-e in the case of a super-slow rotator may change by using an alternative XUV-luminosity prescription.

For this purpose, following the work of Johnstone et al. (2021), we adopted new prescriptions for the computation of the X-ray and EUV luminosities (see Eqs. (9) and (10) in Sect. 2.2) and using the same physics set-up defined in Sect. 4.2, we recomputed the evolution of the planetary mass. In the top panel ofFig. 8, a comparison between the evolution of the different XUV fluxes is presented, the solid red and yellow lines being the XUV flux computed by following Johnstone et al. (2021) and Wright et al. (2011); Sanz-Forcada et al. (2011), respectively, while the horizontal black-dashed line represents the estimation of the XUV flux at the age of the system as reported in Bourrier et al. (2017). For the sake of completion, we also show the evolution of the XUV flux computed as in Jackson et al. (2012), reminding us that in this case it does not directly depend on the rotational history of the host star. One of the main features emerging from this plot is that the slope of the red line is flatter compared to the one of the yellow line, especially in the unsaturated regime, leading to an estimation of the XUV flux at the age of the system that is almost one order of magnitude larger than the one indicated by the yellow line and by Bourrier et al. (2017). As a consequence, we expect to get a much more efficient atmospheric mass loss when using the new XUV-luminosity prescription in the evaporation routine. Indeed, in this case, we predict that the initial mass of the planet should have ranged between roughly 0.28 and 0.35 M_⊕ in order to be able to reproduce its mass at the age of the system. This range of initial masses implies total mass losses that are above the maximum percentage for being compatible with an initial Ganymede-like composition.

In the light of this result, it seems that there is no chance for the planet to retain a fraction of its initial water-ice content with the new XUV-luminosity prescription (even in the case of a super-slow rotator) if we assume it formed with a Ganymede-like composition, unless some physical mechanisms intervene reducing or inhibiting the mass loss.

Fig. 8 Top panel: comparison among the evolutions of XUV fluxes computed by following the work of Wright et al. (2011; solid yellowline), Jackson et al. (2012; solid green line), and Johnstone et al. (2021; solid red line). The horizontal black-dashed line represents the estimation of the XUV flux of Bourrier et al. (2017). Bottom panel: evolution of the planetary mass computed by considering different XUV-luminosity prescriptions. The colour code corresponds to the one described in the top panel.

4.4 Impact of the effective XUV-absorption radius

It is worth noting that the estimations of the atmospheric mass loss presented so far rely on the assumption that the factors (β = R_XUV∕R_pl) and (K_tide) in Eq. (4.1) are set to unity, which is consistent with Bourrier et al. (2017). However, accounting for their variations as a function of time might have a significant impact on the final results. In particular, accounting for the presence of an absorption radius located at larger distances with respect to the surface of the planet increases the effective area heated by the stellar high-energy radiation, consequently leading to a more efficient evaporation.

Here, we consider the contribution of a factor of β > 1 in the estimation of the planetary lifetime mass loss. For the computation of this factor, we refer readers to the work of Salz et al. (2016).

In the light of the results obtained in the previous sections, we aim to study whether the particular set of conditions that we found to favour the hypothesis of an initial Ganymede-like composition for β = 1 are still valid for β > 1. These conditions require that the host star evolved as a super-slow rotator, emitting XUV radiation as predicted by Eqs. (7) and (8), and that the planet had an initial mass in the range of (0.18−0.22) M_⊕.

We thus recomputed the evaporation of the planetary mass, this time including the contribution of the effective XUV-absorption radius. The results related to these computations are shown in Fig. 9 (solid red lines). We notice that the impact of the R_XUV radius is significant, leading to much stronger evaporation rates that would remove the totality of the planetary mass within the first ~50 Myr of evolution. Such strong mass losses are incompatible with the hypothesis that Kepler-444-e could have retained any atmosphere after 11 Gyr of evolution if it formed with a Ganymede-like composition. The simulations obtained for β > 1 thus show that, even in the case of a super-slow rotating star, the evaporation rates are strong enough to efficiently remove a mass content as large as the mass of Kepler-444-e in relatively short timescales.

Fig. 9 Evolution of the planetary mass computed with β > 1, Ωin = Ω⊙ and XUV formalism as in Eqs. (7) and (8). The solid red and black lines show the evolution for η = 0.1 and eta = 0.01, respectively. The black shaded area indicates the range of initial masses compatible with a lifetime mass loss below $70\% M_{{\text{in}}_{\text{pl}}}$.

4.5 Impact of the heating efficiency

As a last point, we would like to consider the role of the heating efficiency. In our study, we treated the evaporation using a value for the heating efficiency of 10% to compare our results with the ones of Bourrier et al. (2017). Lopez (2017) indicates that while this value is usually used for planets with solar composition atmospheres, in the case of higher metallicity atmospheres it should be significantly lower. This is due to the higher mean-molecular weight that would reduce the scaleheight in the evaporative wind and increase the number of scaleheights between the XUV photosphere and the sonic point; additionally, this is due to the fact that in this regime the metal atomic lines are dominant coolants (Lopez 2017). Applying the scaling relation found by Ercolano & Clarke (2010) between heating efficiency and atmospheric metallicity (η ∝ Z^−0.77) to pure water atmospheres, Lopez (2017) finds a value of η closer to 1%. The same scaling relation is also used in Mordasini (2020). In this work, the author points out that the results based on this scaling should be considered with some caution, recalling that the relation of Ercolano & Clarke (2010) is derived from studies of X-ray heating of protoplanetary disc atmospheres rather than planetary atmospheres. Moreover, he argues that while the evaporation rate seems to likely decrease with increasing metal content (given the dominant role of atomic lines as coolants), the exact dependencies have not been explored yet (Mordasini 2020).

Given the uncertainties related to the value of the heating efficiency, we recomputed the evolution of the planetary mass for a slow and super-slow rotator with β > 1, but this time forη = 0.01. As a result, we found that in the case of the slow rotator, we predict initial masses that are not compatible with a Ganymede-like composition, even for such a low value of η. Instead, for the super-slow rotator case, we found that the lower value of the heating efficiency rather compensates for the impact of R_XUV, leading us to recover a range of initial masses (0.175~0.20 M_⊕) that reproduce the mass of Kepler-444-e at the age of the system (see Fig. 9, black hatched area) for a percentage of the total mass loss lower than 70%.

As stated before, some caution is required when using a heating efficiency as low as 0.01 in the case of metal rich atmospheres. Moreover, we recall that these simple evaporation computations tend to underestimate the mass loss rate, especially at the early stages of the evolution when the radius of the planet would be much more inflated (Lopez & Fortney 2014; Chen & Rogers 2016; Zeng et al. 2016), allowing for a larger effective absorption area of the stellar high-energy radiation together with a decrease in the gravitational well for escaping particles. More importantly, we note that if we compute the value of the restricted Jeans escape parameter⁵ of Kepler-444-e at the beginning of the evolution, assuming a maximal value for the initial mass M_pl = 0.22 M_⊕, a minimum value for the radius R_pl = 0.559 R_⊕, and a Bond albedo A = 0.5, we obtain Λ ≃4. For such a low value of the restricted Jeans escape parameter (Λ ≲ 20, Fossati et al. 2017; Kubyshkina et al. 2018), the atmospheric escape occurs in the so called boil-off regime, in which the escape is driven by a combination of intrinsic thermal energy and low gravity. During this phase, the energy-limited formula has been shown to significantly underestimate the mass loss rates (Lammer et al. 2016; Owen & Wu 2016; Stökl et al. 2016; Fossati et al. 2017). We thus see that, even in the most favourable case of a low heating efficiency of 0.01, the Ganymede hypothesis is disfavoured by evaporation computations for Kepler-444-e once the evolution of the stellar rotation and the related emitted XUV flux is takeninto account.

5 Conclusion

In this paper we aimed to study the past history of the planets orbiting Kepler-444-A, using all the constraints available for this system.

We computed rotating stellar models of Kepler-444-A using the stellar parameters provided in Paper I and selecting three different values for the initial surface rotation rate, which are representative of initially slow, moderate, and fast rotators. We correctly reproduced the observed surface rotation rate of Kepler-444-A, finding that this star does not exhibit a particularly fast rotation compared to similar old stars. This result suggests that there is no reason for this star to be exceptionally active, favouring the hypothesis that the HI-Lyα variations observed for this system are not related to stellar activity.

Coupling the rotating stellar models to our orbital evolution code, we studied the impact of dynamical and equilibrium tides on the evolution of Kepler-444-d and Kepler-444-e. Given the low mass ratio between planets and the host star and their relatively large orbital distances, we found that tides do not play a significant role in shaping the architecture of this system.

Afterwards, we studied the impact of atmospheric evaporation on Kepler-444-e. We investigated whether our estimated lifetime mass losses are compatible with the formation hypothesis of the planet with a Ganymede-like composition retaining a fraction of the initial water-ice content (representing up to the 70% of the initial mass) at its current age, allowing the presence of a H-rich exosphere trailing the planet.

We computed XUV-flux evolutionary tracks for a slow rotator and estimated the mass loss rates using the energy-limited formalism. Starting with the initial planetary mass derived in Bourrier et al. (2017), we found that our mass loss rates would remove the entire mass within the first 50 Myr of evolution. Integrating the mass loss rates obtained for a slow rotator backwards in time, we then predicted the range of initial masses that allowed us to reproduce the current mass of Kepler-444-e. However, we find total mass losses above the maximum percentage of the initial water-ice content allowed for a Ganymede-like planet. The same result applies to the case of medium and fast rotators to which larger total mass losses correspond.

We tested whether the same result holds for a star having a very low initial surface rotation rate (super-slow rotator). We found that for such a peculiar rotational history, it is possible to predict initial masses compatible with an initial Ganymede-like composition, but only for a particular set of conditions, requiring that the X-ray and EUV luminosities are computed following the works of Wright et al. (2011) and Sanz-Forcada et al. (2011) and that the effective XUV-absorption radius is equal to the radius of the planet (β = 1). This solution does not hold any more when using an alternative XUV-luminosity prescription, such as the one derived by following the work of Johnstone et al. (2021); similarly, when including the impact of the effective XUV-absorption radius (β > 1), we estimatetotal mass losses that are above the maximum percentage allowed for such a planet to be compatible with an initial Ganymede-like composition. Only when adopting a value of η as low as 1%, in the case of β > 1, we do retrieve a solution for the initial planetary mass compatible with our assumption on the initial composition, but only for a super-slow rotator and for XUV luminosities computed following Wright et al. (2011); Sanz-Forcada et al. (2011).

We note that the need to adopt an exceptionally super-slow rotating star and/or a very low value for the heating efficiency to be in agreement with the initial composition assumed for Kepler-444-e point towards disfavouring the hypothesis that the Ly-α variations observed for this system are related with the presence of H-rich escaping atmospheres today, arising from the oceans at the planet surface. We recall that our rotating models also do not support a stellar origin for the observed HI-Lyα variations, thus the origin of this detected variability remains an open question.

Acknowledgements

We thank the referees for their useful comments that helped to improve the quality of the manuscript. This projecthas been supported by the Swiss National Science Foundation grant 200020-172505. G.B. acknowledges funding from the SNF AMBIZIONE grant No 185805 (Seismic inversions and modelling of transport processes in stars). P.E. and G.M. have received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No 833925, project STAREX). V.B. acknowledges support by the Swiss National Science Foundation (SNSF) in the frame of the National Centre for Competence in Research “PlanetS”. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project Four Aces grant agreement No 724427; project Spice Dune grant agreement No 947634).

References

All Tables

All Figures

	Fig. 1 Evolutionary tracks in the HR diagram obtained with the GENEC (in blue) and CLES (in salmon) evolution codes together with the location of the optimal asteroseismic solution (red dot). The tracks show the evolution from the PMS to part of the MS phase. The GENEC track in particular corresponds to the moderate rotation case.
In the text

Fig. 2 Top panel: surface angular velocity evolution for Kepler-444 in the case of a fast rotator (Ωini = 18 × Ω⊙, solid blue line), a medium rotator (Ωini = 5 × Ω⊙, solid magenta line), and a slow rotator (Ωini = 3.2 × Ω⊙, solid black line). The red dot represents the surface angular velocity of Kepler-444-A at the age of the system (11 Gyr). The red-dashed line shows the critical velocity limit. The vertical black-dotted line indicates the age of the system, while the horizontal red-dotted line shows the value of the surface rotation at the age of the system. Bottom panel: evolution of the EUV (solid lines) and X-ray luminosities (dashed lines) computed according to Eqs. (7) and (8) relative to the different rotational histories considered. The vertical black-dotted line indicates the age of the system.

In the text

Fig. 3 Top panel: orbital evolution of Kepler-444-d and e (solid red and green line, respectively), starting with the current value of the orbital distance ($a_{{\text{ini}}_{\text{d}}}=0.06$, $a_{{\text{ini}}_{\text{e}}}=0.0696$), at a fixed mass (Md = 0.1016 M⊕, Me = 0.0921 M⊕) and host star initial angular rotation Ωin = 18 × Ω⊙. The vertical red-dotted line corresponds to the age of the system (age = 11 Gyr). The blue-dashed lines indicate the evolution of the corotation radii, the blue-dotted ones indicate the evolution of the minimum initial orbital distances below which dynamical tides are no longer active. The magenta area represents the extension of the stellar convective envelope, while the cyan area represents the stellar radiative core. Bottom panel: orbital angular velocity of the planets and the double of the host star angular velocity (blue-dashed line) as a function of time. Dynamical tides are active as long as the orbital frequency of the planets is lower than double the host star surface angular rotation.

In the text

Fig. 4 Impact of a change in the planetary mass on the orbital evolution for a planet located at the same distance as Kepler-444-d (initial value fixed to $a_{{\text{ini}}_{\text{d}}}=0.06$ AU). The four red lines span a range of masses from 10 to 1000 M⊕, with larger masses corresponding to more significant increases in the distances. The dashed lines, dotted lines, and shaded regions are the same as in Fig. 3 (top panel).

In the text

	Fig. 5 Impact of a change in the initial orbital distance for a planet with the maximal value of the mass allowed for Kepler-444-d (mass of the planet fixed to 0.1016 M⊕). The dashed lines, dotted lines, and shaded regions are the same as in Fig. 3 (top panel).
In the text

	Fig. 6 Evolution of the EUV and X-ray luminosity (solid and dashed line, respectively) computed using Eqs. (7) and (8) for the super-slow rotating case (Ωin = Ω⊙). The verticalline indicates the age of the system.
In the text

Fig. 7 Evolution of the planetary mass estimated by considering different rotational histories for η = 0.1 and RXUV∕Rpl = 1. The yellow shaded area shows the evolution for Ωin = 1 Ω⊙, the cyan shaded area for Ωin = 3.2 Ω⊙, and the green shaded area for Ωin = 5 Ω⊙. The black hatched area shows the evolution for Ωin = Ω⊙, which is the range of initial masses compatible with an initial Ganymede-like composition.

In the text

Fig. 8 Top panel: comparison among the evolutions of XUV fluxes computed by following the work of Wright et al. (2011; solid yellowline), Jackson et al. (2012; solid green line), and Johnstone et al. (2021; solid red line). The horizontal black-dashed line represents the estimation of the XUV flux of Bourrier et al. (2017). Bottom panel: evolution of the planetary mass computed by considering different XUV-luminosity prescriptions. The colour code corresponds to the one described in the top panel.

In the text

	Fig. 9 Evolution of the planetary mass computed with β > 1, Ωin = Ω⊙ and XUV formalism as in Eqs. (7) and (8). The solid red and black lines show the evolution for η = 0.1 and eta = 0.01, respectively. The black shaded area indicates the range of initial masses compatible with a lifetime mass loss below $70\% M_{{\text{in}}_{\text{pl}}}$.
In the text