© ESO 2021

1 Introduction

Numerous exoplanets, ranging from small rocky objects to planets more massive than Jupiter, have been discovered orbiting very close to their host¹. The vicinity between the planet and the star both affects and determines the various interactions between them, leading to significant consequences for the evolution of the star–planet system: a planet can orbit so close to its host star (planet grazing) that it can stimulate gas ejections from the star (Stephan et al. 2020). Such a scenario has been reported for the red giant star V Hydrae (Sahai et al. 2016; Salas et al. 2019). The star may be interacting with the magnetic field of the planet, giving way to hot spots on its surface (see e.g., Strugarek et al. 2019) and also inducing some magnetic activity with regard to the planet. Stellar rotation is affected by transfers of angular momentum between the planetary orbit and the star through equilibrium and dynamical tides (Zahn 1977; Ogilvie & Lin 2007; Ogilvie 2013; Mathis 2015). The irradiation received by the planet from the star impacts its equilibrium temperature and thus the hydrostatic structure of its envelope (Yelle 2004; Baraffe et al. 2004; Bear & Soker 2011; Debrecht et al. 2020). The energy deposited by high-energy photons from the star may lead to inflation of the planet atmosphere (as with e.g., WASP-121b Delrez et al. 2016; Bourrier et al. 2020) and even to the partial or total loss of the planetary atmosphere (Ehrenreich et al. 2015). Interestingly, in the future, using the JWST (James Webb Space Telescope), observations of the planetary evaporation process through transmission spectroscopy during transits might allow us to infer the composition of the interior of small rocky planets (Bodman et al. 2018). Some observed features, such as the existence of the radius valley and of the Neptunian desert in the plane of planet radius versus orbital period (see e.g., Van Eylen et al. 2018) have been interpreted as having been shaped, to a great extent, by the planetary evaporation process (Lopez & Fortney 2013a,b; Owen & Wu 2013; Jin et al. 2014; Jin & Mordasini 2018; Owen & Lai 2018; Owen 2019; Owen & Adams 2019). The distance between the star and planet is modified by the tides and also by changes of the planetary and stellar masses. This may cause a planetary engulfment at some time during the evolution of the system (Siess & Livio 1999a,b; Villaver & Livio 2007; Nordhaus et al. 2010; Kunitomo et al. 2011; Mustill & Villaver 2012; Nordhaus & Spiegel 2013; Privitera et al. 2016c; Stephan et al. 2020). A planetary engulfment can be accompanied by transient phenomena such as strong stellar winds as well as emissions in the visible, UV and X-ray radiation (Metzger et al. 2012; Stephan et al. 2020) that may last for durations of 100–1000 yr. Some authors (Morris 1981; Soker 1998; Kim et al. 2017) have also suggested that planetary engulfment can cause non-spherical, dipole-shaped planetary nebulae. Planetary engulfment can also cause the star to spin up (Privitera et al. 2016c; Qureshi et al. 2018; Stephan et al. 2020).

A number of previous works have studied the impact of tides on planet orbits (see e.g., Livio & Soker 1984; Soker et al. 1984; Sackmann et al. 1993; Rasio et al. 1996; Siess & Livio 1999a,b; Villaver & Livio 2007; Sato et al. 2008; Villaver & Livio 2009; Carlberg et al. 2009; Nordhaus et al. 2010; Kunitomo et al. 2011; Bear & Soker 2011; Mustill & Villaver 2012; Nordhaus & Spiegel 2013; Villaver et al. 2014; Privitera et al. 2016b,c,a; Meynet et al. 2017; Bolmont et al. 2017; Gallet et al. 2018) without considering – for the most close-in planets – the effects of planetary evaporation. Other studies haveconsidered planetary evaporation but do not account for the impact of tides induced by the planet on the star (the only tidal effect accounted for is the Roche lobe overflow from the planet; see e.g., Penz et al. 2008; Kurokawa & Nakamoto 2014; Ionov et al. 2018). We note that in most of the cases, there is a reasonable justification for neglecting tidal interactions.Indeed, as is well known (see e.g., a review in Owen 2019, and references therein), evaporation mostly affects low-mass planets whose evaporation timescale is generally much shorter than the timescale of their orbit being affected by tides. For the most massive planets, however, this neglect is less justified. In a few studies, both the processes of evaporation and tides have been accounted for (see e.g., Jackson et al. 2016; Collier Cameron & Jardine 2018; Owen & Lai 2018). However, in these works, it is only equilibrium tides or Darwin’s classical theory of tides that have been considered, even thoughthese are not valid for fast rotating stars where dynamical tides may dominate. In addition, most studies haveonly considered the average X-ray and Ultra-Violet radiation (XUV) luminosity coming from the star, not accounting for the possibility that this XUV luminosity can vary with the initial rotation of the star, which itself followsa natural evolution that can also be modified due to tidal interaction in the case of more massive planets.

In the present work, we take our first steps toward a more coherent and consistent approach, in simultaneously considering the evolution of the planet orbit and mass, together with the evolution of the star. At the moment, to our knowledge, there are very few works, if any, that account, in a consistent way, for the effects of tides and evaporation, including the evolution of the star, whose rotational evolution is impacted by these processes. In this first step, however, we consider a simple model for the planet, which is likely to be too schematic but sufficient in our view to explore possible feedbacks between tides, evaporation,and stellar rotation. Moreover, comparisons with observations and with more sophisticated models for the planet indicate that despite the simple planetary models we use here, the results we have obtained do not appear unrealistic.

Since this paper is the sixth in a series, we briefly compare the physics studied here with that of our past works. In the first set of papers (Privitera et al. 2016a,b,c), we considered the case of planets that are sufficiently distant from their host star to have a tidal interaction only during the red giant phase. Here, only the case of equilibrium tides was considered. We showed that planetary engulfment may cause an increase in the surface velocity of red giants well above what can be achieved even by considering the most extreme cases of initial rotation and angular momentum transport efficiency. We confirmed that planetary engulfment could be responsible for increasing significantly the lithium abundance in the envelope of red giant stars (Meynet et al. 2017) and made predictions of the surface magnetic field of red giant stars after spin-up by planetary engulfment (Meynet et al. 2017). Then, in Rao et al. (2018), we studied the cases of close-in planets to the host star during its pre-main-sequence (PMS) phase, accounting for dynamical tides along with equilibrium tides. We concluded that the region where the planet can be engulfed depends sensitively on the rotation of the star.

In the present work, we add a new aspect in our approach by accounting for the possible evaporation of the planet. This would mainly affect, of course, planets at a close distance to their host star and mainly concern the PMS phase and the early main-sequence (MS) phase, when the star is still rotating rapidly enough to emit high-luminosity XUV radiation. As indicated above, we study the interactions between the following processes: changes of the planetary orbit due to equilibrium and dynamical tides induced by the planet on the star; changes of the planetary mass due to the XUV irradiation of the planet by the star; and changes of the stellar rotation rate due to the stellar interior evolution, magnetic wind braking, tides, and planetary engulfment. These three aspects are tightly linked. To illustrate this, let us consider a planet sufficiently close to its star so that its distance from the host star shrinks due to tidal interactions (i.e., implying that the orbital distance is inferior to the corotation radius). This is expected to lead to at least two obvious effects: First, it will increase the irradiation flux received by the planet closer to the star and, secondly, since the star will receive some angular momentum, it will rotate faster and, therefore, its activity will vary from what it would be had it been evolving in isolation. This would also change the XUV irradiation and, therefore, the evaporation rate of the planetary atmosphere. The mass loss of the planet impacts, in turn, the evolution of the orbit, thus impacting as well the amplitudes of the tidal interactions between the star and the planet.

Within the limits of the hypotheses made in the present work (see Sect. 2), in this paper we aim to consider the following questions of: (1) how the stellar rotation and stellar activity react to the changes of the orbit and how these changes affect the irradiation of the planet and the photo-evaporation of its atmosphere; (2) in the plane of initial orbital distance versus initial mass of the planet, how to locate those regions where tides dominate, where planetary evaporation dominates, and where both aspects play a significant role in impacting the planetary orbit and mass; (3) how the initial rotation of the star impacts the fate of the planets.

This paper is organized as follows: Sect. 2 introduces the physics of the models. We discuss the interplay between tides and the evaporation process in Sect. 3. We summarize the different types of evolution of the star–planet systems considered in this study, depending on the initial conditions, in Sect. 4. A comparison of the predictions of our present simulations with observations is presented in Sect. 5. We discuss the advantages and limitations of the present approach in Sect. 6. The main results are synthesized in Sect. 7.

2 Physics of our computations

2.1 Stellar models

