© ESO 2019

Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 Introduction

The planetary systems discovered during the last two decades show a wide and unexpected diversity. Indeed, the detection of 51 Pegasi b (Mayor & Queloz 1995), a Jupiter-like planet orbiting around its star in less than five days, questioned the theories of planetary system formation which were based on observations of the solar system. Recently, the discovery of Proxima b (Anglada-Escudé et al. 2016) and the planetary system of TRAPPIST-1 (Gillon et al. 2017) paved the way for research of habitable planets around low-mass stars. Understanding how such systems form and evolve is one of the most challenging questions in astrophysics.

A large proportion of systems where one planet or more is orbiting closer to its host star than Mercury to the Sun have been observed. Tidal interactions play a key role in the orbital configuration of these very compact systems since they are likely to circularize orbits, align spins, and synchronize periods (Zahn 1977; Mathis & Remus 2013; Ogilvie 2014). These interactions consists in an exchange of angular momentum between the orbit and the spins of the celestial bodies. This exchange is the consequence of the dissipation of tidal flows. The kinetic energy of these tidal flows is converted into heat through tidal dissipation. Since the planet is synchronized within a timescale of a few thousand years, the stellar tide drives the secular orbital evolution (Guillot et al. 1996; Rasio et al. 1996; Leconte et al. 2010). In this work, we neglect the impact of the dissipation in the radiative zone. In stellar convection zones, there are two kinds of tides and both are dissipated by the turbulent friction applied by convective eddies. On the one hand, the equilibrium tide is the large-scale velocity field associated with tidal deformation, the so-called tidal bulge. This nonwave-like entity corresponds to the hydrostatic adjustment of the star to the gravitational perturbation (Zahn 1966; Remus et al. 2012). The friction applied by convective motions delays the response of the star to the perturbation (e.g., Zahn 1989; Ogilvie & Lesur 2012; Mathis et al. 2016). This results in a lag angle between the axes of the tidal bulge and the line of centers. This angle increases with dissipation magnitude. Hansen (2012) calibrated its value for several stellar masses by constraining the dissipation using observations of planetary systems. Since lower-mass stars have deeper convective envelopes, they dissipate more energy than higher-mass stars. On the other hand, in rotating bodies such as stars, at low tidal frequencies, the Coriolis acceleration acting on this equlibrium tide excites inertial modes (Ogilvie & Lin 2007). Their ensemble, the dynamical tide, constitutes the wavelike part of the tidal response. Its dissipation strongly depends on internal structure since it arises from the reflection of inertial modes on the radiative, stably stratified core (Ogilvie 2013; Mathis 2015). The dynamical tide may also vary over several orders of magnitude with rotation since inertial waves are restored by the Coriolis force. At low frequencies, dissipation of the dynamical tide is several orders of magnitude higher than the dissipation of the equilibrium tide (Ogilvie & Lin 2007).

Orbital evolution occurs simultaneously with variations of stellar structure and rotation (Gallet & Bouvier 2015; Amard et al. 2016; Bolmont et al. 2012; Bolmont & Mathis 2016). In the case of solar-like stars, the latter is slowed down over most of the main sequence (MS) by magnetic braking. This phenomenon occurs because of the wind carrying angular momentum away from the star (Schatzman 1962; Weber & Davis 1967). Models of secular evolution of stellar rotation generally consider that it undergoes three main phases. First, before the disk dissipates, stellar rotation remains constant. The physical processes that balance accretion and contraction have long been thought to have a magnetic origin. Matt & Pudritz (2005) and Matt et al. (2012) investigated the braking caused by accretion-powered stellar winds while Zanni & Ferreira (2013) studied the effect of magnetic ejections on stellar spin. Recently, Bouvier & Cébron (2015) explored the possibility that tidal and magnetic interaction with a close-in planet embedded in the disk could compete with accretion and contraction. After the disk dissipates, the star spins up due to its contraction during the pre-main sequence (PMS; see Amard et al. 2016, and references therein). Finally, once on the MS, the star spins down under magnetic braking. Observations show that rotational velocities of young stars range from one to one hundred times the solar velocity whereas evolved stars tend to converge to the solar rate on the Skumanich sequence (Skumanich 1972). Gallet & Bouvier (2015) investigated the mass-dependence of rotational evolution and showed that the braking torque and the core-envelope coupling timescale strongly depend on stellar mass. Matt et al. (2015) use the Rossby number to disentangle solar-like star populations in two groups: fast, saturated rotators and slow, unsaturated rotators. They showed that the spin-down timescale was decreasing with stellar mass in the former group and increasing in the latter, in agreement with the observations of Barnes (2010).

In star–planet systems, the interactions between the central body and its companion result in intricate phenomena that involve the structure of the star, rotation, and magnetism in addition to the orbital parameters of the planet. For example, the spin-down of a solar-like star over the MS increases its corotation radius and this may occur until the latter becomes larger than the orbital semi-major axis, causing a change of sign of the tidal torque. Moreover, the planet spiralling inward may spin up its host star, thus impacting its magnetism through dynamo processes (e.g., Brun et al. 2004, 2015). Recent numerical and theoretical works have allowed significant advances in our comprehension of these mechanisms. Dobbs-Dixon et al. (2004) studied the combined effects of tidal dissipation and magnetic braking to explain the observed distributions of orbital eccentricities. Barker & Ogilvie (2009) investigated the influence of these effects on spin alignment. Ferraz-Mello et al. (2015) used the creep tide theory from Ferraz-Mello (2013) and a semi-analytical wind model (see Bouvier 2013, and references therein) to compute the past evolution of observed star–planet systems. Damiani & Lanza (2015) assumed a constant tidal efficiency and a Skumanich-type wind-braking law, in which the torque is proportional to the stellar rotation cubed. These authors demonstrated that a pseudo-stable equilibrium state can exist for star–planet systems, in which corotation is not achieved and the ratio of the orbital mean motion divided by the stellar rotation rate is determined by the angular momentum loss rate due to magnetic braking. Zhang & Penev (2014) implemented a star–planet system secular evolution code based on the two-layer rotational model of MacGregor & Brenner (1991) and the constant quality factor framework from Goldreich (1963) and performed a statistical analysis on their numerical simulations to constrain the tidal theory. Bolmont & Mathis (2016) adapted the frequency-averaged results of Mathis (2015) to the time-lag framework (Mignard 1979) and studied its impact on secular evolution of star–planet systems when coupled with a one-layer rotational model. Despite the progress they have made toward the comprehension of planetary systems, all the above cited studies relied on a parametrized description of tidal dissipation or wind braking.

In this first study, we consider systems constituted by a solar-like star and a planet. We characterize the secular evolution of stellar rotation and orbital parameters under magnetic braking and tidal dissipation. Following Zhang & Penev (2014) and Bolmont & Mathis (2016), we aim at developing a tidal model that takes into account the dissipation of both equilibrium and dynamical tides in the convective envelope and their dependence on stellar structure and rotation. We use a two-layer model of stellar interior (MacGregor & Brenner 1991; MacGregor 1991) to study the evolution of rotation. In Sect. 2, we present the hypotheses of our study and detail the interactions that take place in the considered system. In Sect. 3, we explain how we computed the torque of Réville et al. (2015a) as a function of structural and dynamical parameters of the star and compare our results with the model of Matt et al. (2015). As discussed in Sect. 4, to quantify the tidal torque, we used the empirically calibrated results of Hansen (2012) for the equilibrium tide and the theoretical prediction of Mathis (2015) for the dynamical tide. In Sect. 5, we use ESPEM (French acronym for Evolution of Planetary Systems and Magnetism) to compare our model to those of the literature and explore the influence of the fundamental characteristics of systems on their fate, to determine the survival of the exoplanet and its influence on stellar rotation. Finally, in Sect. 6, we summarize our results and detail the perspectives opened by this work.

2 Star–planet interaction model

We consider systems formed by a planet orbiting a solar-like star. We used the two-layer star model of MacGregor & Brenner (1991) to study stellar rotation. In this framework, both the radiative core and the convective envelope rotate uniformly. They exchange angular momentum and only the envelope is directly braked by the wind. Thus, the core is also braked but this depends on the choice of the coupling timescale, since here we do not explicitly solve for the physical mechanisms (magnetic field, waves, turbulence) that are likely to cause this exchange (see for instance Rudiger & Kitchatinov 1997; Spada et al. 2010; Brun et al. 2011; Strugarek et al. 2011; Talon & Charbonnel 2005; Spiegel & Zahn 1992, and references therein). The basic idea of this model is that, in the absence of any perturbation, the core and the envelope evolve towards synchronization of their spins. The amount of angular momentum ΔL that the core andthe envelope should exchange to achieve this final state, positive if the core gives it to the envelope, negative in the opposite case, is (1)