We used a 1 M_⊙ star model computed by the GENEC code (Eggenberger et al. 2008). The computation of the model started from the PMS phase when its radius was about 1.8 R_⊙. The time toreach the zero age main sequence (ZAMS) from this point is 40 My. A full description of this 1 M_⊙ star model in a graphical form is given in Fig. B.1 of Rao et al. (2018).

We consider that the star exhibits, at all times throughout the PMS and MS phases, a solid body rotation. This seems to be a reasonable hypothesis in view of the rotation profile deduced for the Sun from helioseismology (García et al. 2007; Eggenberger et al. 2019)and seismic analyses of solar-like MS stars using Kepler data (Nielsen et al. 2015; Benomar et al. 2015, 2018) and γDor stars (Saio et al. 2021). Hence, the angular momentum of the star, L_st, in terms of its moment of inertia, I_st, and stellar rotation, ω_st, is given by, $$L_{\text{st}}=I_{\text{st}}{\omega }_{\text{st}}.$$

We decouple the evolution of the structure of the star and, therefore, of I_st (by computing it with our stellar model independently for a given choice of the initial rotation) from the evolution of its rotation, ω_st, which depends on the changes in the stellar interior structure, tidal interactions, and magnetic interactions. We assume, therefore, that the feedback of rotation on the internal structure of the star is negligible. This is very well justified as long as the rotation is below about 70% of the critical velocity at any point in the star. It happens that this condition is never violated in our case, or if so, only at the surface and for a very short time. Thus, we can confidently say that indeed the structure evolution of the star can be taken to be the same for a very wide range of rotations.

Overall, the equations describing the evolution of the angular momentum of the star, L_st, as a function of time are as follows: $${\dot{L}}_{\text{st}}={\dot{L}}_{\text{wb}}+{\dot{L}}_{\text{ti}},$$(1)

where ${\dot{L}}_{\text{wb}}$ is the rate of change of angular momentum resulting from magnetized stellar winds according to the prescription from Matt et al. (2015, 2019). As in Eggenberger et al. (2019), we fixed the constant defining the transition between the unsaturated to the saturated regime to 10.

The quantity ${\dot{L}}_{\text{ti}}$ expresses the rate of change of the stellar angular momentum due to the tides induced by the planet in the star: $${\dot{L}}_{\text{ti}}=-\left[\frac{1}{2}{m}_{\text{pl}}{\left(\frac{\dot{a}}{a}\right)}_{\text{ti}}\right]\cdot \sqrt{G\left(M_{\text{st}}+m_{\text{pl}}\right)a},$$

where m_pl is the mass of the planet, a is the orbital distance, G is the gravitational constant, and M_st is the mass of the star. The expression for ${(\dot{a}/a)}_{\text{ti}}$, namely, the secular change in orbit due to dynamical and equilibrium tides, is discussed in Sect. 2.3.

Our stellar model will be used to infer the X-ray (λ ~ 5–100 Å) and Ultra-Violet (λ ~ 100–920 Å) luminosity of our star. Each consists of two contributions: the first is the thermal contribution due to the black-body radiation of the star which can be calculated knowing the total luminosity by multiplying it with a factor which is a Bose-Einstein integral that comes from integrating the Planck black-body radiation law in the XUV regime. The second is a non-thermal contribution composed of an X-ray and ultra-violet (UV) component that arises due to the radiation emitted by the charged stellar-wind particles, which get accelerated by the star’s magnetic field.

The non-thermal X-ray luminosity, L_X, is computed as in Tu et al. (2015), L_X = R_XL_tot, with $$R_{\text{X}}={10}^{-3.13}{\left(\frac{R_{\text{o}}}{0.196}\right)}^{-3.25}$$

if R_o < 0.196. For R_o ≥ 0.196, R_X = 10^−3.13. Here, R_o is the Rossby number². The expressions above reflect the fact that when the rotation increases, i.e., R_o decreases), then R_X increases and thus, so does the X-ray luminosity. This increase however, saturates below some critical Rossby number, here equal to 0.196.

The non-thermal UV contribution can be obtained via an empirical relation between the UV and X-ray contributions (Sanz-Forcada et al. 2011).

In the left panel of Fig. 1, the evolution of the Rossby number is shown as a function of time for 1 M_⊙ stellar models with different initial rotations. We see that for most of the present models, the Rossby number is below 0.196. This means that the X-ray luminosity will follow the total luminosity scaled with a constant factor equal to 0.00074. The luminosities in X-rays, UV, and their sum is shown in the right panel of Fig. 1. Saturation is always reached for the fast rotating model, while for the slow rotating one, it occurs only for ages below about 27 Myr. During the saturation period, the red and blue curves are superposed, so only the red curve is shown in Fig. 1.

2.2 Planetary mass evaporation rates

The evaporation rate is taken as in Erkaev et al. (2007): $${\dot{m}}_{\text{pl}}=\frac{ϵ\pi {(\beta r_{\text{pl}})}^2{F}_{\text{XUV}}}{\Delta \Phi },$$(2)

where ϵ is the heating efficiency. This quantity, ϵ, is smaller than unity because part of the energy can be used for ionisation. In addition, various processes, such as recombinations, collisional excitations, etc., can cause radiative cooling that will reduce the energy available for the heating. Here we assume ϵ = 0.4, as suggested by Valencia et al. (2010). These authors showed that this value is adapted to the cases of hot Jupiters and strongly irradiated rocky planets. β is a number larger than unity, accounting for the inflation of the envelope due to its heating by the irradiation of the star, and β is taken as in Salz et al. (2016), while r_pl is the planet radius. We assume a spherical planet of uniform and constant density. The initial planet radius is computed from the initial mass of the planet using the mass-radius relation found by Bashi et al. (2017). This relation is a piece-wise power law relation obtained by using advanced statistical techniques to fit the power law with direct observational exoplanet data and with numerical simulations of planet formation (see Fig. B.1). The calculation of the irradiation flux, F_XUV, is given by L_XUV∕(4πa²), where a is the orbital distance of the assumed quasi-circular orbit of the planet; ΔΦ is the energy required to lift one unit of mass up to the Roche lobe and is taken as in Erkaev et al. (2007). We applied this evaporation law for planet masses in the range from 6 up to 800 Earth masses (0.02 and 2.5 Jupiter masses) as in Salz et al. (2016). We note that under the present evaporation law, planets more massive than 2.5 M_Jup would lose only a negligible amount of mass and, thus, their evolution would correspond to a constant mass evolution.

The change of the planetary radius (under the assumption of constant bulk density of the interior structure in our simple planet model) is obtained from: $${\dot{r}}_{\text{pl}}=\frac{1}{3}\frac{{\dot{m}}_{\text{pl}}}{m_{\text{pl}}}{r}_{\text{pl}}.$$(3)

Fig. 1 Left panel: evolution of the Rossby number for an isolated 1 M⊙ stellar model during the PMS phase and the very early MS phase for different values of the initial rotation. The initial rotation periods and velocities can be seen in Fig. A.1. The Rossby number corresponding to the slowest initial rotation is the upper curve. Right panel: evolution of the X-ray, EUV luminosities and of their sum for a slow (in blue) and a fast (in red) rotating 1 M⊙ model. The slow and fast rotations correspond, respectively, to the slowest and fastest model shown in Fig. A.1.

2.3 Orbital evolution

Here, we consider the case of a planet orbiting in a quasi-circular orbit around its host star. The secular change of the orbital distance due to equilibrium and dynamical tides produced by the planet on the star is considered. Tides are accounted for only when convective zones are present in the star. We note that since we are considering the case of a 1 M_⊙ star, an outer convective envelope is present during the whole evolution. We also account for the secular change of the orbital distance dueto change in planetary mass from mass evaporation. The expression for the net secular change in the planetary orbit is here written as: $$\left(\frac{\dot{a}}{a}\right)=-\frac{{\dot{m}}_{\text{pl}}}{m_{\text{pl}}+M_{\text{st}}}+{\left(\frac{\dot{a}}{a}\right)}_{\text{eq}}+{\left(\frac{\dot{a}}{a}\right)}_{\text{dy}},$$(4)

where the first term is due to planetary mass evaporation and the last two terms, as these are the contributions attributed to the equilibrium tides and dynamical tides, respectively; together they are equal to the quantity ${(\dot{a}/a)}_{\text{ti}}$ mentioned in Sect. 2.1. All these quantities are computed as in Rao et al. (2018)³. In the present work, we neglect changes of the stellar mass due to stellar winds and, therefore, the planetary drag and the planetary mass accretion.

2.4 Numerical procedure

In comparison to our previous work (Rao et al. 2018), the present version of the code mainly differs as per the following points: Here, we account for the planetary evaporation process and we consider that the star rotates as a solid body during the PMS phase and MS phases. The magnetic braking law has also been updated (see Sect. 2). The time-step variations have been improved to allow the code to evaluate more accurately the evolution of very close-in planets. The present code solves a system of four coupled ordinary differential equations: for the stellar angular momentum, the mass and radius of the planet, and the orbital distance, represented by Eqs. (1), (2), (3), and (4), respectively. Since time does not appear explicitly in any of these equations, we use the explicit classic fourth-order Runge-Kutta method for solving this system of equations.

Fig. 2 Evolution as a function of time of the orbits of 1 MJup mass planetsaround a 1 M⊙ star, obtained using different approaches. The blue continuous lines are obtained in the present work by switching off the planetary evaporation process. The red dashed lines are obtained switching on the evaporation process. The orange dashed lines (only two cases shown) are the orbits obtained in Rao et al. (2018). The blue dashed and dotted lines show respectively, the corotation radius (acorot) and the critical distance below which dynamical tides are no longer active (amin). They are obtained for the orbit beginning its evolution at a distance equal to 0.034 au. The upper-limit of the pink zone corresponds to the radius of the star in au. Its lower limit shows the radius at the base of the stellar convective envelope. The blue zone indicates the size of the stellar radiative interior.

3 Photoevaporation and tides

3.1 Consistency check of the new version of our code

In order to see how the results obtained with the new version of our code compare to the results presented in Rao et al. (2018), we plotted a few orbits from Rao et al. (2018) for a 1 M_Jup mass planet around a 1 M_⊙ star together with those computed with the present version, shown in Fig. 2.

Comparing the dashed orange lines (models by Rao et al. 2018) and the continuous lines (current models without evaporation), we see only very small differences linked to the models used for the stellar rotation (non-solid body rotation in Rao et al. 2018, and solid body rotation in the present work) and also due to the different numerical procedure used in the present version. Concerning the latter point, we see that the orange curve for the planet beginning its evolution at 0.03 au stops at a much earlier stage than the blue one. This is because in the old version, due to a less sophisticated way of handling the variation of the time steps, the orbit plunged into the star very early on. Because of this, as was already explained in Rao et al. (2018), we did not show that part of the orbit (see Fig. 2 in Rao et al. 2018). In the new version, this point has been improved and that is why the orbit can be computed on a longer timescale. Earlier, we were unsure of the fate of the planet corresponding to this computation, but the present calculation shows that in the absence of any evaporation process, the planet would be engulfed after a few 100 Myr. Based on this discussion above, we conclude that there is consistency between the results we obtained in Rao et al. (2018) using the old version of the code and the present results obtained using the new version.

Fig. 3 Evolution for a given planetary mass and for different initial radii of the orbit of the orbital distance as a function of time (first column panel), of the planetary mass as a function of time (second column), of the planetary radius as a function of the orbital period (third column), and of the evaporation rate as a function of orbital period (fourth column). Identical colors in a set of row panels correspond to models with the same initial distance to the star (see the left panel in each row). We note that the panels of the fourth row in the second and third column have a different ordinate.

3.2 Impact of evaporation and tides on a Jupiter-mass planet

The magenta lines in Fig. 2 shows the evolution of the orbits of a Jupiter-mass planet when the evaporation processes are accounted for. At first glance, there are not many differences between the evolution of the planetary orbits with and without evaporation. However, as shown below (see the panel in the second column, second row of Fig. 3), depending on the initial distance of the planet to the star, we can see a significant change of the mass of the planet for distances below 0.04 au. As a numerical example, we find that the planet, having begun its evolution at a distance of 0.03 au, would be obliterated at an age of 20 Myr due to evaporation and not be engulfed at 100 Myr as predicted by non-evaporating model. We assume here a completely gaseous planet, which is not very realistic, but we may interpret that the photo-evaporation process will be efficient enough to completely remove the planetary atmosphere at least. Afterwards, the remaining rocky core of the planet may pursue its evolution. Evaporation may have thus a strong impact for Jupiter-mass planets whose orbits shrink as a result of the dynamical tides. The present simulations indicate that all planets of this mass beginning their evolution at a distance less than 0.034 au are evaporated during the PMS phase (this number is valid for the considered initial rotation of our 1 M_⊙ star model).

The evaporation process can be accelerated by the action of the dynamical tides, which, on the one hand, shrinks the orbit and thus imposes a stronger high-energy flux received by the planet and, on the other hand, leads to the angular momentum that is lost by the orbit of the planet being transferred to the star, which, therefore, rotates faster and has an increased radiation activity. The latter point, however, does not play a significant role here because during this phase, the X-ray and UV luminosities are independent of the stellar rotation because they are saturated (see Fig. 4). However the duration of the phase during which saturation occurs depends on the initial rotation. This duration is around 170 and 570 Myr if the initial rotation is respectively equal to 3.5 and 70 km s⁻¹.

In cases where saturation was not achieved, we further investigate whether the change of the orbital angular momentum is sufficient to significantly spin-up the star. The spin angular momentum of the star is on the order of ~ 1.65 × 10⁵⁰ g cm² s⁻¹, while the initial orbital angular momentum of the planet is less than 10% of this value. Therefore, even if the totality of the orbital angular momentum were to be transferred to the star, it would produce an increase of the surface velocity by around 10% (if the stellar structure is kept constant). In reality, only a fraction of the initial orbital angular momentum can be transferred to the star since, in our computation, the planet is completely evaporated at some distance from the host star. This causes the surface velocity of the star (at the time when the planet disappears due to the evaporation process) to be only 5% larger than the velocity it would have had in the absence of any planet. Outside the saturation regime, such a stellar spin-up would cause an increase of the XUV flux by about 15%, which would be non-negligible.

For the cases where the planets survive, we find that photo-evaporation decreases their total mass by a fraction between 11 and 16% at the end of the MS phase for the initial distances equal to 0.05 and 0.04 au, respectively. Most of the mass loss will occur during the first 300 Myr, when the Rossby number is below 0.196 (see the blue continuous line in the left panel ofFig. 4). As soon as the Rossby number becomes larger than this limit, we can see that the XUV flux decreases rapidly. This occurs in our model when the surface velocity is reduced below about 16 km s⁻¹ (see the right panel of Fig. 4). We note that starting with the same angular momentum content, a model where differential rotation occurs shows a very different evolution for the surface velocity than a model rotating as a solid body. The differences are small during the contracting phase towards the ZAMS (corresponding to the peak, occurring around 30 Myr), but they are important during the MS phase. The solid body rotating model shows larger surface velocities⁴ that keep the Rossby number below 0.196 for a much longer period (compare the continuous blue line for the solid body rotating model with the dotted one for the differentially rotating one). This illustrates well the fact that starting from the same initial angular momentum content, depending on the way angular momentum is transported inside the star, very different results can be obtained.

Fig. 4 Evolution of the XUV irradiation flux and of the Rossby number. Left panel: red continuous line shows the evolution as a functionof time of the XUV flux on the planet beginning its evolution at a distance of 0.05 au. The blue continuous line shows the evolution of the Rossby number obtained from the solid body stellar rotation model. The dotted blue line is the Rossby number obtained from the differentially rotating model. Below the horizontal line at 0.196, the X-ray luminosity is a constant fraction (0.00074) of the bolometric luminosity. When the curve is above the horizontal line, the X-ray luminosity decreases rapidly as the Rossby number increases (see text). Right panel: evolution of the surface velocity of the 1 M⊙ model. The red line shows the surface velocity assuming solid body rotation. The blue line is the surface velocity obtained from the model used in Rao et al. (2018) where differential rotation occurs due to the internal angular momentum redistribution resulting from shear instability and meridional currents.

3.3 Other planet masses

Figure 3 presents time evolutions of the planetary orbit and mass during the first 40 Myr, as well as the variations of the planetary radius and of the evaporation rate as a function of the orbital period. We considered here a rotation of the star equal to about 40 km s⁻¹, 40 Myr before the star reaches the ZAMS.

As has been found by many previous works, evaporation will impact mostly low mass planets (see e.g., references in the recent review by Owen 2019). From the present simulations, we distinguish two planetary mass ranges: the intermediate and high-mass planets, showing different behaviours under the action of tides and evaporation⁵ :

Intermediate-mass planets (~ 0.02 M_Jup < m_pl < ~2.5 M_Jup)