where I_c and Ω_c (respectively I_r and Ω_r) are the moment of inertia and rotational velocity of the convective (respectively radiative) zone. The timescale of angular momentum redistribution, τ_int, is a free parameter of the model.

Another process occurring within the star causes angular momentum exchange between the core and the envelope. During the PMS, while the radiative zone is growing in size and mass, shells of matter rotating at the same velocity as the envelope are deposited on the core, which results in an angular momentum transfer from the envelope to the core. Thus, the internal torque Γ_int corresponding to the angular momentum exchange per time unit given by the radiative zone to the convective zone is given by (2)

where M_r and R_r are the mass and radius at the base of the envelope, respectively. The first term of the equation is due to simple parametrized angular momentum redistribution while the second is caused by stellar evolution (MacGregor 1991).

We used the STAREVOL evolutionary tracks (Amard et al. 2016) to get the one-dimensional (1D) internal structure of the star at each ESPEM time step. This model computes the stellar mass and radius, the mass and radius at the base of the envelope, and the moments of inertia of the radiative and convective zones and their variations along the evolution. We studied stars of mass ranging from 0.5 to 1.2 M_⊙.

The orbit is assumed to be circular and coplanar and we focus on variations of the semi-major axis. Moreover, we consider the rotation of the planet to be synchronized. Indeed, Guillot et al. (1996) showed that the synchronization of a hot Jupiter occurs within a typical timescale of 10⁶ yr, which is negligible compared to the lifetime of a star–planet system. Thus, the planet can be considered as synchronized from the beginning of the simulation. Consequently, in our model, only tidal dissipation within the star impacts orbital evolution. As a simplification for this first study, we assume that this process takes place in the envelope, which involves that the tidal torque is applied on the convective zone only, and not on the core. In our future works, the dynamical tide in the radiative zone, which is constituted of gravity waves dissipated by thermal diffusion and breaking, will have to be implemented in our code (Zahn 1975; Goldreich & Nicholson 1989; Terquem et al. 1998; Ogilvie & Lin 2007; Barker & Ogilvie 2010). In this first step, we focus on the convective zone, which allows consistent comparison with Zhang & Penev (2014) and Bolmont & Mathis (2016).

The core is interacting with the envelope through internal coupling as detailed in Eqs. (1) and (2). The orbit is exchanging angular momentum with the envelope. The latter acts as an interface between the planet and the core and loses angular momentum through magnetic braking. Consequently, the semi-major axis of the orbit a, the angular momentum of the convective zone L_c, and of the radiative zone L_r are determined by the following equations:

where G is the gravitational constant, m_p is the planetary mass, M_* is the stellar mass, and Γ_wind and Γ_tide are the wind braking and the tidal torques, respectively. The latter two correspond to the rates of change of angular momentum associated with wind braking and tidal dissipation. Their calculation is presented in the following sections. We used the ODEX solver (Hairer et al. 2000) implementing the Bulirsch–Stoer algorithm to solve these equations and compute the secular evolution of stellar rotation and planetary orbit. In the following, we detail how the braking and tidal torques were computed. The different entities and their interactions are shown in a simplified view in Fig. 1.

When planets come too close to their host star, they are destroyed either by tidal forces or after having plunged into the stellar atmosphere. To model this phenomenon in ESPEM, we followed Zhang & Penev (2014) who used the ratio of the density of the planet divided by that of the star ρ_p∕ρ_* to determine which scenario occurs. Their approach was based on the observational results of Metzger et al. (2012), who showed thatthree scenarios were likely to take place when a planet inspirals towards its host star. In the first case, ρ_p ∕ρ_* > 5 and the planet spirals until it plunges in the stellar atmosphere. In the second case, 1 < ρ_p∕ρ_* < 5 and the planet overflows its Roche lobe before reaching the stellar surface. This results in an unstable mass transfer from the planet to the star and the former is torn apart within a timescale of several hours. In the third case, ρ_p ∕ρ_* < 1, the planet also overflows its Roche lobe before coming into contact with its star but the mass transfer is stable. Consequently, the planet is disrupted over a much longer timescale, typically several thousand years. Even if these three cases lead to various timescales, planet destruction lasts only a short time compared to stellar lifetime. This is why we consider that the entire process occurs instantaneously. If ρ_p ∕ρ_* > 5, the planet isremoved from the simulation if the semi-major axis becomes smaller than the stellar radius. In the other two cases, the planet is destroyed if it starts orbiting the star below the Roche limit, which we calculated with the formula of Zhang & Penev (2014): (6)

where r_p is the planetary radius, M_* is the stellar mass, and m_p is the planetary mass. We adopted the assumption that the planets have the same mean density as Jupiter, that is, 1.33 g cm⁻³. This is a reasonable choice for planets more massive than 3 × 10⁻²M_Jup that are known to be gaseous (Baraffe et al. 2014). We only considered planets with a mass above this value, which corresponds to the upper mass limit of super Earths (~10 M_⊕). We made this simplifying choice because the planets that raise the most significant tides within their host star are the most massive ones.

3 Stellar rotation

In this section, we compute Γ_wind. Magnetic braking occurs because of the stellar wind carrying angular momentum away (Schatzman 1962). Parker (1958) showed that the particles of the wind were accelerated along their way through the corona. The Alfvén radius r_A is the distance at which their speed reaches the Alfvén speed. In a simplified 1D model, Weber & Davis (1967) showed that the resulting angular momentum loss was equal to the product of the stellar convective zone rotation Ω_c, the mass loss and the Alfvén radius squared. After integration on a sphere, the equality becomes a proportionality relation: (7)

In this product, the only easily measurable factor is the stellar rotation. This can be done either by spectroscopy (Reiners & Schmitt 2003) or by photometry (McQuillan et al. 2013a; García et al. 2014). Such measurements of the other factors are however more difficult. This is why we used the model of Réville et al. (2015a) to express them as a function of the structure and dynamics of the star. We note however that other options have been used in the literature (e.g., Brown 2014; Gallet & Bouvier 2015; Johnstone et al. 2015a,b; van Saders et al. 2016; Sadeghi Ardestani et al. 2017, and references therein).

Réville et al. showed that the Alfvén radius was only dependent on the magnetic flux through the open field lines, Φ_open, for any given topology: (8)

where R_* is the stellar radius, is the escape velocity, is the ratio of rotational velocity of the envelope over the brake-up rotation rate, and the constants K₃, K₄, and m were set by the authors using numerical simulations (Réville et al. 2015a). There are two kinds of magnetic field lines: the closed ones, also known as magnetic loops, and the open ones. The latter coincide with the wind streamlines in the interplanetary medium. Thus, the open flux measures the magnetic flux in the wind; it depends on topology, which is the repartition of magnetic energy in the different spherical harmonics. Higher-order topology magnetic fields decrease more steeply with distance from the star than lower-order ones: a multipolar field of degree l decreases as 1∕r^l+2 where r is the distance to the star. This results in smaller magnetic fluxes at equal distances. Therefore, the open flux decreases with the order of the magnetic multipole considered.

Equation (8) has been obtained with 2D numerical simulations. Réville et al. (2016) verified this result and revised the list of dimensionless constants (K₃, K₄, m) with 3D simulations constrained by realistic magnetic mappings obtained by spectropolarimetric measurements and concluded that K₃ = 0.55, K₄ = 0.06, and m = 0.3.

Formulating the torque of the wind as a function of the open magnetic flux allows us to write a topology-independent law. This idea was also used in the observational study of See et al. (2017), who estimated the open flux of 66 solar-like stars from Zeeman–Doppler magnetograms. After Réville et al. (2015a,b) and Garraffo et al. (2015), Finley & Matt (2017, 2018) used the open flux to study the wind of stars with complex magnetic topologies involving interactions between the dipolar, quadrupolar, and octupolar components. In the following subsections, we describe the steps of our method to estimate the wind-braking torque from the stellar structure and rotation.

3.1 Mass-loss rate and open magnetic flux