In this domain, both tides and evaporation have significant effects. Interestingly, this is also the mass domain where the evaporation rate goes through a maximum. Indeed, as shown in the right panel of Fig. B.2, the evaporation rate becomes maximum for a mass around 124 M_Earth, or 0.39 M_Jup, where the mass radius relation (see the left panel of Fig. B.1) changes slopes. For the 0.1 M_Jup planet, we see that regardless of whether the orbit shrinks or expands due to tides, as long as the initial distance is below 0.05 au, the planet is completely evaporated on timescales shorter than a few 100 Myr. A main difference between the 1 M_Jup planet and the 0.1 M_Jup planet is the fact that this more massive planet escapes complete evaporation if its initial distance is larger than around 0.034 au. As seen previously, it is only when the planetary orbit shrinks that we have a complete evaporation. The cases where the orbits expand display very little mass loss and very small changes of the radius even after a few Gyr. The sharp drop in the evaporation rate shown in the panel of the last column in the second row of Fig. 3 occurs at the transition between the saturated and unsaturated regime for computing the X-ray flux. This occurs at an age around 300 Myr (see the left panel of Fig. 4). The evaporation rate drops to a value of 10⁹ g s⁻¹, which is 1.7 × 10⁻¹¹ M_Jup yr⁻¹. The case of the 2 M_Jup planet shows many qualitative similarities to the case of the 1 M_Jup planet. The evaporation rates are smaller since the gravitational potential well is deeper for more massive planets.

High-mass planets (m_pl > ~2.5 M_Jup)

We recall here that above 2.5 M_Jup, the present evaporation laws cannot be applied, in principle. However, these laws lead to very small evaporation rates, as can be seen lookingat the panels of the fourth row in Fig. 3. Hence, the evolution in these cases is similar to models without evaporation and this is the domain where tides dominate. Planets with an initial distance that is below around 0.03 au havetheir orbits strongly affected by dynamical tides and these planets are engulfed by the star. Such an engulfment process produces a significant spin-up of the star. We can see in Fig. 5, that for the 30 M_Jup mass planetstarting its evolution at 0.03 au, the surface velocity is increased by a factor of two⁶.

In Fig. 6, we compare the fate of planets during the PMS phase, with and without evaporation. We see, as expected, that the region that is the most affected by the evaporation process is the low-intermediate planetary mass range (planets with a mass below about 0.04 M_Jup). In this massdomain, a main consequence of the evaporation process is that it significantly expands the domain of initial conditions leading to the disappearance of the planet (or at least of the planetary atmosphere) as compared to the corresponding domain obtained by considering only tides and no evaporation. Typically, for a planet with an initial mass of 0.02 M_Jup, with no evaporation process, the planet avoids engulfment if its initial distance is larger than about 0.025 au. With evaporation, the planet avoids evaporation if the initial distance is above 0.035 au. In Sect. 5, a comparison of the present results (of planets or at least of planets with an atmosphere) is made with those obtained by Kurokawa & Nakamoto (2014).

3.4 Evolution beyond the PMS phase

In Fig. 7, we show how the mass of the planets changes as a function of time depending on the initial distance between theplanet and the star and for two different initial-mass planets. A much larger range of initial distances is considered, compared to Figs. 3 and 6, as well as a much larger time domain (the PMS phase ends approximately at an age equal to 40 Myr). The domain of disappearance of the planet (the limit between thegrey and the purple area) increases significantly with time for planets with an initial mass equal to 0.1 M_Jup planets. This is, of course, a consequence of the fact that more time allows for more evaporation, which implies that the planet, to avoid becoming completely evaporated by a given time, has to begin its evolution at a greater initial orbital distance.

We see that for this mass, the limits between the diverse colored areas increases significantly up to around 300 Myr. After this period, as mentioned previously, the X-ray flux decreases, and hence the evaporation rate decreases as well, causing the limits to vary less rapidly over time than in the preceding period.

In the case of the 1 M_Jup mass planet, all the cases beginning their evolution at a distance below 0.03 au are completely evaporated at more or less the same age (see the second row, second column panel in Fig. 3). Interestingly, for this specific planetary mass, complete evaporation occurs only during the PMS phase – otherwise, it does not occur at all. This is reflected by the fact that the upper limit of the grey zone is a nearly horizontal line. If a planet begins its evolution at a distance a bit larger than about 0.03 au, then it will escape complete destruction during the whole MS phase. We see that starting at an orbital distance above 0.03 au will make the time during which tides can shrink the orbit very brief. Indeed, at this distance, the orbit rapidly crosses the corotation radius. After this crossing, the orbital distance expands, thus saving the planet from the loss of its atmosphere as it spirals away from the star.

Comparing the two panels of Fig. 7 shows that increasing the initial mass of the planet reduces the zones where the planetary mass changes. This simply reflects the fact that a more massive planet, in order to be evaporated at a given age, must start its evolution nearer to its host star. This effect is particularly significant when going from 0.1 to 1 M_Jup mass, that is, through the mass range where the mass-radius relation for the planet changes slope. Beyond a mass of around 0.4 M_Jup, the gravitational potential increases more rapidly with the mass than in the lower-mass domain. This makes the planets more resistant to an evaporation process and, thus, significantly restricts the conditions for losing a given amount of mass.

Fig. 5 Time evolution of the stellar surface rotation velocity when the orbits of planets with 10 (blue dotted curves) and 30 MJup (red dotted curves) shrink, compared to the surface velocity evolution of the star without any planet (continuous green curve). Three dotted curves are shown per planet for three initial distances from the star: the lower curves correspond to 0.01 au, and these follow the same path whatever the mass of the planet, coinciding with the green curve. The middle and upper curves correspond to initial distances equal to 0.02 and 0.03 au, respectively.

Fig. 6 Fate of planets of different masses (in Jupiter masses), starting their evolution at various distances (in au) from their host star during the PMS phase. The left and right panels show respectively, the fate of planets without accounting for the evaporation process and with account of the evaporation process. Planets in the salmon-colored region are engulfed or evaporated. For low-mass planets, evaporation usually occurs before engulfment, and vice-versa for high-mass planets. Therefore, the transition of the fate happens somewhere in the intermediate-mass region. Planets in sky-blue regions and in turquoise regions have their orbits widened by tides.

Fig. 7 Evolution as a function of time of the minimum initial distance between the planet and its host star above which the planet survives (limit between the grey and pink areas), the planet loses less than 20%of its initial mass (limit between the purple and salmon area), and the planet loses less than 3% of its initial mass (limit between the salmon and light blue areas). The different panels show these limits for different values of the initial mass of the planet. The initial rotation of the star is 40 km−1 at 30 Myr before the ZAMS.

3.5 Impact of the initial stellar rotation

In Fig. 8, we compare, for three different initial rotations of the star, how the minimum distance for an evaporation or an engulfment during the PMS phase (or both) varies as a function of the planetary mass.

We note that for planetary masses below about 0.1 M_Jup, changing rotation has little effect. This finding is in line with two facts: first, for this planetary range, the evaporation process dominates. Second, the XUV flux is saturated (see the left panel of Fig. 4), and thus does not depend on rotation.

In the high-mass planetary range (for masses above about 0.1 M_Jup), the initial stellar rotation has a significant impact. As a numerical example, a 1 M_Jup planet orbiting a 1 M_⊙ star with an initial rotation of 4 km s⁻¹, 30 Myr before reaching the ZAMS, can survive if its initial distance is larger than about 0.01 au. When the initial stellar rotation is 40 km s⁻¹, the planet should have an initial distance larger than 0.03 au. This is a consequence of the dependence of the dynamical tides with the rotation rate of the star.

If we compare the situation where the planet orbits a moderately fast rotating star (40 km s⁻¹) to the one where it orbits a fast-rotating model (80 km s⁻¹, in comparing the middle and lower panels of Fig. 8), we observe for the faster-rotating star, a decrease of the minimum distance for engulfment in the planetary mass range above 0.1 M_Jup (engulfment due to tides dominates over evaporation in the high-mass planet range). This may appear surprising at first because faster rotation implies stronger dynamical tides and, thus, more efficient shrinking of the orbit. But at the same time, faster rotation also implies that the corotation radius is pushed inward, reducing thus the zone where dynamical tides shrink the orbit and favor engulfment.

4 Different types of evolution

Very generally, close-in planets show different types of evolution depending on the relative magnitudes of three timescales. They are the evaporation timescale, τ_evap = M_pl∕Ṁ_pl, the tidal timescale, ${\tau }_{\text{tides}}={(a/\dot{a})}_{\text{ti}}$ with a being the orbital distance, and the stellar evolutionary timescale τ_star, which is more rapid during the PMS phase (mainly the Kelvin–Helmholtz timescale) than during the MS phase (the nuclear timescale, although as explained below for the specific question addressed here, the magnetic braking timescale is certainly more relevant). We briefly summarise the different types of evolution with emphasis on how stellar rotation and feedback affect them since these effects are the focus of the present study. The different cases described below apply to different initial planet mass at orbital distance domains. We order their description beginning with the domain containing the most massive, distant planets and ending with the least massive, closest ones.