Computing the torque of the wind requires to first calculate the mass-loss rate and the open magnetic flux. Since these quantities are determined by the properties of the open magnetic field lines, a consistent model of the corona is needed to infer their global properties. To that end, we used the starAML routine developed by Réville et al. (2015b), who presented a method to determine the maximal size of the main coronal streamer for a given star with known mass, radius, effective temperature, surface magnetic field, base coronal density, and temperature (see Appendix A for more details). The result of this calculation corresponds to the radius of the surface beyond which all field lines are opened by the wind, the so-called source surface (Schatten et al. 1969; Altschuler & Newkirk 1969; Schrijver & DeRosa 2003). Figure 2 illustrates the coronal structure of a solar-type star with its magnetic structures.

The method consists in extrapolating the surface magnetic field to compute the source surface radius from the equilibrium between the ram and thermal pressures of the gas and the magnetic pressure. The mass-loss rate and open magnetic flux are then inferred from this calculation. First, the radial profiles of the temperature and density of the wind are computed assuming a polytropic equation of state of index γ = 1.05. These profiles allow us to compute the velocity field in the corona. The mass-loss rate is deduced from these computations as Ṁ =4 πr²v²ρ², where r, v, and ρ correspond to the distance to the star, velocity, and mass density at the outer boundary of our calculation grid, respectively.

Subsequently, the magnetic field is extrapolated from the surface to the rest of the corona. From these quantities, the thermal and ram pressures of the wind can be computed in the whole corona as well as the magnetic pressure. Close to the stellar surface, the magnetic pressure dominates the other two and is sufficient to hold the magnetic field lines closed. Far from the star, the pressure difference changes sign and the ram and thermal pressures of the wind open the magnetic field lines. The surface at which this inversion occurs is defined as the source surface. If we approximate this surface by a sphere, it is possible to define its radius, r_ss. We computed this radius from the radial profiles of the thermal, ram, and magnetic pressures and deduced the open magnetic flux, which is equal to the flux of the magnetic field through this surface (see Fig. 2). This method requires assumptions on the density n_c and temperatureT_c at the baseof the corona. These boundary conditions are detailed below.

3.2 Assumptions at the base of the corona

For the solar values, we followed Réville et al. (2016) and used n_⊙ = 10⁸ cm⁻³ and T_⊙ = 1.5 × 10⁶ K, which accurately reproduce the speed of the solar wind measured on Earth. For other stars, we followed Holzwarth & Jardine (2007) and Réville et al. (2016), who extended these hypotheses assuming that the temperature T_c and density n_c at the base of the corona depended mostly on the surface rotation rate:

A more recent model on mass-loss rates of cool MS stars (Cranmer & Saar 2011) predicts a slightly steeper dependence of n_c on Ω. However, the scatter in their results still allows for a dependence like our Eq. (9). Please note that alternative laws for the variations of coronal temperature with respect to global stellar parameters have recently been proposed (e.g., Johnstone & Güdel 2015; Wood et al. 2018). These formulations are discussed in Sect. 6.

At fast rotations, we introduced saturation by considering that n_c and T_c remain constant.

3.3 Saturation

To determine the saturation rate, we followed the approach of Matt et al. (2015), who used the stellar Rossby number to define it: (11)

where τ_conv is the convective turnover timescale, which we computed with the formula of Cranmer & Saar (2011) that is mostly valid on the MS, using the value of the effective temperature at the ZAMS: (12)

We kept the effective temperature at the ZAMS in this expression because it allowed us to simplify the calculations with starAML. This is a reasonable assumption since the variations of T_eff are negligible during the MS phase (<5%). This point is further discussed at the end of Sect. 6.

The saturation value of the Rossby number was set to R_o,sat = 0.1 R_o,⊙. The convective turnover timescale is indirectly linked to stellar mass. Indeed, among solar-type stars, lower-mass stars have thicker envelopes and slower flows. Therefore, their convective cells travel on longer timescales, which is why the turnover time is a decreasing function of stellar mass. Consequently, the Rossby number contains information about rotation and mass at the same time. It is also a good indicator of magnetic activity, as shown by Noyes et al. (1984) who emphasized the correlation between chromospheric activity and Rossby number, and Pizzolato et al. (2003) and Wright et al. (2011) who found evidence of saturation of X-ray emissions at R_o = 0.1 R_o,⊙.

3.4 Global surface magnetic field

Solar-type magnetic fields are believed to be generated by a dynamo in the envelope (Brun et al. 2004, 2015). Therefore, their strength depends on the depth of the convective zone, which decreases with stellar mass, and the rotation of the star (Noyes et al. 1984). This leads us to consider the following law for the magnetic energy at the surface:

The solar value was set to B_⊙ = 3 G, which is in good agreement with observations (Vidotto et al. 2014).

We assumed a dipolar magnetic topology to compute the braking torque. The Alfvén radius from our calculations was 5.25 R_* for a solar-mass star rotating at the solar rate. However, this value was not sufficient to brake the 1 M_⊙ star efficiently and reproduce the convergence of rotational velocities at solar age. This could be related to the “Open Flux Problem”, which reveals that magnetohydrodynamics wind models consistent with surface observations systematically underestimate the interplanetary magnetic flux (see Linker et al. 2017). To fix this lack of braking, we artificially multiplied the Alfvén radius by a factor 3.6. This correction allows us to reconcile our theory with observations of open clusters.

3.5 Braking torque

Figure 3 shows mass-loss rate and braking torque as functions of rotational velocity for several stellar masses. As can be seen in the upper panel, the mass-loss rate increases with rotation until saturation. This is similar to the behavior of density and temperature at the base of the corona. At high rotation, a dependence on stellar mass appears because stars have different saturation velocities. Despite the fact that our model predicts a mass-loss rate similar to that of Cranmer & Saar (2011) for a solar-mass star at solar rotation rate, our results suggest a weaker dependence on stellar mass and rotation rate. This discrepancy stems from the additional constraint we put on the spin-down timescale. Barnes (2010) and Matt et al. (2015) suggested that this quantity decreases with stellar mass for saturated rotators and increases for unsaturated ones (see Sect. 3.6). Since we aimed at reproducing this behavior, we had to impose a different mass-loss rate from that of Cranmer & Saar (2011).

The lower panel shows the braking torque resulting from our calculations. As for mass-loss rate, the unsaturated and saturated regimes are clearly visible. The dependence on rotation can be approximated by two power laws:

where p = 2.11. These exponents are close to those found by Kawaler (1988). We point out that the dependence on mass is negligible in the unsaturated regime and becomes important in the saturated regime. At 10 times the solar rotation rate, the torque is one order of magnitude higher for a 1 M_⊙ star than for its 0.5 M_⊙ counterparts.

3.6 Spin-down time

From these properties of the stellar coronae, we can now discuss the characteristic timescale of magnetic braking. Barnes (2010) showed that the dependence of magnetic breaking on mass was not the same for fast and slow rotators. Among the former, more massive stars tend to spin down more quickly, whereas this trend reverses for slow rotators. Matt et al. (2015) defined two different spin-down times for the saturated and unsaturated regimes and reproduced this result. To investigate the mass dependence predicted by ESPEM, we calculated the equivalent expression of these timescales as functions of the angular momentum of the star and the braking torque:

where I_* is the stellar moment of inertia. Equations (15) and (16) indicate that, with these definitions, τ_sat and τ_unsat do not depend on rotational velocity. Indeed, in both saturated and unsaturated regimes, the braking torque increasing with rotation rate compensates the explicit dependence on Ω_c. Thus, they allow a comparison of spin-down timescale for different stellar masses, different ages, and the same state of evolution.

Figure 4 shows the spin-down timescale as a function of stellar mass for both saturated and unsaturated regimes. The slopes of both curves are determined by the competition between two quantities increasing with stellar mass: the moment of inertia on the one hand, and the braking torque on the other hand. Since the latter does not vary significantly with mass for slow rotators, the corresponding timescale increases proportionally to the moment of inertia. On the contrary, for fast rotators, this competition results in a spin-down timescale decreasing with mass (). The figure also shows a very good agreement of our results with the model of Matt et al. (2015).

3.7 Secular evolution

We applied this braking torque to stellar evolutionary tracks for masses 0.5, 0.8, and 1 M_⊙ computed with the STAREVOL code (Amard et al. 2016) to assess the variations of the rotational velocity of solar-type stars on secular timescales. For each star, we considered three different initial conditions: a slow, a median, and a fast rotator, corresponding respectively to the 25th percentile, median, and 90th percentile of the observed distribution of rotational velocities of young stars. We used the parameters given by Gallet & Bouvier (2015) to set the free parameters of our model: initial rotation period, coupling constant, and disk lifetime (see Table 1).

The rotation of the envelope was kept constant from the beginning of the simulation until a time corresponding to the dissipation of the disk. During this phase, the radiative zone is free to spin up under contraction. However, it is linked to the external layer of the star through internal coupling. After the disk dissipation, rotational velocities of both envelope and core evolve under the action of stellar evolution, coupling, and magnetic braking, as described in Sect. 2.

Figure 5 shows the evolution of rotation of stars of masses 0.5, 0.8, and 1 M_⊙ from the PMS to the end of the MS. For each mass, the initial spread in rotational velocities is reduced over the MS and all stars converge on a unique sequence. Our model predicts that rotation of evolved stars is proportional to t^−0.473 in accordance with the Skumanich law (Skumanich 1972). As is visible in the top panel, at the age of the Sun, solar-mass stars rotate at the solar rate. Less-massive stars are slower at this age, which had been found by Gallet & Bouvier (2015). This result is also in agreement with those of Matt et al. (2015), who noted that the braking timescale in the unsaturated regime is shorter for lower-mass stars.

4 Tidal dissipation

In this section, we detail how the tidal torque Γ_tide is computed. We work within the tidal quality factor formalism (Goldreich 1963; Kaula 1964; MacDonald 1964). In this framework, the response of the star to tidal perturbation is measured by the modified tidal quality factor, Q′, which is the ratio of the total energy stored in the tidal velocity field divided by the energy dissipated over one revolution of the planet. This convenient formalism has been used to account for stellar structure and rotation variations over the secular evolution of the system (Mathis 2015; Bolmont & Mathis 2016; Gallet et al. 2017b). Mathis & Le Poncin-Lafitte (2009) showed that when the planet is far enough (a > 5 R_*) and weakly deformed, it could be considered as a point-mass disturber for calculating tidal dissipation within the star. In this case, the deformation of the star due to gravitational perturbation is quadrupolar and only the corresponding frequency, ω_tide = 2 (Ω_c − n), is excited in the circular coplanar case. As the planet gets closer to its star, which corresponds to a semi-major axis inferior to 0.018 AU in the case of a solar-radius star, this hypothesis is no longer valid. However, since the planets orbiting their star below this limit are destroyed after a few thousand years in our calculations, we do not take this complication into account in this work. Nevertheless, differences arising from the fact that the planet is an extended body should be studied in a future work. Tidal frequency depends on the rotation rate of the envelope because we only consider the dissipation in the convective envelope of the star. The torque depends on the dissipation rate and the structural parameters of the star and the planet (Murray & Dermott 1999):

(19)

where m_p is the mass of the planet, assumed to be constant over the evolution of the system. We interpolated the STAREVOL evolutionary tracks (Amard et al. 2016) to obtain the values of stellar structure at each ESPEM time step. Since tidal dissipation originates from two different physical mechanisms, the equilibrium tide and the dynamical tide, each of them have to be studied separately. The total tidal dissipation is proportional to the inverse of the equivalent quality factor. It is the sum of the two contributions: (20)

where and are the modified quality factor related to the equilibrium and the dynamical tide, respectively.

4.1 Equilibrium tide

To compute the equivalent quality factor related to the equilibrium tide, we used the observational results of Hansen (2012), who calibrated the value of the dissipation for stellar masses ranging from 0.3 to 1.5 M_⊙. He worked in the constant time lag framework (Mignard 1979; Hut 1981; Eggleton et al. 1998), which is different from ours. Therefore, we had to reconcile both formalisms. The equilibrium tide corresponds to the hydrostatic adjustment of the star to the perturbation of the planet. In this work, we identify it to its quadrupolar moment. In the adiabatic case, the axis of the tidal bulge would be aligned with the one joining the centers of the star and of the planet. However, dissipation induces an angle of 2δ between these axes. In the quality factor framework, dissipation is proportional to this angle and Q′ is related to the lag angle through the relation 2δ = 3∕(2Q′). The time lag Δτ = 2δ∕|ω_tide| corresponds to the time delay between the positions of the bulge and the axis joining the centers. In the constant time lag framework, dissipation is quantified by the constant σ_*, which measures the ratio of energy loss due to tidal friction divided by the magnitude of the quadrupolar moment of the deformation. Hansen (2012) found an empirical law between the normalized dissipation , where , and the stellar parameters: (21)

Following Bolmont & Mathis (2016), we inverted the relation to obtain the corresponding modified tidal quality factor: (22)

We point out that the quality factor associated with the equilibrium tide is inversely proportional to tidal frequency. Thus, we expect different predictions from those of models assuming a constant quality factor (Zhang & Penev 2014). This dependence might cause computational problems close to corotation. In practice, we replaced |ω_tide| in Eq. (22) by |ω_tide| + ε_tide where ε_tide = 10⁻¹⁰ rad s⁻¹. This technique was introduced by Bolmont & Mathis (2016), who showed that this does not affect the predictions for the final states of the system. Moreover, since , that is, Γ_tide ∝sign(ω_tide)ω_tide with Eq. (19), there is no actual singularity of the torque here.

4.2 Dynamical tide

At low frequency, the time-varying gravitational potential excites inertial waves in the convective envelope of the rotating star. This phenomenon occurs for ω_tide ∈ [−2 Ω_c,  2 Ω_c] and constitutes the wavelike part of the tidal response. Its dissipation may be several orders of magnitude larger than that of the equilibrium tide. Ogilvie & Lin (2007) first calculated the dissipation of the dynamical tide within the interiors of solar-type stars. Their result is highly dependent on tidal frequency because of resonances induced by the wavelike nature of the dynamical tide. Moreover, the properties of the resonances strongly depend on the assumed values of the eddy coefficient applied to tidal waves to account for the turbulent friction applied by convection. These latter values are still poorly known, as pointed by Ogilvie & Lesur (2012), Ogilvie & Lin (2007), Guenel et al. (2016), and Mathis et al. (2016), for example. This would imply heavy calculations to study secular orbital evolution with complex hydrodynamical simulation for each step of the structural and rotational evolution of the star (Witte & Savonije 2002; Auclair-Desrotour et al. 2014). To address this situation, Ogilvie (2013) computed the frequency-averaged value of the tidal waves dissipation using an impulsional method. Over the range [−2 Ω_c,  2 Ω_c], this leads to a result that depends only on the structure and rotation of the star. Mathis (2015) then applied this method to the envelopes of solar-like stars to compute the equivalent modified quality factor associated with the dynamical tide: (23)

In this formula, α = R_r∕R_* and β =M_r∕M_* are the mass and aspect ratios of the core, respectively, γ = α³(1 − β)∕(β(1 − α³)) is the ratio of the density of the envelope to that of the core, and . We computed these ratios from the stellar evolution tracks of the STAREVOL code (Amard et al. 2016), which is well adapted to our bilayer stellar rotation model. Mathis (2015) showed that, for a given rotation rate, the dissipation could vary over several orders of magnitude with aspect and mass ratios. Moreover, a solar-like star rotation rate may vary during its life on the MS in the range [Ω_⊙,  100 Ω_⊙] (Gallet & Bouvier 2015). Since the dissipation is proportional to the rotation rate squared, a given star may dissipate significantly more at its arrival on the MS than at solar age, as discussed by Gallet et al. (2017a). In this paper, we focused on stars with solar metallicity. The impact of this parameter should be carefully studied in a further work, as done by Bolmont et al. (2017), who showed that tidal dissipation in the convective zone of solar-like stars was higher for metal poor stars on the PMS and that this trend was inversed on the MS.

Figure 6 shows the tidal quality factor Q′ from Eq. (20) for a star of mass 1 M_⊙ at solar age as a function of rotation rate and orbital mean motion. It is visible that tidal dissipation is the sum of two contributions. The white line with dark blue contour delimits the domain of application of inertial waves. Below, the tidal frequency is in the range [−2 Ω_c,  2 Ω_c] because n <2 Ω_c. In this domain, the dissipation of the dynamical tide dominates that of the equilibrium tide. It is visible on the figure that the dynamical tide contributes to significantly enhance the tidal dissipation. Indeed, the quality factor in the region where they are raised is orders of magnitude lower than in the upper part of the figure. The dependence of the dynamical tide on rotation also clearly appears. Indeed, in the lower region, the contour lines are vertical since the dependence on tidal frequency of the dissipation of the dynamical tide was averaged in this work. On the contrary, the contour lines in the upper region are actually determined by ω_tide (see Eq. (22)).