Fig. 8 Evolution as a function of the initial mass of the planet of the minimum initial distance between the planet and its host star above which the planet survive the PMS phase (limit between the grey and pink areas), the planet loses less than 20% of its initial mass (limit between the purple and salmon area), and the planet loses less than 3% of its initial mass (limit between the salmon and light blue areas). Upper panel: initial rotation of the star, taken at an age equal to 30 Myr before the ZAMS is 4 km s−1. Middle panel: same as upper panel, but for an initial rotation of the star equal to 40 km s−1. Lower panel: same as upper panel, but for an initial rotation of the star equal to 80 km s−1.

4.1 Constant orbit and mass evolution (τ_evap & τ_tides ≫ τ_star)

This evolution occurs for distant and high-mass planets, so that neither the orbit nor the mass of the planet are significantly affected by tides and evaporation. The evolution of the planet in such cases is not affected by the rotation of the star either, although the minimum mass (given an initial distance) or the minimum distance (given an initial planet mass) for such evolution to occur is shifted to larger values when the stellar rotation increases.

4.2 Constant mass evolution (τ_evap ≫ τ_star > τ_tides)

This occurs for sufficiently large mass planets having an initial orbital distance below the corotation radius, so that tides shrink the orbit. The upper limit of the initial orbital distance for this type of evolution to occur is a bit below the initial corotation radius, where the orbit plunges just fast enough such that it always remains below the corotation radius, which itself also shrinks with time during the PMS phase. Thus, the orbit keeps plunging till it reaches the minimum orbital distance above which dynamical tides become active, a_min. In other cases, if the planet begins its evolution below a_min, then for a while at first, the orbit is not changed. At a given point, however, the orbital distance will become equal to a_min, because thestar contracts during the PMS phase, so its rotation increases and $a_{\text{min}}\propto {\Omega }_{\text{star}}^2$ decreases. From this stage on, in all cases, the orbital distance will remain in the vicinity of this limit (see, as general examples of this effect, the orbits beginning at 0.02, 0.03 au in the first panels of the second, third and fourth rows, in Fig. 3⁷). This happens because every time the orbital distance drops a bit below a_min due to the tides, the dynamical tides switch off and the orbital distance momentarily stops changing until the a_min limit reduces and crosses the orbit, switching the dynamical tides back on, and the process thus continues.

Engulfment will occur for these planets during the PMS phase. In those cases, tides may spin-up the star significantly and feedback may be important enough to be accounted for. The domain for this type of evolution depends on rotation in two ways: first, through the strength of the dynamical tides that increases when rotation increases and second, through the determination of the upper (corotation radius) and lower distances (a_min) where dynamical tides can shrink the orbit. Both these distances are shifted to closer-in distances when rotation increases. Inorder to describe the evolution in this case, we need to have this region sufficiently close to the star so that the tidal timescale is short enough. In general, increasing the initial rotation will increase the domain in the plane ofinitial mass versus the initial distance where this type of evolution occurs. However, beyond some limiting rotation, we expect to have a reduction of this domain due to the fact that tides will still be high enough to produce an engulfment during the PMS phase, but the domain where the tides are active becomes smaller (see Fig. 8). Here, feedback effects can also be important since, as the orbit shrinks, the star will spin-up and thus, the strength of tides increases, while the domain where tides are active shrinks and shifts closer to the star.

Fig. 9 Data points are masses (Msin i) and semi-major axis for confirmed planets orbiting stars with masses between 0.9 and 1.1 M⊙. The light-grey shaded area indicates the Neptunian desert zone after 40 Myr as predicted by our models for our slow-rotating case (Vini = 4 km s−1). The light pink zone shows the enlargement of that Neptunian desert region just after 40 Myr as predicted by our models for our moderate rotating case (Vini = 40 km s−1). The magenta shaded area shows the enlargement of that zone after 300 Myr. For comparison, also plotted are the lines corresponding to the minimum survival masses of planets with core masses of 0.03 MJup (solid line), 0.06 MJup (dashed line) and 0.09 MJup (fuzzy line) at an age of 10 Gyr obtained by Kurokawa & Nakamoto (2014). The corresponding core masses are shown by horizontal lines.

4.3 Orbit and mass evolution (τ_evap & τ_tides < τ_star)

This evolution characterizes the domain of sub-Neptune and hot Jupiter planets. In these cases, evaporation can occur either during the PMS phase or during the early MS phase when the star still has a high enough rotation to emit a saturated XUV luminosity. This is an important domain because it impacts a large part of the survival zone shown in Fig. 9.

If the planet begins its evolution at an initial orbital distance which is above the upper limit mentioned in the previous subsection, then the orbital distance will eventually increase due to tides after the corotation radius shrinks below it. Therefore, the star willspin-down, which will not impact the evaporation efficiency during the PMS phase since the star contracts rapidly during this time, keeping its XUV photon emission saturated even if some angular momentum is lost at the advantage of the expanding planetary orbit. The irradiation flux does, however, decrease because the orbital distance increases, which would therefore also increase the evaporation timescale. If they are not too massive and not too distant from their host stars, these planets can still lose significant amounts of their initial mass over periods of a few 100 Myr. Examples of such an evolution in the present simulations are represented by the cases of the 0.1 M_Jup planets withinitial orbital distances of 0.03 and 0.04 au that are evaporated around 100 Myr (calculated for an extrapolation of the cases shown in the first row of Fig. 3). For higher values of the initial rotation of the star, the PMS phase described above will not be affected much. However, the period during which the star emits high energy photons after the PMS phase will grow. This decreases the evaporation timescale. Not only is the initial rotation of the star important, but also the efficiency of the internal angular momentum transport mechanism. A less efficient transport mechanism than the one considered here may reduce, by a factor of three, the duration of the XUV radiation emissions.

If, instead, the planet begins its evolution below the orbital distance upper limit mentioned in the previous subsection, then the orbit shrinks. The shrinking of the orbit causes the planet to experience a higher irradiation flux, and thus, both timescales, τ_tides and τ_evap, decrease. Inour computations, we find that the evaporation effect dominates over tides in these cases, obliterating the planet before engulfment. This is seen in, for instance, the cases of the 1 M_Jup mass planet with initial orbital distances of 0.02 au and 0.03 au, in the third row of Fig. 3.

Changing the initial rotation has modest effects on this type of evolution for the reasons already mentioned above. Feedback effects of the tides on the stellar rotation are modest because of the fact that this type of evolution involves planetary masses and distances that provide modest amounts of angular momentum and only part of it is transmitted to the star through tides because of the evaporation process.

After the early MS phase that lasts a few 100 Myr for a star like the Sun, the planets that have survived no longer evolve much, at least for those whose mass is below the critical mass for activating tidal dissipation in radiative zones. It is only when the star becomes a red giant that new important changes can occur, mainly due to equilibrium tides. Thus, the surviving planet here will follow a constant orbit and mass evolution during the MS phase.

4.4 Mass evolution (τ_evap ≪ τ_tides < τ_star)

This evolution occurs for planets that are the most close-in and have the lowest mass in the ranges considered in our work. The planet loses its mass before tides can modify its orbit. For instance, in the present models, we find that a 0.1 M_Jup planet with an initial distance of 0.01 au would have such an evolution (see first row of Fig. 3). Stellar rotation has little impact on this type of evolution. Even starting from a very small initial rotation, since the star is contracting during the PMS phase, it will rapidly spin-up and emit high-energy photons at the saturation limit. Since tides are not important, rotation will also have no impact on their action.

5 Comparisons with the observations

5.1 Neptunian desert