This figure also reveals three possible regimes of tidal interaction. Below the light blue line, the orbital motion is slower than the stellar rotation rate. This region corresponds to planets beyond the corotation radius. The planets are pushed away as a consequenceof the dissipation of the equilibrium and dynamical tides. Between the blue lines, the dynamical tide is excited in the stellar envelope but the planets are situated below the corotation radius. Consequently, they spiral inward under the effect of both the equilibrium and dynamical tides. Above the dark blue line, inertial waves are no longer raised. In this region, planets spiral inward under the sole influence of the dissipation of the equilibrium tide. The secular evolution of a star–planet system generally involves the succession of different phases during which the system is in one of these three states.

5 Results

In this section, we analyze the coupled influence of magnetic braking and tidal dissipation on the secular evolution of star–planet systems. Webegin by showing that the rotation of a star can be significantly impacted by the presence of a planet, especially when the latter falls into the former. The fate of a star–planet system can follow very different scenarios, depending on parameters such as the mass of the star and initial rotation rate, the planetary mass, and the initial semi-major axis of the orbit. To investigate these dependences, we browse the parameter space to study which combinations lead to the planet falling into the star. We then define a criterion allowing to say, for a given system, whether the planet will survive or be destroyed.

5.1 Impact of a planet on the rotation rate of its host star

We first investigated the consequences of tidal dissipation in the convective envelope on stellar rotation. To that end, we computed the secular evolution of a star–planet system composed of a solar-mass star and a Jupiter-mass planet with an initial semi-major axis equal to 0.025 AU. In order to compare the results with those of Fig. 5, we calculated this evolution for three different initial rotation periods: 1.4, 5, and 8 days.

As visible in Fig. 7, the stellar rotation can be strongly impacted by the presence of a planet. The most noticeable differences with Fig. 5 are the two peaks in the curves of the slow rotator (dark blue) and the median rotator (medium blue). They are both caused by the destruction of the planet. Indeed, as is further analyzed later in this section, when the planet orbits the star below its corotation radius, the orbit gives angular momentum to the stellar spin, causing the central body to spin up and the planet to migrate closer to the star. This process is maintained until either the semi-major axis crosses the corotation radius, in which case the sense of the tidal exchange of angular momentum is reversed, or the planet is disrupted after having approached too close to the star. In the latter case, the angular momentum transfer from the orbit to the stellar spin during the inward spiral motion of the planet is important enough to increase the rotation rate of the envelope by one order of magnitude. After the planet has been destroyed, the wind brakes the envelope rapidly so that a peak is formed in the rotation rate evolution curve.

The bottom panel shows the difference δP between the stellar rotation period in the case where the star hosts a planet and the case where it does not: (24)

This difference quantifies the influence of the planet on stellar rotation. For the medium- and dark-blue curves, it is negative because the star has been spun up. We note that the rotation remains fast after the destruction of the planet. For the dark-blue curve, δP is on the order of 12 days for the last 5 Gyr of the simulation, which means that the star spins twice faster than what gyrochronology predicts. In the case of the light-blue curve, the star has been spun down because its planet was initially located beyond its corotation radius. This is why δP is positive. Since the effect in this case is significantly smaller than in the other two simulations, the scale of the plot does not allow for a detailed analysis of this curve. We refer the reader to Appendix B.2 for a detailed analysis of positive-δP curves.

The fast rotator (light blue curve) seems unaffected by the presence of the planet. Indeed, contrary to the two other simulations, the planet in this case starts orbiting the star beyond its corotation radius. This leads to a fast outward migration which has weak effects on the stellar spin. This phenomenon, further discussed in Sect. 5.2.3, illustrates that the evolution of a star–planet system strongly depends on the parameters characterizing it. In the following, we study how the secular evolutions of the stellar rotation and the orbit differ when varying the initial conditions and parameters.

5.2 Planet survival time for a solar-mass star

We now investigate the influence of the characteristics of the system on its evolution. We identified four main parameters for this study: the mass of the star, its initial rotation, the mass of the planet, and the initial semi-major axis of the orbit. In this subsection, we only consider 1 M_⊙ stars and focus on how the fate of the system differs when the three other parameters vary. We first analyze the influence of each characteristic individually by comparing simulations done with different values of this parameter. Then, we study how their crossed variations impact the fate of the system by computing, for a large set of parameters, the time at which the planet is destroyed.

5.2.1 Initial semi-major axis

We started by looking at the impact of the initial semi-major axis on the fate of a system. To that end, we computed the secular evolution of a star–planet couple constituted by a 1 M_⊙ star with an initial rotation period P_*,ini = 5 days and a 1 M_Jup mass planet for three different values of the semi-major axis: 0.025, 0.035, and 0.045 AU. In the following, we compare the results obtained in each case and discuss the differences observed from one scenario to another.

It appears in Fig. 8 that the farthest planet (light green curves) is too far from the star and not heavy enough to undergo a significant migration. The closest planet (dark green curves) starts orbiting its star at a distance where only the equilibrium tide is raised. The orbital evolution in this case arises from the variations of the limit of excitation of inertial waves, which is represented by the dotted dark green line in the middle panel of the figure (confounded with the other two dotted curves during the first billion years). As the star spins up during its contraction, this limit decreases until it equals the semi-major axis. Then, the system achieves a resonance state in which n = 2 Ω_c. This resonance, characterized on the figure by the solid dark green curve sticking to the dotted dark green curve, is well visible in the bottom panel. Indeed, during this phase, the curve representing the quality factor is noisy, which is the consequenceof the dynamical tide being alternatively excited and switched off in the convective envelope. This state is maintained until the end of the contraction and as a consequence the planet is moved significantly closer to the central body within a short timescale (a few million years) compared to usual tidal evolution timescales. On the MS, the planet undergoes orbital evolution due to dissipation of the equilibrium tide and eventually spirals inward to its destruction. The planet in the middle (medium green curves) experiences a similar scenario. However, unlike the previous system, it does not undergo a resonance at the limit of application of inertial waves. This is because the dynamical tide is not as strong as in the former case. As in the case with a_ini = 0.025 AU, the planet with a_ini = 0.035 AU moves progressively closer to its star during the MS until its demise. Figure 8 shows that the initial semi-major axis directly influences the speed of tidal evolution: when it decreases, the system evolves faster.

5.2.2 Planetary mass

We now show how the mass of the planet influences the system. To achieve this, we performed a similar calculation to that of the previous paragraph. We computed the evolution of a star–planet couple with a solar-mass star initially rotating with a period P_*,ini = 5 days, orbited by a planet at an initial distance a_ini = 0.035 AU for three different planetary masses: 0.1, 1, and 5 M_Jup. The systemwith the Jupiter-mass planet is actually the same as the one in the previous paragraph for which a_ini = 0.035 AU. Therefore, Sects. 5.2.1 and 5.2.2 correspond to variations around the same point in different directions of the parameter space.

Figure 9 shows the evolution of three star–planet systems for a fixed semi-major axis and varying planetary mass. Similar to the farthest planet in Fig. 8, the lightest planet (black lines) is not influenced enough by tidal dissipation to undergo significant orbital evolution. The heaviest planet (orange lines) is disrupted during the early phases of the life of the system. As the semi-major axis of its orbit crosses the limit of excitation of inertial waves, the tidal torque is strong enough to raise a resonance between the orbit and the stellar spin. Unlike the closest planet of Fig. 8, the planet here is massive enough to significantly spin up its host star, which results in a decrease of the dynamical tide limit, represented by the orange dotted line in the middle panel. The resonance then imposes the planet to follow this evolution and migrate closer to the star. Thus, the semi-major axis of the orbit collapses in a few million years and the planet is destroyed. The medium-mass planet here corresponds to the same system as the middle planet of Fig. 8 (the one with a_ini = 0.025 AU). Its evolution is an intermediate case of the two previous cases. Figure 9 shows that tidal evolution is faster for systems with more massive planets. Systems with massive planets behave like systems with close planets since they evolve faster than others. However, a massive planet may spin up its star, which is not possible for a lighter and closer planet.