In Fig. 9, the masses (Msin i) and distances of confirmed exoplanets orbiting solar mass stars are shown⁸. We overplotted shaded areas indicating the regions that, according to our simulations, should be devoid of planets either because of evaporation or engulfment. We have focused here on the planetary mass domain between 0.03 and 1 M_Jup (i.e., between 10 and 310 Earth masses), which lies well inside the domain where the evaporation law adopted here has been applied in previous works (Salz et al. (2016).

The grey- and salmon-shaded areas show the Neptunian desert after 40 Myr when the initial stellar rotation is 4 and 40 km s⁻¹ respectively.Increasing the initial stellar rotation enlarges the desert zone, especially for the most massive planets of the mass domain considered here. This is mainly an effect linked to the fact that when the rotation of the star increases, tides are stronger (given a fixed initial distance). For those planets below the corotation radius, this will favor planetary engulfment. For the lower mass end considered here, increasing rotation also enhances the tides which, for the orbit shrinking planets, favors evaporation.

When systems older than about 40 Myr are considered (here, up to 300 Myr and for an initial rotation of 40 km s⁻¹), the desert domain is enlarged, mainly because it gives more time for the planet to be irradiated by the high energy photons of the star. Changing rotation has an impact on this region, since it modifies the duration of the period during which the XUV flux is saturated.

A very interesting feature to note is that the enlargement of the desert zone when the age increases is much larger for the 0.1 M_Jup planet than for the 1 M_Jup one. This is because for a low-mass planet, the evaporation timescale at a given distance is much shorter than for a high-mass one. This effect imposes thus that the right frontier of the magenta zone goes towards higher initial distances when the mass of the planet decreases. We see that overall, the zone devoid of planets predicted by the present simulations is reasonable with respect to the observed distribution of planets in this plane. Of course, much more refined comparisons are needed to incorporate more realistic planetary models.

For comparison, the minimum survival masses of planets with different core masses observed at an age of 10 Gyr, as predicted by Kurokawa & Nakamoto (2014), are shown. These authors considered more sophisticated interior planet models, but they did not consider any effects of tides. This is mainly the reason why their survival limit is at very small distances for the 1 M_Jup model. They adopted the XUV evolution models given in Ribas et al. (2005), based on observations of solar analog stars, thus not accounting for the impact of different initial rotations and for exchanges of angular momentum between the planet orbit and the stellar rotation. We see that our results significantly enlarge the zones devoid of planets. On the one hand, this is a result of our accounting for the impacts of tides here. On the other hand, we may here somewhat overestimate the desert zone because we do not allow for the existence of rocky cores.

We note that a more comprehensive interpretation of the observed Neptunian desert also requires accounting for high-eccentricity orbital migration (Owen 2019). Although the present simulations do not account for eccentric planet orbits, we show here that simply accounting for tides, evaporation, and stellar evolution provides a good approximation for explaining the observed results.

5.2 Orbital period upper limit for bare-core planets

From the right panel of Fig. 6, we can deduce the trend that planets with the smallest initial mass in the range we consider canbegin their evolution at the farthest initial orbital distance from the star, such that at the end of the evolution, that is, at a given final age of the system, the planets lose their entire atmosphere due to evaporation, while still being as far away from the star as possible. In other words, for star–planet systems of a given final age, computations of the final orbital period of the planets beginning their evolution with the smallest initial mass considered in our study, which just lose their atmosphere due to evaporation at that age, would be an upper limit on the orbital period as predicted by the present simulations, above which small-radii or bare-core planets cannot be observed.

In Fig. 10, we make this computation, and find that the upper limit of the orbital period for observing bare-core planets around solar-type stars of a few Gyr age, is around 17 (Earth) days. This number is in surprisingly good agreement with the locus of the small-radii or bare-core population constituting in part the “radius valley" feature, namely, the bimodal distribution of exoplanets observed by the Kepler mission, in a scatter plot of planet radius versus the orbital period (see Van Eylen et al. 2018). We find that both in observations, as well as in planet population simulations accounting for more detailed planet physics than we do (Owen & Wu 2017), the bulk of the small-radii or bare-core population lies just below the upper limit predicted by our code, indicating this limit to mark, as we would expect in principle, a boundary of the radius valley feature.

However, since we do not model the interior structure of planets in the present work, we cannot say anything about the radii of the bare-core planets left over after complete evaporation of their atmospheres. Hence, we can only find numerical limits/boundaries along the orbital period (horizontal) axis, but not along the planet radius (vertical) axis of the planet radius versus orbital period scatter plot.

It would appear that the small fraction of observed bare-core planets that do lie above our predicted orbital period upper limit must have followed an evolution or initial conditions different from the ones covered by the present simulations, implying, unsurprisingly, that more physics is involved in star–planet systems than is considered here. On the other hand, the very good agreements, taken at face-value, imply that many star–planet systems in the observationsmay have followed an evolution similar to the ones predicted by our simple, albeit comprehensive models.

Fig. 10 Evolution of the planet masses as a function of the age. The initial mass of the planets is 0.02 MJup, being the smallest in the mass range considered in our study. The stellar model is for a 1 M⊙ star beginning its evolution with a surface velocity of 40 km s−1, 30 Myr before the ZAMS. Each curve is labeled by the orbital period of the planet at the end of the evolution (see text) in units of Earth days. From the shortest to the longest orbital period, the initial distance has the following respective values in au: 0.0100, 0.0150, 0.0250, 0.0450, 0.0650, 0.0850, 0.1100, 0.1294.

6 Advantages and limitations of the present approach

The main advantages of the present consistent treatment are the possibilities to study the impact of different initial rotations on the evolution of the planetary mass and orbit, of accounting for the feedback of these changes on the stellar rotation, and thus on the tides and evaporation rates.

A main present limitation is the rather schematic interior model used for the planets. In this work, we use a mass-radius relation that fits direct observational exoplanet data and numerical simulations of planetary formation (Bashi et al. 2017), therefore serving as a valid first approximation. What is shown here is that, at least in some planetary mass domains, the impact of changing the initial rotation of the star can have effects as large, or even larger than those observed when different planetary interior models are considered (see Fig. 9).

We did not consider multiple planets around the star and, therefore, any interactions between them, nor the impacts of possible magnetic interactions between the star and planet. Neither did we account for eccentric or inclined planetary orbits, consistently assuming a quasi-circular orbit of the planet in the rotation plane of the star.

The limits of the survival zone depend on when the planet arrived at the point defined here as the initial orbital distance to the star. We recall that here we consider a starting point that corresponds to a time of about 40 Myr before the ZAMS. Would we have considered that the protoplanetary disk dissipates at later time, then it would impact the strength of the tides, making them weaker. Indeed, the extent of the convective envelope and the radius of the star decrease as a function of time during the PMS phase. The planet may also become a close-in planet while beginning its evolution at much larger distances due to migration or dynamical interactions. A fraction of close-in planets could indeed have migrated long after their formation via Kozai migration, thus delaying their atmospheric evolution and evaporation by several Gyr (Attia et al. 2021).

Some uncertainties pertain the XUV stellar luminosity. Recently, King & Wheatley (2021) argued that the decline in the EUV emission is much slower than for the X-rays and thus might dominate the irradiation of the planetary atmosphere after the saturation phase.

We neglected the mass loss by the star and the possible drag forces acting on the planetary orbit. These aspects do not appear to be very relevant, however. In Fig. 2, the models from Rao et al. (2018) were computed while accounting for these drag forces, and as we can see, they compare extremely well with the present models that neglected these forces.

In a previous section, we briefly mentioned the fact that not only the initial rotation has an impact on the findings, but also the efficiency of the internal angular momentum transport inside the star. This is certainly a point worthwhile to be studied in a future work.

We focused in the present study on solar mass star models at solar metallicity. Different types of stars have different radiative and convective regions, which change the way tides act. They also exhibit different XUV radiation emission, modifying the photo-evaporation process. Also, the evolution of their rotation can be different depending on the importance of the outer convective zone. We plan an exploration of these effects in a future work.

We consider here two types of tides, namely, the equilibrium and dynamical tides in convective regions. The expressions for these tides are still subject to discussion (see e.g., the discussion in Barker 2020). Also, the dissipation of tides in radiative zones due to internal gravity waves has not been accounted for in the present work. However, according to the work cited above (Barker 2020), this mechanism is efficient only for planetary masses above a certain critical limit. For planets orbiting a 1 M_⊙ star, this critical limit is larger than 10 M_Jup for ages below about 3.5 Gyr (see Fig. 8 in Barker 2020). The results presented above concern lower mass planets in an age domain below the one indicated above. Thus, they should not be affected by the neglect of the tides in radiative zones. In Fig. 7, we show an evolution over a longer period, until an age of 8.4 Gyr, which corresponds to the end of the core H-burning phase for a 1 M_⊙ star. The minimum planetary mass at thatage is 0.4 M_Jup according to Barker (2020), hence being above 0.1 M_Jup.

There are also the uncertainties pertaining to the evaporation rate. Different choices for the heating efficiency ϵ, the factor β describing the inflation of the planet’s outer layers, the distance to the Roche radius etc. will change the results. Here, we have adopted the given set of values because we wanted to study the impacts of both tides and the evaporation process, as well as the impact of the initial rotation of the star. Of course, different evaporation laws may bring significant changes to the results.

7 Conclusions

As mentioned in the introduction, very few studies follow, in a consistent way, the change of the planetary orbit with that of the stellar rotation, despite the fact that these quantities are tightly connected and impact each other as well as the planetary irradiation and, thus, its evaporation rate.

In the present work, using complete stellar models for describing the evolution of the structure, assuming solid body rotation for the star (which appears as a very reasonable assumption, see e.g. García et al. 2007; Benomar et al. 2015; Eggenberger et al. 2019), we explore the consequences of these two effects in a 4D-space domain whose axes are the planetary mass (in the range between 0.02 up to 2.5 M_Jup mass), the initial orbital distance (from 0.01 to 1 au), the stellar rotation (from 4 to 80 km⁻¹, 40 Myr before the ZAMS), and the age of the system (from the PMS until the end of the early MS phase).

We checked the continuity of the present results with those obtained using a previous version of our code where the evaporation process was not accounted for. Despite the simple models we use for the planet’s interiors, comparisons with observations shows that at least for sub-Neptune planets, the results obtained here are reasonable. This gives us the confidence that our study of the impact of the initial rotation of the star and of the interactions between tides, evaporation, and stellar rotations have some validity.

We distinguish two effects of rotation: those linked to different values of the initial rotation of the star and those arising from the change of the stellar rotation due star–planet interactions. Changing the initial stellar rotation has the greatest effects on two families of evolution: the constant-mass-evolution, and the orbit-and-mass-evolution (Sect. 4.3). In both cases, overall, increasing the initial rotation tends to increase the domain in initial planetary masses and orbital distances where this type of evolution occurs. As mentioned earlier, however, beyond some rotation limit, the domain may decrease due to the dependence on rotation of the limits where dynamical tides are active. The effects of changing the initial rotation are large on the populations of planets expected just after the PMS phase or after the early MS phase. At larger ages, the efficient wind magnetic braking makes the system forget the initial rotation.

Interaction between tides and stellar rotation are expected to occur in the case of constant-mass-evolution and to a lesser extent in the case of orbit-and-mass-evolution families. This may produce some significant increase of the surface velocity of the star, as has also been shown by Gallet et al. (2018). The spin-up of the star due to an engulfment may have a significant impact on other, less massive planets that are also present in the system, which, after the engulfment of a massive planet, may experience stronger tides. This is a point of particular interest that ought to be studied within the framework of the evolution of planetary systems.

We have seen that the spin-up of the star due to tides has little effect in terms of increase of the XUV luminosity since, in the present simulations, when irradiation is efficient, it is saturated. The effect of tides (when shrinking the orbit) implies an increase of the irradiation because the star–planet distance becomes smaller and, thus, the planet lifetime is shortened.

Although we did not study this effect here, we expect that a change in the efficiency of the angular momentum transport inside the star may significantly change the results for those planets whose evaporation occurs during the early MS phase. Adopting a less efficient angular momentum transport than the one in this study may shorten the duration of the phase during which the star emits a high XUV radiation luminosity. One general conclusion of the present work is that the huge efforts that have been carried out thus far to better describe the planet physics, evaporation process (etc.) should also be done within the framework of models where the physics of the star is better described.

Acknowledgements

We acknowledge the anonymous referee’s critical comments which motivated us to significantly revise and improve the manuscript. This research has made use of the Extrasolar Planet Encyclopedia, http://exoplanet.eu. This project has been supported by the Swiss National Science Foundation grant 200020-172505, and has been carried out in part within the frame of the National Centre for Competence in Research (NCCR) PlanetS.

Appendix A Stellar models

Figure A.1 shows the evolution of the rotation period and of the surface velocity (at the equator) of a 1 M_⊙ stellar model from the PMS phase to the early MS phase (only the first 40 Myr are shown). The star is here assumed to rotate as a solid body. During the PMS phase, the star is contracting rapidly. This decreases the rotation period (increasing the surface equatorial velocity). The position of the ZAMS, that may be defined here as the time at which the minimum radius is reached, corresponds to the peak in velocity shown on the right panel. During the MS phase, the surface velocity decreases due to the wind magnetic braking. Whatever the rotation of the star on the ZAMS, after 4.6 Gyr, the surface rotation would be less than 2 km s⁻¹, right in the range of values observed for the Sun.

Fig. A.1 Evolution of the rotation period (left panel) and of the surface velocity (right panel) for an isolated 1 M⊙ stellar model during the PMS phase and the very early MS phase for different values of the initial rotation. The green line shows respectively the critical period (left panel) and the critical velocity (right panel). We note that the bottom line in the left panel corresponds to the upper line in the right panel.

Appendix B Planetary models

The variation of the planetary radius as a function of the planetary mass, from Bashi et al. (2017), is shown in the left panel of Fig. B.1 (see the green curve). The piece-wise description provides a distinction between “small" and “large” planets. The transition mass separating the small planets from the large planets is m_t = 124 m_Earth and the transition radius is r_t = 12.1 r_Earth.

When the planet is very close to the star, its atmosphere may be inflated due to heating caused by the XUV irradiation from the star. In that case, the radius of the planet has to be corrected. The inflated radius is given by βr_pl, where r_pl is the original planet radius, and β taken as in Salz et al. (2016). The expression for β is valid only in a limited domain of gravitational potential, between 10^12.2 and 10^13.6 erg g⁻¹. Using the mass-radius relation by Bashi et al. (2017), this corresponds to masses between 6 and 800 Earth masses.

The left panel of Fig. B.1 shows the variation of the original radius (green solid line) and of the inflated radius (red solid lines) as a function of the mass of the planet for different irradiation values (see caption).

If a blob of planetary matter receives sufficient energy to enable it to reach the limit of the Roche lobe between the planet and the star, then it will escape from the effective gravitational potential well of the planet orbiting its star, and will be lost. The Roche lobe boundary, R_Roche, is taken from Erkaev et al. (2007). Its variation as a function of the planetary mass for different orbital distances is shown in the right panel of Fig. B.1. We see that in general, the Roche lobe is larger than the (non-inflated) planetary radius, except for cases when the orbital distance becomes equal to about 0.01 au or less. In the case XUV absorption radii are considered, orbital distances larger than about 0.02 au are needed to avoid that the Roche lobe becomes smaller than the planet radius. The present treatment is not sophisticated enough to follow the evolution in a very rigorous way when the planetary radius exceeds the Roche lobe, because this effect involves a mass loss mechanism different from the one triggered by the XUV luminosity of the star. It is likely that in such cases, the planet will lose its whole atmosphere in a very short timescale, and for gaseous planets, nearly complete evaporation may occur. This is the scenario considered to occur in the present work.

In the left panel of Fig. B.2, the energy required for a unit mass to reach the Roche lobe, as given by Erkaev et al. (2007), is shown for different values of the mass of the planet and for different distances to the star (distances decreases from the green line to the lowest red curve).

In the right panel of Fig. B.2, the timescale needed for the evaporation of 1 Earth mass according to the evaporation law given by Eq. (2) is shown as a function of the mass of the planet and for different XUV irradiation fluxes. The most striking feature is the fact that the evaporation timescale has a minimum where the planetary mass-radius relation (shown in the left panel of Fig. B.1) changes slope. For masses below the mass corresponding to that transition, the increase of the XUV absorption radius (βr_pl) overcomes the relatively gentle increase with the planetary mass of K|ϕ_G|. For masses above this mass limit, βr_pl decreases, while ΔΦ increases, therefore the evaporation timescale increases. Hence, we see that the mass-radius relation plays a key role here. We note that the masses of the planets that show the largest evaporation rates are masses around that transition (here ~120 M_Earth,  ~0.38M_Jup).

The process of planetary evaporation during the PMS phase has a rather limited effect when the distance is equal or above 0.1 au. Indeed, the PMS phase lasts 30–40 Myr, while the time for evaporating 1 earth mass from a 100 M_Earth is around 250 Myr. When the distance is decreased by a factor of ten, the situation is very different. Here we have that whatever the irradiation flux provided superior to 10^4.5 erg s⁻¹ cm⁻², the low-mass planets with a few M_Earth will be significantly affected by the evaporation process.

Fig. B.1 Variation of the planetary radius and of the Roche lobe as a function of the planetary mass. Left panel: variation of the planetary radius as a function of the planetary mass (green curve; see text). The red curves indicate the effective radius accounting for the fact the atmosphere of the planet can be inflated. From bottom to top, the curves corresponds to values of the XUV flux between 104.5 and 105.5 in steps of 0.2 dex. These XUV fluxes corresponds to XUV luminosities equal to 10−3 L⊙ and orbital distances between about 0.05 and 0.2 au. Right panel: variation of the Roche lobe as a function of the planetary mass (red curves) compared to the planetary radius (thick green curve). The red curves corresponds to distances between the planet and the star equal to 0.01, 0.02, 0.03, 0.04, 0.05, and 0.1,au, from bottom to top. The mass of the star is 1 M⊙.

Fig. B.2 Left panel: variation of the energy needed to bring a unit-mass planetary material up to the Roche lobe as a function of the mass of the planet and of the orbital distance. The green line corresponds to the case when the planet is in isolation, showing the energy needed to bring a unit mass to infinity. The red lines corresponds to cases where the planet is sufficiently near the star so that the energy would be lesser than the former case, which is needed for bringing a unit mass only up to the Roche lobe between the planet and the star. From top to bottom, red lines correspond to orbital distances equal to 0.1, 0.05, 0.04, 0.03, 0.02, and 0.01 au. Right panel: variation of the initial evaporation timescale of planets of different masses irradiated by different XUV fluxes. The green line corresponds to a flux equal to 104.5 erg s−1. The lowest red curve corresponds to a flux equal to 105.5 erg s−1. The lines in between correspond to intermediate values, from 104.6 to 105.4 in steps of 0.1 dex. The distance between the planet and the star has been assumed to be 0.1 au for the continuous lines and 0.01 au for the dashed lines.

References

All Figures

Fig. 1 Left panel: evolution of the Rossby number for an isolated 1 M⊙ stellar model during the PMS phase and the very early MS phase for different values of the initial rotation. The initial rotation periods and velocities can be seen in Fig. A.1. The Rossby number corresponding to the slowest initial rotation is the upper curve. Right panel: evolution of the X-ray, EUV luminosities and of their sum for a slow (in blue) and a fast (in red) rotating 1 M⊙ model. The slow and fast rotations correspond, respectively, to the slowest and fastest model shown in Fig. A.1.

In the text

Fig. 2 Evolution as a function of time of the orbits of 1 MJup mass planetsaround a 1 M⊙ star, obtained using different approaches. The blue continuous lines are obtained in the present work by switching off the planetary evaporation process. The red dashed lines are obtained switching on the evaporation process. The orange dashed lines (only two cases shown) are the orbits obtained in Rao et al. (2018). The blue dashed and dotted lines show respectively, the corotation radius (acorot) and the critical distance below which dynamical tides are no longer active (amin). They are obtained for the orbit beginning its evolution at a distance equal to 0.034 au. The upper-limit of the pink zone corresponds to the radius of the star in au. Its lower limit shows the radius at the base of the stellar convective envelope. The blue zone indicates the size of the stellar radiative interior.

In the text

Fig. 3 Evolution for a given planetary mass and for different initial radii of the orbit of the orbital distance as a function of time (first column panel), of the planetary mass as a function of time (second column), of the planetary radius as a function of the orbital period (third column), and of the evaporation rate as a function of orbital period (fourth column). Identical colors in a set of row panels correspond to models with the same initial distance to the star (see the left panel in each row). We note that the panels of the fourth row in the second and third column have a different ordinate.

In the text

Fig. 4 Evolution of the XUV irradiation flux and of the Rossby number. Left panel: red continuous line shows the evolution as a functionof time of the XUV flux on the planet beginning its evolution at a distance of 0.05 au. The blue continuous line shows the evolution of the Rossby number obtained from the solid body stellar rotation model. The dotted blue line is the Rossby number obtained from the differentially rotating model. Below the horizontal line at 0.196, the X-ray luminosity is a constant fraction (0.00074) of the bolometric luminosity. When the curve is above the horizontal line, the X-ray luminosity decreases rapidly as the Rossby number increases (see text). Right panel: evolution of the surface velocity of the 1 M⊙ model. The red line shows the surface velocity assuming solid body rotation. The blue line is the surface velocity obtained from the model used in Rao et al. (2018) where differential rotation occurs due to the internal angular momentum redistribution resulting from shear instability and meridional currents.

In the text

Fig. 5 Time evolution of the stellar surface rotation velocity when the orbits of planets with 10 (blue dotted curves) and 30 MJup (red dotted curves) shrink, compared to the surface velocity evolution of the star without any planet (continuous green curve). Three dotted curves are shown per planet for three initial distances from the star: the lower curves correspond to 0.01 au, and these follow the same path whatever the mass of the planet, coinciding with the green curve. The middle and upper curves correspond to initial distances equal to 0.02 and 0.03 au, respectively.

In the text

Fig. 6 Fate of planets of different masses (in Jupiter masses), starting their evolution at various distances (in au) from their host star during the PMS phase. The left and right panels show respectively, the fate of planets without accounting for the evaporation process and with account of the evaporation process. Planets in the salmon-colored region are engulfed or evaporated. For low-mass planets, evaporation usually occurs before engulfment, and vice-versa for high-mass planets. Therefore, the transition of the fate happens somewhere in the intermediate-mass region. Planets in sky-blue regions and in turquoise regions have their orbits widened by tides.

In the text

Fig. 7 Evolution as a function of time of the minimum initial distance between the planet and its host star above which the planet survives (limit between the grey and pink areas), the planet loses less than 20%of its initial mass (limit between the purple and salmon area), and the planet loses less than 3% of its initial mass (limit between the salmon and light blue areas). The different panels show these limits for different values of the initial mass of the planet. The initial rotation of the star is 40 km−1 at 30 Myr before the ZAMS.

In the text

Fig. 8 Evolution as a function of the initial mass of the planet of the minimum initial distance between the planet and its host star above which the planet survive the PMS phase (limit between the grey and pink areas), the planet loses less than 20% of its initial mass (limit between the purple and salmon area), and the planet loses less than 3% of its initial mass (limit between the salmon and light blue areas). Upper panel: initial rotation of the star, taken at an age equal to 30 Myr before the ZAMS is 4 km s−1. Middle panel: same as upper panel, but for an initial rotation of the star equal to 40 km s−1. Lower panel: same as upper panel, but for an initial rotation of the star equal to 80 km s−1.

In the text

Fig. 9 Data points are masses (Msin i) and semi-major axis for confirmed planets orbiting stars with masses between 0.9 and 1.1 M⊙. The light-grey shaded area indicates the Neptunian desert zone after 40 Myr as predicted by our models for our slow-rotating case (Vini = 4 km s−1). The light pink zone shows the enlargement of that Neptunian desert region just after 40 Myr as predicted by our models for our moderate rotating case (Vini = 40 km s−1). The magenta shaded area shows the enlargement of that zone after 300 Myr. For comparison, also plotted are the lines corresponding to the minimum survival masses of planets with core masses of 0.03 MJup (solid line), 0.06 MJup (dashed line) and 0.09 MJup (fuzzy line) at an age of 10 Gyr obtained by Kurokawa & Nakamoto (2014). The corresponding core masses are shown by horizontal lines.

In the text

Fig. 10 Evolution of the planet masses as a function of the age. The initial mass of the planets is 0.02 MJup, being the smallest in the mass range considered in our study. The stellar model is for a 1 M⊙ star beginning its evolution with a surface velocity of 40 km s−1, 30 Myr before the ZAMS. Each curve is labeled by the orbital period of the planet at the end of the evolution (see text) in units of Earth days. From the shortest to the longest orbital period, the initial distance has the following respective values in au: 0.0100, 0.0150, 0.0250, 0.0450, 0.0650, 0.0850, 0.1100, 0.1294.

In the text

Fig. A.1 Evolution of the rotation period (left panel) and of the surface velocity (right panel) for an isolated 1 M⊙ stellar model during the PMS phase and the very early MS phase for different values of the initial rotation. The green line shows respectively the critical period (left panel) and the critical velocity (right panel). We note that the bottom line in the left panel corresponds to the upper line in the right panel.

In the text

Fig. B.1 Variation of the planetary radius and of the Roche lobe as a function of the planetary mass. Left panel: variation of the planetary radius as a function of the planetary mass (green curve; see text). The red curves indicate the effective radius accounting for the fact the atmosphere of the planet can be inflated. From bottom to top, the curves corresponds to values of the XUV flux between 104.5 and 105.5 in steps of 0.2 dex. These XUV fluxes corresponds to XUV luminosities equal to 10−3 L⊙ and orbital distances between about 0.05 and 0.2 au. Right panel: variation of the Roche lobe as a function of the planetary mass (red curves) compared to the planetary radius (thick green curve). The red curves corresponds to distances between the planet and the star equal to 0.01, 0.02, 0.03, 0.04, 0.05, and 0.1,au, from bottom to top. The mass of the star is 1 M⊙.

In the text

Fig. B.2 Left panel: variation of the energy needed to bring a unit-mass planetary material up to the Roche lobe as a function of the mass of the planet and of the orbital distance. The green line corresponds to the case when the planet is in isolation, showing the energy needed to bring a unit mass to infinity. The red lines corresponds to cases where the planet is sufficiently near the star so that the energy would be lesser than the former case, which is needed for bringing a unit mass only up to the Roche lobe between the planet and the star. From top to bottom, red lines correspond to orbital distances equal to 0.1, 0.05, 0.04, 0.03, 0.02, and 0.01 au. Right panel: variation of the initial evaporation timescale of planets of different masses irradiated by different XUV fluxes. The green line corresponds to a flux equal to 104.5 erg s−1. The lowest red curve corresponds to a flux equal to 105.5 erg s−1. The lines in between correspond to intermediate values, from 104.6 to 105.4 in steps of 0.1 dex. The distance between the planet and the star has been assumed to be 0.1 au for the continuous lines and 0.01 au for the dashed lines.

In the text