5.2.3 Initial stellar rotation rate

We finally reproduced the calculations of Sects. 5.2.1 and 5.2.2to study the influence of the initial stellar rotation rate. We set the planetary mass to 1 M_Jup and the initial semi-major axis to a_ini = 0.035 AU and computed the evolution of a system formed by this planet and a solar-mass star with an initial rotation period P_*,ini = 1.4, 5, and 8days, which correspond to the fast, median, and slow rotators defined by Gallet & Bouvier (2015). The results of these simulations are shown on Fig. 10.

In the system with the fastest rotating star, the planet starts its orbit beyond the corotation radius and is rapidly pushed away from the star as a consequence of the dissipation of the dynamical tide. When the star spins down during the MS phase, the corotation radius increases and eventually exceeds the semi-major axis of the orbit. Consequently, the tidal torque changes signs but, at this distance, tidal dissipation is not strong enough to cause significant orbital evolution and the planet eventually survives. The system with the median rotator, shown in regular blue, corresponds to the same system as the middle planet from Fig. 8 and the medium-mass planet from Fig. 9. Its evolution is discussed in the previous paragraphs. It is expected that higher initial stellar rotation rates lead to longer planet lifetimes if we consider the planet lifetime as the duration before it spirals towards its host star. However, the planet orbiting the slow rotator, shown in dark blue, lives longer than the one orbiting the median rotator. This can be explained by considering the contraction phase of thestar. For the slow rotator, the limit of excitation of inertial waves is farther away from the star than for the median rotator. Consequently, inertial waves, which result in an enhanced tidal dissipation, are raised over a longer duration in the latter case. Moreover, inertial waves in the median rotator are excited at a higher rotation rate than in the slow rotator. This leads to lower quality factors in the former case. To summarize, the dynamical tide is raised over a longer duration and induces higher tidal torques in the median rotator, which is why its planet is destroyed before the one orbiting the slow rotator. The results of Fig. 10 emphasize the impact of the contraction phase on the rest of the orbital evolution. Unlike the two precedent parameters, m_p and a_ini, the initial stellar rotation rate does not have a monotonic influence on the planet lifetime.

5.2.4 Browsing the parameter space

Determining how the characteristics of the system influence the planet lifetime is crucial to achieving a better understanding of orbital evolution. This is why we repeated the previous calculations on a wider and more detailed sample. We set the stellar mass to 1 M_⊙ and considered three initial rotation periods: 1.4, 5, and 8 days. For each period, we computed the secular evolution of systems with varying planetary mass and initial semi-major axis. For the former, the bin was equally spaced in logarithm with ten points between 10 M_⊕ and 10 M_Jup. For the latter, the bin was linearly spaced between 0.02 AU and 0.05 AU. We used the parameters of Table 2 from Gallet & Bouvier (2015) for the disk lifetime and internal coupling constant. For each simulation, we computed the planet lifetime as either the duration before the star terminates its MS phase in the case where the planet survives, or the duration before the planet is disrupted in the alternative case.

Figure 11 shows the planet lifetime as a function of planetary mass and initial semi-major axis for three different initial stellar rotation rates. Its value spans from a few million years (planet decayed during the PMS) to 10 billion years (terminal age on the MS). As can be seen in the figure, in this two-bodied simplified approach, for this mass range, planets which formed above 0.045 AU are very likely to survive with their host star. These results confirm the monotonic dependence of the planet lifetime on planetary mass and initial semi-major axis: apart from an outlier for massive close-in planets in the left panel, colors always become darker in the bottom-right corner, that is, closer and heavier planets are destroyed earlier. The blue dots represent the situation on these maps of the cases treated in detail in Figs. 8–10.

The thick black line represents the demarcation between the region where planets survive and the one where they eventually spiral inward to their demise; it behaves differently for the fast rotator compared to the two other cases. Indeed, in the right panel, the line is horizontal, which means that the survival of the planet does not depend on its mass, whereas in the others, it is oblique, which indicates that the survival of the planet results from a trade-off between its mass and distance from the host star.

This different behavior can be explained by considering the corotation radius associated with the initial stellar rotation period, as was analyzed in Sect. 5.2.3. For the slow and the median rotators, its value is higher than 0.05 AU, meaning that all planets shown on the left and middle panels of Fig. 11 start their orbit below the corotation radius of their host star. On the contrary, the initial corotation radius of the fast rotator is about 0.025 AU. As shown in the right panel of Fig. 11, all planets that start orbiting the fast rotator below this limit are rapidly destroyed, leading to a disentangling of the population into two groups: planets which formed above the corotation radius and rapidly pushed away from their stars, which allowed them to eventually survive, and those which formed below and were disrupted within the first million years.

This implies a depopulation of the region close to fast rotating stars, as observed by McQuillan et al. (2013b). Teitler & Königl (2014), who used a secular evolution code to explain this dearth, concluded that it was due to tidal engulfment. Our results suggest that this depopulation is actually due to two mechanisms: tidal disruption for the closest planets and fast outward migration inthe early stages of the evolution for the others, as shown for instance in Fig. 10. This second alternative solution was pointed out by Lanza & Shkolnik (2014). They found that neither the tidal quality factor framework nor the constant viscous time model (Eggleton et al. 1998) could reproduce the distribution of planetary orbital periods and favored a scenario in which close-in planets form at a distance of about 1 AU from their star and migrate closer to their star because of chaotic dynamic evolution caused by the other planets of the system.

In the cases of the slow and median rotators, the thick line involves more complex, intricate dependences on the characteristics of the system. The slope of the line in the (m_p,  a_ini) plane results from the competition between the planetary mass and semi-major axis in the expression of the tidal torque. This competition is analyzed below.

5.3 Survival-rate criterion

Here we determine a criterion allowing us to predict whether the planet will survive knowing the characteristics of the system (M_*, P_*,ini, m_p, and a_ini); we call this the survival-rate criterion. To this end, we repeated the simulations of Fig. 11 for stars of mass 0.5, 0.6, …, 1.1 M_⊙. For each star, we defined a slow, median, and fast rotator, as Gallet & Bouvier (2015). We then calculated the planet lifetime as a function of m_p and a_ini for each rotator. Finally, we defined two regions, the survival region, where planets survive with their host star, and the disruption region, where planets are eventually disrupted.

Figure 12 shows these regions for various stellar masses and initial stellar rotations. On each panel, the orange region indicates the planetary masses and initial semi-major axes for which the planet does not survive to tidal interaction; it is always located in the bottom right corner, which illustrates that the more a planet is massive andclose to its host star, the more it is likely to eventually decay. The figure also emphasizes the radically different behavior of fast rotators. Indeed, in the right column, the line delimiting the blue and orange regions is horizontal, suggesting as in Fig. 11 that the survival of the planet does not depend on its mass. For the median and slow rotators, behaviors are similar and the border between the two regions is a straight line, the slope of which is the same for all stellar masses and initial rotation rates. The decayed-planets region is wider for higher-mass stars, which is in agreement with the result of Bolmont & Mathis (2016) who found that tidal evolution has a stronger impact on systems with higher-mass stars despite the fact that the bulk dissipation is less important.

In each panel of Fig. 12, the thick black line separates the image into two half spaces. Obtaining the equation of this border as a function of the characteristics (M_*, P_*,ini, m_p, and a_ini) is sufficient to define the survival-rate criterion. For initially fast rotating stars, the survival of the planet is determined by its initial semi-major axis. If the latter is greater than the corotation radius, the planet survives; in the opposite case, it is rapidly decayed. For other initial stellar rotations, the other parameters play a more significant role in the fate of the system. Figure 10illustrates that the initial stellar rotation does have a clear influence on the survival of the planet. Therefore, we sought a relation involving only M_*, m_p, and a_ini to determine the equation of the border. For each stellar mass, we collected the coordinates of the points on the border in the slow and median rotator cases. Since increasing (resp. decreasing) the stellar or planetary mass (resp. the semi-major axis) reduces the chances of planet survival, we expected that all these points arrange according to the following equation: (25)

where α, β, and K are constants that we fitted with our simulated data.

Figure 13 shows the result of our fit of the border between the survival and disruption regions. The stellar mass is expressed in solar masses, the planetary mass in jovian masses and the semi-major axis in astronomical units. The points are scattered because the relation between the parameters on the border is not an ideal power law. Despite this dispersion, a trend clearly appears in the figure. We used a multidimensional least-squares method based on the Levenberg–Marquardt algorithm to fit the results of our simulations and found that the points follow the law with K = 43.5, α =1.12, and β =0.142. This equation allows us to define the survival-rate criterion:

6 Conclusion

We presented ESPEM, a code implementing a model of secular evolution of star–planet systems under magnetic braking and tidal interaction. Our wind model based on the work of Réville et al. (2015a) allows a fine analysis of magnetic phenomena occurring in the stellar corona. Its results are in agreement with recent observations (Barnes 2010; Gallet & Bouvier 2015) and it has the advantage of explaining the mass dependence of the braking torque (see Matt et al. 2015). We emphasized the difference between the ab initio tidal prescription of Mathis (2015) and the constant quality factor model used by Zhang & Penev (2014). The former may predict that a planet will survive whereas the latter predicts its decay. Therefore, understanding the physical processes at stake in tidal interactions is of prime importance to correctly interpret the various architectures of planetary systems discovered during the last decade. We showed some examples of complex orbital evolution scenarios that could occur as a consequence of tidal dissipation within the star. We showed that the spin-down caused by a planet orbiting a fast rotator could be detected only for planets more massive than Jupiter. We showed that, for hot Jupiters orbiting solar-like stars at a distance superior to 0.05 AU (~ 10.7 R_⊙), tidal interaction had no significant impact on the system. Finally, we defined a survival-rate criterion, allowing the prediction of the fate of a star–planet system (survival or decay of the planet) knowing its main parameters.

It is now of prime importance to consider the impact of tides within the stellar radiative zone, which are constituted by gravity waves and dissipated by turbulent friction or breaking (see e.g., Guillot et al. 2014). Indeed, this phenomenon is likely to enhance the tidal torque or compete with the dissipation in the envelope and modify the rotation rate of the stellar radiative core. It is therefore necessary to build a consistent model of tidal evolution. The present work can be considered as a step towards this aim. Implementing a physical description of internal angular momentum transport in the star is also among the perspectives of this paper (e.g., Zahn 1992; Mathis 2013, and references therein). Even if the simplified two-layer model implemented in ESPEM is already able to provide good insight into stellar rotation evolution, a more detailed view of the interior as developed by Amard et al. (2016) and Gallet et al. (2017b) would deeply improve our comprehension.

After the physical phenomena in the star are better treated in the code, considering tidal dissipation within the planet will be the next step. Indeed, whereas the spin angular momentum of the planet is by several orders of magnitude smaller than that of the star, tides raised in the former may lead to comparable tidal torques (see Mathis & Remus 2013; Ogilvie 2014, and references therein). Such theoretical studies could be compared to the recent results of O’Connor & Hansen (2018), who studied hot-Jupiter populations to constrain the dissipation in such planets. Taking this dissipation into account is required to study planet spin synchronization and alignment and orbital circularization and migration in the early stages of the life of a system (see also Heller 2018). The frequency dependence of the dynamical tide (Auclair-Desrotour et al. 2014), as well as the impact of eccentricity, inclination, and dynamical instabilities arising in multi-planet systems, have to be taken into account to compare theoretical predictions with the observed distributions (e.g., Bolmont et al. 2015; Damiani & Mathis 2018).

We will also consider the consequences of Roche-lobe overflow in a future work. Indeed, this process may cause mass transfer from the planet to the star which would result in an outward migration to conserve angular momentum (Trilling et al. 1998). In their theoretical study, Gu et al. (2003) showed that this scenario occurs in eccentric systems. More recently, Jackson et al. (2017) developed a consistent model of Roche-lobe overflow and showed that some observed systems were experiencing significant mass transfer. Such developments will be useful to improve our model of star–planet interaction.

In this work, we only considered stars with dipolar magnetic fields. However, the strength and topology of the global surface magnetic field of a star evolve with its age and structure (Gregory et al. 2012; Vidotto et al. 2014; Emeriau-Viard & Brun 2017). This may have an impact on the wind and stellar spin-down during the MS (Brun & Browning 2017; Metcalfe et al. 2016). Convection is able to influence magnetism via dynamo processes (Brun et al. 2004, 2015; Brown et al. 2011; Strugarek et al. 2017a; Augustson et al. 2015; Käpylä et al. 2014) and differential rotation (Brun et al. 2017). Determining the impact of these mechanisms is necessary to ensure realistic rotational evolutions of stars.

This work also opens a discussion on the temperature and density at the base of the corona of solar-like stars (see Eqs. (9) and (10)). An empirical relation between the coronal temperature, the mass, and the stellar rotation rate was proposed by Johnstone & Güdel (2015), who studied the correlation between the coronal temperature and the X-ray flux of low-mass MS stars. From scaling laws and the relation between the X-ray luminosity and the stellar rotation period found by Reiners et al. (2014), they derived the expression which differs from what we used (Holzwarth & Jardine 2007). Concerning the density at the base of the corona, the hypothesis of Holzwarth & Jardine (2007) was justified by coronal electron density measurements from Ivanova & Taam (2003). Both the reliability ofthese observations and the possibility of inferring the base density of the wind from such quantities were criticized (on the former point, see Güdel 2004). However, more recent stellar wind models do not contradict this law (Cranmer & Saar 2011; Suzuki et al. 2013). In a future work, we will consider changing the assumptions we made on coronal parameters for relations taking into account the stellar structure and evolution. However, such developments are beyond the scope of the present paper.

Another important point of improvement is the Rossby number. In this paper, we used the expression developed by Cranmer & Saar (2011) for their models of MS solar-like stars to express the Rossby number and used the value of T_eff at the ZAMS in the formulation. However, convective turnover times are significantly longer for young stars than for their MS counterparts. Thus, even slowly rotating PMS stars may be in a saturated regime, which has been observed by recent X-ray surveys on young clusters (Getman et al. 2005; Güdel et al. 2007). This suggests that only stars older than a certain age can reach the unsaturated regime. The limit age has to be of the order of 13 Myr, the age of the h Per cluster which was studied by Argiroffi et al. (2016). They found that it was the youngest known cluster that behaves as clusters of MS stars.

Recent works have started accounting for the variations of the Rossby number with stellar structure using stellar evolution tracks (Folsom et al. 2016; Sadeghi Ardestani et al. 2017). Both computed the convective turnover time at each time step as the ratio of the local pressure scale height H_p divided by the local convective velocity at one H_p over the base of the convective zone. These approaches are close to that of Spada et al. (2013), who published evolutionary tracks of solar-like stars for various compositions and mixing length parameters, allowing one to study the variations of the convective turnover timescale along stellar evolution. The peculiarity of these formulations is that they allow consistent treatment of the Rossby number over the PMS phase, which to our knowledge is not taken into account in other secular evolution models of rotation of solar-like stars. The question of the Rossby number was also extensively studied by Brun et al. (2017). In particular, they highlighted that the fluid Rossby number R_of, defined as the ratio of the advection term divided by the Coriolis force, was a good indicator to assess the large-scale rotational behavior of solar-like stars. Based on numerical simulations performed with the ASH code, they gave an expression of this number as a function of stellar mass and surface rotation rate.

In a future work, we will take into account the studies mentioned above to improve this part of our model. A better understanding of the magnetic and coronal properties of solar-like stars will also be useful in order to implement star–planet magnetic interactions in the code (Strugarek et al. 2014, 2015). Recently, Strugarek et al. (2017b) proved that this effect could compete with tidal dissipation in shaping star–planet systems.

Finally, in this work we only considered isolated star–planet systems. Future developments of our code will include implementing dynamical interactions in multiplanet systems (e.g., Laskar et al. 2012; Bolmont et al. 2015). Developing such physical models is important to extract the maximum amount of information from the data from upcoming space missions such as PLATO (Rauer et al. 2014).

Acknowledgements

The authors acknowledge Florian Gallet, Emeline Bolmont, Antoine Strugarek, and Rafael A. García for helpful discussion during the project. We also thank Jérémy Ahuir for providing Figs. C.3 and C.4. S.M. acknowledges funding by the European Research Council through ERC grant SPIRE 647383. We also thank the ERC STARS2 207430, FP7 SPACEINN, the ANR Toupies SIMI5-6 020 01, the INSU/PNST and CNES support via our PLATO and Solar Orbiter fundings. We thank the referee and editor for useful comments that improved the content of the paper.

Appendix A Description of the code

This appendix summarizes the key equations in the paper that are solved by the code. The purpose of ESPEM is to compute three quantities of interest over the life of the considered star–planet system: the semi-major axis a, the angular momentum of the shellular convective envelope L_c, and that of the inner radiative zone L_r. At each time step, the torques Γ_int, Γ_wind, and Γ_tide are computed from the state of the system. The time derivatives of the quantities of interest are then computed from these torques.

The computation of Γ_wind is different from that of the other two torques. To compute the stellar wind, we used the starAML code which implements the method of Réville et al. (2015b). We refer the reader to Sect. 3 for more details about this method. Since their technique implies to compute the radial profile of the wind and magnetic field in the whole stellar corona, we decided to precompute the torque of the wind rather than calling the starAML routine at each ESPEM time step because the latter option would have significantly slowed down our simulations.

The other two torques, Γ_int and Γ_tide, are directly computed from the state of the system and the stellar evolution variables. The latter are read from a precomputed STAREVOL evolutionary track (Amard et al. 2016).

Appendix B Comparison to other models

This work is similar to other studies of star–planet systems. In this section, we discuss the differences brought by the tidal prescription and the two-layer internal rotation model we implemented in ESPEM. To that end, we compare the results of these frameworks with those of state-of-the-art models.

B.1 Comparison with the constant quality factor model

In this paragraph, we emphasize the consequences of changing the tidal prescription. We compare the predictions made by the model described in Sect. 4 with those of the constant quality factor model, used for instance by Zhang & Penev (2014). For this purpose, we computed for each tidal prescription, the evolution of a system composed by a solar-mass star and a Jupiter-mass planet with an initial semi-major axis equal to 0.03 AU. The initial rotation period of the star was set to 1.4 days, which corresponds to the fast rotator described by Gallet & Bouvier (2015). We tested three values of Q′ for the constant quality factor prescription, 10⁶, 10⁷, and 10⁸. Here, we detail the case where Q′ = 10⁶. The two other scenarios are discussed at the end of this subsection.

Figure B.1 shows the comparison between the constant Q′ and variable Q′ models. As can be seen in the middle panel, the former case leads to the destruction of the planet whereas in the latter the planet survives. What allows the survival in the variable Q′ model is thefast outward migration during the first million years of the simulation. This phenomenon is due to the low value of the tidal quality factor in this period.

The color backgrounds correspond to the different evolution phases of the system in the variable Q′ case. The orange background stands for the disk-locking phase. The green background represents the outward-migration phase. During this stage, both the equilibrium and the dynamical tides are active. The pink background corresponds to the inward-migration phase during which both tides are raised. During this stage, the orbital period P_orb is such that P_*∕2 < P_orb < P_*, where P_* is the surface rotation of the star. Finally, the light pink background depicts the inward-migration phase during which only the equilibrium tide is active. During the latter two phases, the planet is too far from the star to undergo significant orbital migration.

This figure illustrates the difference between our tidal model and the constant Q′. Bolmont & Mathis (2016) performed simulations on a extensive parameter range to show the consistency of this difference.

B.2 Comparison with the one-layer rotation model

One of the main differences between our code and that of Bolmont & Mathis (2016) is the internal rotation model implemented. We used the two-layer model of MacGregor & Brenner (1991), in which the tidal and wind braking torques are applied on the convective envelope, which relaysthem to the core through a parametrized internal coupling torque with a constant timescale, whereas they assumethat the star is in solid body rotation. Consequently, stellar rotation in our model is expected to be more sensible to tidal dissipation due to the low inertia of the convective envelope. To investigate this, we measuredthe spin-down of a star caused by a close-in massive planet in orbit beyond the corotation radius. For both the one-layer and the two-layer models, we computed the secular evolution of the rotation period of a 0.6 M_⊙ star such as P_*,ini = 1.2 days, in a first time without planet, then with a 5 M_Jup planet initially orbiting at a distance a_ini = 0.024 AU. We then calculated the difference δP between the two cases as defined in Eq. (24).

Figure B.2 shows that both models predict similar evolutions. As can be seen in the top panel, the planet migrates outward during the first hundreds of millions of years. Consequently, the star is spun down, which is why δP is positive. In both cases, δP increases during the outward-migration phase, reaches a maximum of the order of one day around the age of 1 Gyr, and decreases over the MS phase. The curve of δP in the two-layer model is slightly delayed compared to that of the one-layer model. Studying this delay requires to perform simulations on an extensive parameter space, which is beyond the scope of this work.

Appendix C Discontinuities of the tidal torque

The formulation of the tidal torque in Eq. (19) implies discontinuities. This occurs at two events: when the planet orbits the star at a distance close to the corotation radius and when the orbital period is close to half the rotation period. In the former case, the tidal frequency is close to zero and the sign of the tidal torque may change with the absolute value remaining constant. In the latter case, the tidal frequency approaches − 2 Ω_*, the limit of application of the dynamical tide, where the torque may discontinuously vary by several orders of magnitude, to reach the values of the equilibrium tide.

Figure C.1 illustrates the different phases of the secular evolution of a star–planet system formed by a star of mass 1 M_⊙ and a planet of mass 1 M_Jup initially orbiting its host at 0.04 AU. The events at which the discontinuities of the tidal torque occur are marked by vertical dashed graylines. As is visible on the bottom panel, changes of tidal regime induce strong variations of the tidal torque.

These discontinuities generate numerical problems during the integration of differential equations. The Bulirsch–Stoer algorithm from Hairer et al. (2000) implements an adaptive step size to ensure the stability of the solution. Discontinuities in the derivatives make the step size shrink and prevent the resolution over the evolution. This is why we slightly modified the dependence of the tidal torque on ω_tide in the code.

To fix the null-tidal-frequency situation, we replaced the sign function by a hyperbolic tangent in Eq. (19). We introduced the characteristic pulsation ω₀ = 10⁻⁸ rad s⁻¹ such that the expression used in ESPEM is (C.1)

This modification makes the torque change sign continuously without altering its value for |ω_tide| > |ω₀|. For example, Fig. C.2 shows the tidal torque as a function of the tidal frequency around 0 for a system formed by a star of mass 1 M_⊙ and a planet of mass 1 M_Jup.

We performed a similar correction in the second case. For ω_tide < −2 Ω_*, the dynamical tide does not apply so only the equilibrium tide contributes to the orbital evolution. For ω_tide > −2 Ω_*, inertial waves are raised in the stellar envelope and both tides contribute to the interaction. Consequently, the theoretical dissipation can be written as¹

The discontinuity arises from the addition of the dissipation of the dynamical tide as soon as the tidal frequency passes the limit 2 Ω_*. To fix this problem, we replace this dependence with a more regular function based on a hyperbolic tangent: (C.4)

where is the regularized tidal quality factor and ω₁ = 10⁻¹⁰ rad s⁻¹. The factor multiplying is null when the inertial waves are not raised and it equals 1 when they apply. Consequently, this regularized dissipation is continuous and fits the theoretical dissipation when the tidal frequency is farther than ω₁ from − 2 Ω_*.

The choice of ω₁ is worth commenting upon. To this end, we computed the secular evolution of given star–planet systems for different values of ω₁ ranging from 10⁻¹⁰ to 10⁻⁶ rad s⁻¹. We did not try larger values for ω₁ because this parameter has to be small compared to Ω_*, which is of the order of Ω_⊙ on the MS for a solar-like star. Here, we illustrate the results of this experiment with the cases ω₁ = 10⁻¹⁰ rad s⁻¹ and ω₁ = 10⁻⁶ rad s⁻¹ for a 1 M_⊙ star and a 1 M_Jup planet with P_*,ini = 5 days and a_ini = 0.025 AU.

Figure C.3 shows the tidal torque in each case as a function of tidal frequency. Despite the fact that the torque is significantly smoothed for the larger value of ω₁, the evolution of the system is only weakly influenced by this parameter. This can be seen in Fig. C.4. Indeed, in both cases, the evolution phases of the system are the same. The lifetime of the planet is slightly shorter in the smoothed case but it is not significantly different from the other. We conclude from this study that the parameter ω₁ has a weak influence on the evolution. Therefore, we set its value to ω₁ = 10⁻¹⁰ rad s⁻¹, so as to stayclose to the theoretical tidal calculations. However, the excitation and disappearance of the dynamical tide should be discussed in a further work.

References

---

All Tables

All Figures