© J. Ahuir et al. 2021

1 Introduction

Since the detection of 51 Pegasi b by Mayor & Queloz (1995), more than 4000 exoplanets have been detected. The observed populations show a wide variety of host stars, orbital architectures, planetary sizes, and masses. Moreover, because of the biases of the most successful detection methods, namely, the transit and radial velocity techniques, a majority of the discovered planets orbit close to their host stars, whether they are of mass comparable to that of Jupiter (forming the population of hot Jupiters, e.g., Mayor & Queloz 1995; Henry et al. 2000; Charbonneau et al. 2000) or slightly larger than that of the Earth (the so-called super-Earths, such as 55 Cnc e; see Dawson & Fabrycky 2010). Close-in planets orbit in a dense and magnetized medium, which leads to the emergence of star–planet interactions that can affect the dynamics and evolution of the orbital systems (Cuntz et al. 2000). In particular, angular momentum exchanges can occur between the planetary orbit and the stellar spin, leading to migration of the planet. Potential signatures of these interactions may have been identified in individual systems (e.g., HD 189733 ; see Dowling Jones et al. 2018; Cauley et al. 2018) as well as in the distribution of some planetary populations. More precisely, McQuillan et al. (2013) estimated the rotation period of 737 stars hosting Kepler objects of interest using an auto-correlation method, and identified a possible dearth of planets with orbital periods shorter than 2–3 days around fast rotators (with a rotation period shorter than 10 days). Teitler & Königl (2014) first proposed that such a phenomenon may be attributed to the engulfment of close-in planets by their host stars through tidal interactions. Lanza & Shkolnik (2014) suggested an alternative scenario based on secular perturbations in multiplanet systems. These latter authors showed that remote planets which are excited on a sufficiently eccentric orbit around old stars may be tidally circularized on shorter orbits. Furthermore, Walkowicz & Basri (2013) found a concentration of massive planets with an orbital period equal to either the rotation period of their host star or half that period, which could be the signature of tidal interactions. In view of these different aspects, understanding how compact systems form and evolve is a key astrophysical question to be addressed.

The role of the protoplanetary disk in shaping the observed structure of planetary systems is strongly emphasized in the literature. Indeed, the disk structure and evolution have a significant influence on the mass and semi-major axis distribution of the young planets (e.g., Mordasini et al. 2009a,b). In particular, planet migration in the disk through Lindblad resonances is thought to be efficient in shaping planetary systems (Baruteau et al. 2014; Bouvier & Cébron 2015; Heller 2019). Moreover, population synthesis models have been developed to better understand the interplay between the properties of the disk and the different processes shaping planetary systems (we refer the reader to Mordasini 2018, for an extended review). The predicted distributions were then compared to Kepler observations in order to constrain models of planetary formation and evolution (Mulders et al. 2019). For a large range of multi-planet system properties (e.g., orbital period ratios, mutual inclination, position of the innermost planet), these synthetic populations are in good agreement with Kepler global distributions, if multiple interacting seed planet cores per disk are taken into account.

After the dissipation of the disk, dynamical interactions, in particular Kozai oscillations, may occur in multiplanet systems, leading to an intricate evolution of their orbital architecture (Laskar et al. 2012; Bolmont et al. 2015). However, isolated close-in planets can already suffer efficient migration because of magnetic and tidal interactions with their host star. We consider this simpler case in the present work, where we aim to account for the variety of such star–planet interactions. We therefore consider a simplified system comprising a star and a single planet on a circular orbit perpendicular to the stellar rotation axis. One of the main physical processes acting in such a configuration are tidal interactions. These result from the gravitational response of a star to the presence of a planet and play a key role in the evolution of the orbital configuration of the system. Two components arise from the stellar response: the hydrostatic nonwavelike equilibrium tide, dissipated by turbulent friction (Zahn 1966; Remus et al. 2012; Ogilvie 2013), and the dynamical tide, which consists in the excitation of waves inside the star by the tidal potential as well as their dissipation. A dynamical tide can exist in both radiative zones (through internal gravity waves, Zahn 1975; Goldreich & Nicholson 1989; Goodman & Dickson 1998; Terquem et al. 1998) and convective zones (Ogilvie & Lin 2007; Ogilvie 2013; Mathis 2015). We focus on the latter in this study and leave a detailed investigation of internal gravity waves to future work (Barker 2020; Ahuir et al. 2021). In the stellar convective zone, if the orbital period of the planet is longer than half of the stellar rotation period, inertial waves restored by the Coriolis force are excited and dissipated in the envelope of solar-type stars. Otherwise, inertial waves cannot be excited and no dynamical tide is raised in the star (see the violet area in Fig. 1). The associated dissipation, which depends on stellar internal structure as it arises from the reflection of the waves on the radiative core (Ogilvie 2013; Goodman & Lackner 2009; Mathis 2015), can be several orders of magnitude higher than the dissipation of the equilibrium tide (Ogilvie & Lin 2007; Bolmont & Mathis 2016). When tidal dissipation is taken into account, the stellar response presents a delay and the associated bulge is misaligned with the line joining the centers of the two celestial bodies. This misalignment then induces a lag angle δ, which increases with dissipation magnitude, and is at the origin of an exchange of angular momentum between the star and the planet. Indeed, a net tidal torque is applied both on the planetary orbit and on stellar rotation to reduce the angle δ. The position of the planet with respect to the co-rotation radius (for which n = Ω_⋆, see the black dashed line in Fig. 1) then determines the evolution of the system. If the planet is situated beyond this characteristic distance, an outward migration makes it move away from the host star, which then spins down. Otherwise the planet migrates inward, moving closer to the spinning-up star. In the latter case, the fate of the system is determined by the orbital angular momentum L_orb of the planet and the stellar angular momentum L_⋆: if L_orb ≥ 3 L_⋆, an equilibrium state is reached where the angular velocities Ω_⋆ and n are synchronized. Otherwise, the planet migrates too efficiently and it spirals towards its host star until its disruption at the Roche limit (Hut 1980, 1981). It is important to note that such a condition is based on the conservation of the total angular momentum of the star–planet system. Damiani & Lanza (2015) derived a similar criterion by taking into account magnetic braking. More generally, to provide a realistic and complete equilibrium criterion, it is necessary to take into account magnetic braking, angular momentum redistribution within the star, and the various star–planet interactions simultaneously.

Tidal interactions are modulated by the structural and rotational evolution of the star. The stellar rotation rate generally exhibits a complex evolution due to the initial disk–star interaction, the internal redistribution of angular momentum within the star during contraction phases, and the angular momentum extraction by the stellar wind (e.g. Weber & Davis 1967; Skumanich 1972; Kawaler 1988; Matt et al. 2012; Réville et al. 2015a). As this intricacy may significantly affect the evolution of a given star–planet system, it is necessary to take such processes into account as much as possible. This requires models calibrated on gyrochronology and rotational distributions in open clusters (Gallet & Bouvier 2015; Matt et al. 2015), leading up to general frameworks combining stellar rotation, wind, and magnetism (Johnstone et al. 2015; Ahuir et al. 2020). Past studies have taken these constraints into account to some extent. For instance, Zhang & Penev (2014) relied on the two-layer rotational model of MacGregor & Brenner (1991) and a constant tidal dissipation to deal with the secular evolution of star–planet systems, and subsequently adopted a statistical approach to their numerical simulations in order to apply constrains to tidal theory. Bolmont & Mathis (2016) then first incorporated the effects of dynamical tide in stellar convection zones and studied their impact on the secular evolution of star–planet systems by considering a one-layer rotational model for the central star. More recently, Benbakoura et al. (2019) performed a study based on the amalgamation of the two previous studies, taking a bi-layer structure for the star and both the equilibrium and dynamical tides into account. This allowed them to provide a criterion for planetary engulfment due to tidal effects, taking into account stellar evolution. They were also able to characterize the influence of such a phenomena on the rotation of the host star. In parallel, Gallet et al. (2018) and Gallet & Delorme (2019) relied on a similar model to investigate the rotational evolution of planet-hosting stars in open clusters (in particular in the Pleiades) in more detail, which allowed them to assess some limits of gyrochronology (Barnes 2003). Finally, the evolution of star–planets system under tidal interactions during the red giant phase has also been extensively investigated (Privitera et al. 2016a,b; Meynet et al. 2017; Rao et al. 2018).

However, in past studies, tidal effects and magnetism have not been taken into account together systematically (apart from wind braking). Bouvier & Cébron (2015) first explored the possibility that tidal and magnetic interactions may compete with accretion and contraction in the case of a close-in planet embedded in a disk. Furthermore, after the dissipation of the latter, star–planet magnetic interactions may occur because of the relative motion between the planet and the ambient wind at the planetary orbit (represented with orange arrows in Fig. 1). If the planet is below the Alfvén radius (at which the wind velocity is equal to the local Alfven speed; see the blue line in Fig. 1), the magnetic torque applied to the planet can lead to efficient transport of angular momentum between the planet and the star through the so-called Alfvén wings (Neubauer 1998, see the orange lobes around the planet in Fig. 1). Beyond the Alfvén radius, the wind becomes superalfvénic. In this case, Alfvén wings may still exist but do not connect back the star (see the gray region in Fig. 1). In this case, the planet may transfer energy and angular momentum to the ambient wind instead. In the context of close-in planets, we only consider here the subalfvenic scenario. Several regimes then appear depending on the star–planet configuration (Strugarek 2017). If Alfvén waves have enough time to go back and forth between the star and the planet before the magnetic field lines slip through the planet, the interaction acts as a unipolar generator, leading to the so-called unipolar interaction (Laine et al. 2008; Laine & Lin 2012). In the opposite case, magnetic interactions between the planet and the star still occur and the interaction becomes dipolar (Saur et al. 2013; Strugarek et al. 2015; Strugarek 2016). As the relative motion between the planet and the ambient wind is at the origin of the subsequent magnetic torques, these are then expected to act in the same way as tidal effects in most cases. The co-rotation radius then plays a determining role in planetary migration in both cases. Many other notable effects, such as anomalous emissions or planet inflation, may result from star–planet magnetic interactions (we refer the reader to Lanza 2018, for a recent review). Strugarek et al. (2017) performed a first study on planetary migration taking into account tidal and magnetic torques simultaneously. In particular, they computed the migration timescale of the planet for both contributions, finding that both effects could play a key role depending on the characteristics of the star–planet system considered. Thus, following the orbital evolution of a planet along the rotational and structural evolution of the host star by taking into account the coupled effects of tidal and magnetic torques is essential to better understanding the evolution of star–planet systems. Furthermore, synthesizing planetary populations by taking into account the whole variety of star–planet interactions to explain the observed distributions of exoplanets still has to be performed.

Following Strugarek et al. (2017) and Benbakoura et al. (2019), the main goal of this work is to assess the relative contributionof both tidal and magnetic interactions on the secular evolution of star–planet systems, and to investigate the role of the associated torques in shaping the distributions of planetary populations. In Sect. 2, we present the hypotheses of our study, detail the interactions involved in our modeled star–planet systems, and describe the modeling approach used in this work. In Sect. 3, we investigate the influence of the main characteristics of a star–planet system (e.g., stellar mass, stellar magnetism, semi-major axis, planetary type, etc.) on its secular evolution by assessing the relative contribution of magnetic and tidal torques. All these parameters are then taken into account simultaneously in Sect. 4 in order to classify planetary populations emerging from the action of star–planet interactions and to highlight regions of interest resulting from their evolution. Populations of star–planet systems are then synthesized in Sect. 5 and are confronted with a statistical distribution obtained from Kepler data. All of those results are then summarized, discussed, and put into perspective in Sect. 6.

Fig. 1 Sketch of the main features and locations of interest involved in a star–planet system undergoing tidal and magnetic interactions. The planet (in blue) orbits the star (in orange) with an orbital angular velocity n. As a result of the presence of a planet, the star presents a bulge misaligned with the line joining the centers of the two celestial bodies, which induces a lag angle δ. This angle has been greatly exaggerated here for visualization purposes (it is indeed much smaller than 1 degree in almost all cases). If the planetary mean motion is greater than twice the stellar rotation rate, inertial waves cannot be excited in the stellar convective zone and no dynamical tide is raised within the convective envelope of the star (see the purple area). The relative motion between the planet and the ambient wind (represented with orange arrows) leads to the formation of Alfvén wings (orange lobes around the planet) if the planet is below the Alfvén radius (in blue). Beyond thisdistance, no Alfvén wings can connect the star and the planet (see the gray area). For both tidal and magnetic interactions, a planet situated below the co-rotation radius (see the black dashed line) undergoes an inward migration, and a planet situated beyond this distance migrates outwards. The relative position of the different orbits of interest mayvary depending on the initial configuration of the system considered.

2 Star–planet interaction model

2.1 ESPEM: an overview

ESPEM (French acronym for Evolution of Planetary Systems and Magnetism; see Benbakoura et al. 2019) is a numerical code computing the secular evolution of a star–planet system by following the semi-major axis of the planetary orbit as well as the stellar rotation rate. Furthermore, we assume here a coplanar and circular orbit, and a synchronized planetary rotation, as the reservoir of angular momentum of the planet is less important than the one in its orbit (Guillot et al. 1996). In this model, we consider a two-layer solar-type star composed of a radiative core and a convective envelope. Tidal dissipation is only considered in the stellar envelope in this work. The core is interacting with the envelope through internal coupling, and the latter exchanges angular momentum with the orbit through tidal and magnetic interactions (SPMI). Moreover, the whole system loses angular momentum through magnetic braking by the stellar wind. Hence, the angular momentum of the planetary orbit, L_orb, the stellar convective zone, L_c, and radiative zone, L_r, are evolved by the following system of equations: $$\frac{\text{d}{L}_{\text{orb}}}{\text{d}t}=-{\Gamma }_{\text{tide}}-{\Gamma }_{\text{mag}}$$(1) $$\frac{\text{d}{L}_{\text{c}}}{\text{d}t}={\Gamma }_{\text{int}}+{\Gamma }_{\text{tide}}-{\Gamma }_{\text{wind}}+{\Gamma }_{\text{mag}}$$(2) $$\frac{\text{d}{L}_{\text{r}}}{\text{d}t}=-{\Gamma }_{\text{int}},$$(3)

where Γ_int is the internal torque, coupling the core and the envelope of the star, Γ_wind is the wind-braking torque, and Γ_tide and Γ_mag are the tidal and MHD torques between the star and the planet, respectively. A schematic global view of the system studied by ESPEM is provided in Fig. 2.

2.2 Stellar structure and evolution

Stellar structure and evolution are taken into account in ESPEM during the pre-main sequence (PMS) and the main sequence (MS). The internal structure of the star, especially the radii, the masses, and the moments of inertia of radiative and convective zones are provided at each ESPEM time-step through grids precomputed with the stellar evolution model STAREVOL (Siess et al. 2000; Palacios et al. 2006; Amard et al. 2016, 2019).

The coupling between the radiative zone and the convective zone following the model proposed in MacGregor & Brenner (1991) is taken into account as an exchange of angular momentum allowing the synchronization of their spins on a characteristic timescale τ_c−e, which is determined by internal transport processes in the radiative core (Brun & Zahn 2006; Mathis 2013; Aerts et al. 2019) as well as the effective coupling between the radiative and the convective zones (Brun et al. 2011; Strugarek et al. 2011). The amount of angular momentum to be transferred between the two layers of the star to equilibrate their angular velocities can be expressed as $$\Delta L=\frac{I_{\text{r}}{I}_{\text{c}}}{I_{\text{r}}+I_{\text{c}}}\left({\Omega }_{\text{r}}-{\Omega }_{\text{c}}\right),$$(4)

where I_r, Ω_r are the moments of inertia and the rotation rate of the core, and I_c,  Ω_c are the same quantities assessed in the envelope.

Moreover, the expansion of a radiative core during the PMS involves a rapid conversion of convective state to radiative state (Emeriau-Viard & Brun 2017). During this transition phase, a significant mass transfer occurs, which is accompanied by a transport of angular momentum. This way, the coupling between the core and the envelope can be modeled as a torque with two components applied to the convective zone: $${\Gamma }_{\text{int}}=\frac{\Delta L}{{\tau }_{c-e}}-\frac{2}{3}{R}_{\text{r}}{}^2{\Omega }_{\text{c}}\frac{\text{d}{M}_{\text{r}}}{\text{d}t},$$(5)

where M_r and R_r are the mass and radius of the radiative core. The coupling timescale τ_c−e is a free parameter of our model and has been calibrated with the Gallet & Bouvier (2015) study as follows: $${\tau }_{\text{c--e}}(\text{Myr})=3.05{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{-5.02}{\left(\frac{P_{\text{rot,c}}}{P_{\text{rot,}\odot }}\right)}^{0.67},$$(6)

where M_⋆ is the stellar mass and P_rot,c is the rotation period of the stellar convective zone. Star–disk interaction is taken into account in a simplified way at the beginning of the PMS by assuming a constant surface rotation rate during the disk’s lifetime, which is also fixed by the Gallet & Bouvier (2015) study as $${\tau }_{\text{disk}}[\text{Myr}]=13.41{\left(\frac{P_{\text{rot,c}}}{P_{\text{rot,}\odot }}\right)}^{-0.56}.$$(7)

During this phase, the semi-major axis of the planet is assumed to be constant and the rotation of the radiative zone is only constrained through internal coupling. We focus here on the evolution after the disk dissipation. A more precise treatment of the early phase could be added in future works (Bouvier & Cébron 2015; Gallet et al. 2019).

Fig. 2 Schematic view of the system and its interactions (adapted from Benbakoura et al. 2019). The radiative core (in yellow) and the convective envelope (in orange) exchange angular momentum (green arrows). Stellar wind carries away angular momentum from the envelope and spins the star down (purple arrows). Stellar rotation and planetary orbit are coupled through tidal (red arrows) and magnetic effects (blue arrows).

2.3 Wind braking torque

The wind braking torque is given in our model by Matt et al. (2015): $${\Gamma }_{\text{wind}}={\Gamma }_{\odot }{\left(\frac{R_{\star }}{R_{\odot }}\right)}^{3.1}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{0.5}{\left(\frac{Ro}{R{o}_{\odot }}\right)}^{-2}\left(\frac{{\Omega }_{\text{c}}}{{\Omega }_{\odot }}\right),Ro>R{o}_{\text{sat}},$$(8) $${\Gamma }_{\text{wind}}={\Gamma }_{\odot }{\left(\frac{R_{\star }}{R_{\odot }}\right)}^{3.1}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{0.5}{\left(\frac{R{o}_{\text{sat}}}{R{o}_{\odot }}\right)}^{-2}\left(\frac{{\Omega }_{\text{c}}}{{\Omega }_{\odot }}\right),Ro\le R{o}_{\text{sat}},$$(9)

with Γ_⊙ = 6.3 × 10³⁰ erg. Such a prescription allows us to account for the mass and age dependencies of the distribution of stellar rotation periods in open clusters (Rodríguez-Ledesma et al. 2009; Agüeros et al. 2011) and in the sample of stars observed by the Kepler satellite (McQuillan et al. 2014). We use the stellar Rossby number for simplicity purposes, expressed as $$Ro=\frac{P_{\text{rot,c}}}{{\tau }_{\text{c}}},$$(10)

where P_rot,c is the rotation period of the stellar convective zone (see Landin et al. 2010; Brun et al. 2017, for a discussionon the various definitions found in the literature). The convective turnover time τ_c can be assessed with the Sadeghi Ardestani et al. (2017) prescription: $$\begin{array}{ll}{\log }_{10}{\tau }_{\text{c}}[\text{s}]=\hfill & 8.79-2|{\log }_{10}(m_{\text{CZ}}){|}^{0.349}-0.0194{\log }_{10}^2(m_{\text{CZ}})\hfill \\ \hfill & -1.62\text{min}\left[{\log }_{10}(m_{\text{CZ}})+8.55,0\right],\hfill \\ \hfill & \hfill \end{array}$$(11)

where m_CZ is the mass of the convective envelope normalized to the stellar mass. The formulation obtained by these latter authors has the advantage of being valid during the pre-main sequence and the main sequence for metallicities ranging from [Fe/H] = −0.5 to 0.5 and was obtained with the CESAM stellar evolution code (Morel & Lebreton 2008). Their prescription leads to a solar value of Ro_⊙ = 1.113 and a saturation value of Ro_sat = 0.09.

Figure 3 shows the typical rotational evolution obtained with our model for isolated stars. We show our results for three stellar masses, and for three initial rotation periods (1.4, 5, and 8 days) covering fast, median, and slow rotators from Gallet & Bouvier (2015). The initial rotation spread is reduced over the MS as all the stars converge towards a sequence where their rotation rate is fully determined by their age and mass (Barnes 2003). As seen in the top panel, solar-mass stars spin down to reach the solar rate at the solar age whereas less-massive stars reach a lower rotation at the same age. A steeper evolution of the stellar rotation rate compared to the Skumanich law occurs for each stellar mass because of the core-envelope coupling, in accordance with the Gallet & Bouvier (2015) results, as the angular momentum stored in the radiative zone is redistributed on secular timescales.

Fig. 3 Secular evolution of the rotation rate of the convective envelope (solid lines) and of the radiative core (dashed lines) for stars with a mass M⋆ = {1,  0.8,  0.5} M⊙ (from top to bottom). Slow (orange), median (light blue), and fast (dark blue) rotators are considered. The solar rotation rate at solar age is represented by a black circle in each panel. The dots with error bars correspond to the 25th, 50th, and 90th percentiles of rotational distributions of observed stellar clusters published by Gallet & Bouvier (2015). In black: Skumanich law normalized to the Sun.

2.4 Planetary properties

We consider in our model a planet of mass M_p, ranging between 0.5 and 1589 M_⊕ (corresponding to 5 M_Jup), and a radiusR_p. It is assumed to be punctual and its rotation is synchronized with its orbit. Furthermore, we adopt the probabilistic mass–radius relations proposed by Chen & Kipping (2017) based on a sample of well-constrained planets: $$R_{\text{p}}\propto \{\begin{array}{ll}{M}_{\text{p}}^{0.28},\hfill & M_{\text{p}}<2.0M_{\oplus }(6.29\times {10}^{-3}{M}_{\text{Jup}})\hfill \\ M_{\text{p}}^{0.59},\hfill & 2.0M_{\oplus }\le M_{\text{p}}<0.4M_{\text{Jup}}\hfill \\ M_{\text{p}}^{-0.04},\hfill & M_{\text{p}}\ge 0.4M_{\text{Jup}}.\hfill \end{array}$$(12)

Recent studies show that the distribution of planetary radii in the Kepler sample presents a gap between 1.5 R_⊕ and 2 R_⊕ (Fulton et al. 2017). Such a bimodality in the distribution may be due to photoevaporation which may drive atmospheric mass loss on close-in planets. If so, the gap would originate from a discrepancy between planets with H/He envelopes of small mass and bare rocky cores. In this first statistical study, as such a feature does not appear in the star–planet sample we have considered (we refer the interested reader to Sect. 5.1 for more details), the incorporation of this radius valley in our model is left for future work.

The equatorial field at the planetary surface B_p, which is of prime importance in star–planet interactions, is by default assumed to be constant. Nevertheless, this is varied in Sect. 3.3.

2.5 Tidal effects

Tidal interactions lead to an angular momentum exchange between the star and the planet, which can be translated into a tidal torque applied to the stellar envelope (Murray & Dermott 1999; Benbakoura et al. 2019): $${\Gamma }_{\text{tide}}=-\text{sign}({\omega }_{\text{tide}})\frac{9}{4Q^{\prime }}\frac{G{M}_{\text{p}}^2}{a^6}{R}_{\star }^5,$$(13)

where M_p is the planetary mass, R_⋆ the stellar radius, G the gravitational constant, a the semi-major axis, ω_tide = 2(Ω_c − n), which is the tidal frequency in the case of a planet with a circular orbit and a synchronized rotation, and n is the mean motion of the planetary orbit. The equivalent quality factor Q′ takes into account the nature and the efficiency of tidal dissipation as a function of the internal structure and rotation of the star. Here, it is used to describe the so-called equilibrium (Zahn 1966; Remus et al. 2012) and dynamical tides (Ogilvie 2013; Mathis 2015). We now summarize their treatment (we refer the reader to Benbakoura et al. 2019, for a more detailed description).

The equilibrium tide is taken into account in ESPEM by following the Hansen (2012) prescription, relying on a constant value for the dimensionless dissipation factor ${\overline{\sigma }}_{\star }$ along the evolution of the system. This quantity is calibrated using observations and leads to the quality factor (Bolmont & Mathis 2016) $${Q^{\prime }}_{\text{eq}}=\frac{3}{2}\frac{1}{{\sigma }_0{\overline{\sigma }}_{\star }}\frac{G}{R_{\star }^5\left|{\omega }_{\text{tide}}\right|},$$(14)

where ${\sigma }_0=\sqrt{G/(M_{\odot }{R}_{\odot }^7)}$. In this work, the value of ${\overline{\sigma }}_{\star }$, which decreases with stellar mass, is taken from Fig. 3 of Hansen (2012). Such a formulation provides the same order of magnitude as prescriptions derived from physical models (for two different approaches, see Strugarek et al. 2017; Barker 2020).

The dynamical tide, and more precisely the dissipation of tidal inertial waves (governed by the Coriolis acceleration) within the convective zone, is based on the prescription of Ogilvie (2013) and Mathis (2015), who introduced a frequency-averaged effective constant tidal quality factor. Indeed, when computing the frequency dependence of the tidal torque due to tidal inertial waves in stellar convection zones (see e.g., Ogilvie & Lin 2007), its frequency dependence is highly resonant and erratic. This complex behavior relies on the physics of the friction applied by the turbulent convection on tidal inertial waves (e.g., Ogilvie & Lin 2004; Auclair-Desrotour et al. 2015); works are ongoing to improve its complex modeling (Duguid et al. 2020). Therefore, as explained in Sect. 4.2 of Benbakoura et al. (2019), a consistent treatment of the frequency dependence of the torque induced by tidal inertial waves would require a coupling of ESPEM with 2D hydrodynamical computation of tidal inertial modes at each time-step, which is beyond the scope of this work and would also not allow us to explore the broad parameter space describing the diversity of star–planet systems.

To get an order of magnitude of this dissipation for a given stellar mass, age, and rotation, we perform the frequency average of the dissipation as proposed and described in Ogilvie (2013), Mathis (2015), and Barker (2020) which provides us with good trends when compared with observational constraints. In the case of a two-layer star, the formulation of the frequency-averaged tidal dissipation $\overline{{Q^{\prime }}_{\text{dyn}}}$ ¹ provided by these latter authors gives $$\begin{array}{ll}\frac{3}{2\overline{{Q^{\prime }}_{\text{dyn}}}}=\hfill & \frac{100\pi }{63}{\left(\frac{{\Omega }_{\text{c}}}{{\Omega }_{\text{crit}}}\right)}^2\left(\frac{{\alpha }^5}{1-{\alpha }^5}\right){(1-\gamma )}^2{(1-\alpha )}^4\hfill \\ \hfill & \times \frac{{\left(1+2\alpha +3{\alpha }^2+\frac{3}{2}{\alpha }^3\right)}^2\left[1+\left(\frac{1-\gamma }{\gamma }\right){\alpha }^3\right]}{{\left[1+\frac{3}{2}\gamma +\frac{5}{2\gamma }\text{}\left(1+\frac{1}{2}\gamma -\frac{3}{2}{\gamma }^2\right){\alpha }^3-\frac{9}{4}(1-\gamma ){\alpha }^5\right]}^2},\hfill \end{array}$$(15)

where $\alpha =\frac{R_{\text{r}}}{R_{\star }}$, $\beta =\frac{M_{\text{r}}}{M_{\star }}$, and $\gamma =\frac{{\alpha }^3(1-\beta )}{\beta (1-{\alpha }^3)}$. The latter quantity corresponds to the ratio of the density of the envelope to that of the core. ${\Omega }_{\text{crit}}=\sqrt{G{M}_{\star }/R_{\star }^3}$ is the critical angular velocity of the star. Such a contribution will affect the secular evolution of the system if inertial waves are likely to be excited by the tidal potential, i.e., $P_{\text{orb}}>\frac{1}{2}{P}_{\text{rot}}$ (Bolmont & Mathis 2016). The total quality factor Q′, accounting for the sum of the tidal dissipation, is then given by: $$\frac{1}{Q^{\prime }}=\frac{1}{{Q^{\prime }}_{\text{eq}}}+\frac{1}{\overline{{Q^{\prime }}_{\text{dyn}}}}.$$(16)

2.6 Magnetic star–planet interactions

When a planet orbits in a magnetized medium, MHD disturbances propagate away from the planet vicinity while transporting energy and angular momentum, forming the so-called Alfvén wings (Neubauer 1998; Saur et al. 2013). Two Alfvén wings are always produced. Depending on the magnetic topology and the alfvenic Mach number of the interaction, one, both, or neither of the two may reach back to the star (for more details, see Strugarek et al. 2015). If Alfvén waves do not have enough time to travel back and forth between the star and the planet before the magnetic field lines slip through the planet, the magnetic interaction is dubbed dipolar (Saur et al. 2013; Strugarek et al. 2015; Strugarek 2016). Otherwise the interaction becomes unipolar (Laine et al. 2008; Laine & Lin 2012).

The existence of Alfvén wings results in a magnetic torque applied to the planet. Assuming the planet possesses a magnetosphere, it can be written as a drag torque (Strugarek 2016): $${\Gamma }_{\text{mag}}=-\text{sign}({\omega }_{\text{mag}})c_{\text{d}}\left(A_0{M}_{\text{a}}^{\beta }{\Lambda }_P^{\alpha }\pi R_{\text{p}}^2\right)p_{\text{tot}}a,$$(17)

where $c_{\text{d}}\approx M_{\text{a}}/\sqrt{1+M_{\text{a}}^2}$ is a drag coefficient representing the efficiency of the magnetic reconnection between the wind and the planetary magnetic fields, M_a is the alfvenic Mach number in the frame rotating with the planet, and ${\Lambda }_P=(B_{\text{p}}^2/2{\mu }_0)/p_{\text{tot}}$ is the ratio between the planetary magnetic pressure and the total pressure of the ambient wind at the planetary orbit. ω_mag corresponds to the difference between the rotation rate of the ambient wind and the planetary mean motion. As the wind is in a first approximation co-rotating with the star below the Alfvén radius, one can assume that ω_mag and ω_tide have the same sign in the vast majority of cases. In the case of close-in planets, the total wind pressure can be approximated by the magnetic pressure of the wind (Réville et. al. 2015b). The quantities A₀, α, β are determined from a set of 3D MHD simulations in Strugarek (2016). Only the case of a planetary dipole aligned with the stellar magnetic field at the planetary orbit is considered from now on in our model. Such a configuration, which maximizes the torque, gives A₀ = 10.8,  α = 0.28,  β = −0.56. This way, the considered torque provides an upper bound of the influence of magnetic effects on planetary migration in the dipolar regime.

The dipolar interaction regime considered here is likely to be realized in most compact star–planet systems. Generally, for planets sustaining a magnetosphere against the ambient pressure, alfvenic perturbations do not have enough time to travel back and forth between the star and the planet before the magnetic field line slips around the planet, unless the planetary magnetosphere of the planet is sufficiently large (comparable to the size of the Sun; see Strugarek 2017, for an extensive review of star–planet magnetic interactions). In the case of a weakly magnetized planet, the time-dependent component of the stellar magnetic field is either dissipated in the planetary interior or screened by the magnetic field induced by large surface currents, depending on the planetary resistivity. Therefore, we only consider the steady component of the stellar magnetic field. In this configuration, if the planetary diffusivity is sufficiently high, the time-independent component of the stellar magnetic field is efficiently dissipated inside the planet, creating a true magnetic cavity in the planetary interior. The dipolar regime is found to generally hold in this case. Otherwise, if magnetic diffusivity is sufficiently low, the magnetic field lines are frozen in the planet interior and dragged with the orbital motion of the planet. In that case, propagating Alfvén waves can generally reach back to the planet, and the interaction becomes unipolar. Such a configuration has been extensively treated by Laine et al. (2008), Laine & Lin (2012) and is found to lead to far stronger magnetic torques than for the dipolar regime (typically 4 or 5 orders of magnitude; see Strugarek et al. 2017). More complex situations can occur depending on the conductive properties of the planet material and its degree of ionization. As the transition between the unipolar and the dipolar regimes is still poorly understood, we focus in this work on the dipolar regime, which is more likely to occur in exosystems, and leave the study of the unipolar regime to a future study.

In this context, if the planet is not able to sustain a magnetosphere, we consider in this work that it screens the surrounding wind magnetic field, which leads to a dipolar star–planet interaction. The effective area of the planetary obstacle then corresponds to the geometrical cross-section of the planet. Henceforth, we rely on the following prescription for the dipolar torque: $${\Gamma }_{\text{mag}}=\{\begin{array}{l}\frac{-\text{sign}({\omega }_{\text{mag}})M_{\text{a}}}{\sqrt{1+M_{\text{a}}^2}}\left(10.8M_{\text{a}}^{-0.56}{\Lambda }_P^{0.28}\right)\pi R_{\text{p}}^2{p}_{\text{tot}}a,{\Lambda }_P>1.\hfill \\ \frac{-\text{sign}({\omega }_{\text{mag}})M_{\text{a}}}{\sqrt{1+M_{\text{a}}^2}}\pi R_{\text{p}}^2{p}_{\text{tot}}a,\text{otherwise}.\hfill \\ [-20pt]\hfill \end{array}$$(18)

Star–planet magnetic interactions occur because of the relative motion between the planet and the ambient wind at the planetary orbit. We must therefore estimate the radial profiles of the main characteristics of the wind, such as its velocity or its density. To this end, weincorporate a 1D isothermal magnetized wind model in ESPEM (see Lamers & Cassinelli 1999; Preusse et al. 2005; Johnstone 2017, for an extensive description of the model). Such a modeling requires knowledge of the temperature T_c and density n_c at the base of the wind. For consistency with the observational constraints on stellar rotation, wind, and magnetism, we rely on the Ahuir et al. (2020) prescriptions for those quantities. More precisely, as the stellar magnetic field measured from Zeeman broadening and Zeeman-Doppler imaging (see Montesinos & Jordan 1993; Vidotto et al. 2014; See et al. 2017) have only exhibited linear or super-linear dependencies between the large-scale magnetic field and the Rossby number, we consider for the sake of simplicity the following scaling law to assess the magnetic field at the stellar surface B_⋆ (Ahuir et al. 2020): $$B_{\star }[\text{G}]=2.0{\left(\frac{Ro}{R{o}_{\odot }}\right)}^{-1}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{-1.76},Ro>R{o}_{\text{sat}}.$$(19)

This leads to the following expressions for the coronal properties (Ahuir et al. 2020): $$T_{\text{c}}[\text{MK}]=1.5{\left(\frac{Ro}{R{o}_{\odot }}\right)}^{-0.11}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{0.12},Ro>R{o}_{\text{sat}},$$(20) $$n_{\text{c}}[{\text{cm}}^{-3}]=7.25\times {10}^7{\left(\frac{Ro}{R{o}_{\odot }}\right)}^{-1.07}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{1.97},Ro>R{o}_{\text{sat}}.$$(21)

Stellar magnetic field as well as wind temperature and density are assumed to be independent of the Rossby number in the rotation-saturated regime (Ro ≤ Ro_sat).

For the sake of simplicity as well as to provide an upper bound on magnetic effects in the dipolar regime, we assume that the planet is located in an open field region. The stellar magnetic field is then assumed to be radial. For more complex topologies, the magnetic field decays faster with distance to the star, which reduces the efficiency of the star–planet magnetic interactions accordingly.

3 Tidal and magnetic interactions: an evolutive approach of star–planet systems

3.1 Outline of star–planet secular evolution

3.1.1 Planet migration: reference case

We now aimto investigate the influence of the main properties of a star–planet system on its secular evolution by assessing the relative contribution of tidal and magnetic torques. Such an approach allows us to study star–planet magnetic interactions from a dynamical and evolutive point of view and to compare the associated results to the Strugarek et al. (2017) study, which relied on the instantaneous migration timescale of the planet. For the sake of simplicity, we use a reference case to investigate the influence of each parameter of our model on the secular evolution of star–planet systems. We consider a young star–planet system formed by a fast-rotating K star orbited by a strongly magnetized hot Neptune whose main features are presented in Table 1.

3.1.2 Secular evolution of a reference star–planet system and influence of initial semi-major axis

We summarize the secular evolution of our reference system in Fig. 4. The top panel shows the evolution of the semi-major axis of the planetary orbit and the bottom panel the tidal and magnetic torques applied to the planet (see the thick curves in the figure). Our reference model thus shows an outward migration of the planet after the disk dissipation (gray area on the left), followed by an inward migration after t ~ 350 Myr. This change occurs when the co-rotation radius $r_{\text{corot}}={(G{M}_{\star }/{\Omega }_{\text{c}}^2)}^{1/3}$ (for which the orbital period is equal to the rotation period of the stellar envelope; see the black dashed line in Fig. 4) crosses the orbital distance. Indeed, the co-rotation radius varies throughout the life of the system in a similar way to the stellar rotation rate (see Fig. 3): during the PMS, while the star is contracting, the induced spin-up leads to a decrease in the co-rotation radius; and after the ZAMS, as the stellar structure has stabilized, stellar wind spins the star down, leading to an increase in r_corot. The limit of excitation of the dynamical tide, defined as $P_{\text{orb}}=\frac{1}{2}{P}_{\text{rot}}$, evolves in the same way (see the black dotted line in Fig. 4).

The initial outward migration in our reference model can be attributed to the dynamical tide. As the star–planet system reaches higher semi-major axis (near the ZAMS), the dissipation of inertial waves becomes less and less effective (see the red curves in the bottom panel). The secular evolution of the system is then driven by magnetic torques (in blue), which leads to a more efficient inward migration after 100 Myr compared to a system evolving through tidal effects only (in red in the top panel of Fig. 4), no longer evolving significantly at older ages. Interestingly, the tidal torque drops by two orders of magnitude when the planet then crosses the dynamical tide excitation limit (see bottom panel of Fig. 4) as the equilibrium tide then remains the only contributor. However, because the magnetic torques already dominate during that phase, this does not affect the overall evolution of the system.

We finally also track the limit of formation of a planetary magnetosphere, corresponding to Λ_p = 1, as a limit orbital distance a_cav below which, in our model, a magnetic cavity is formed around the planet. Within our model hypotheses (see Sect. 2.6), this limit can be expressed as $$a_{\text{cav}}=R_{\star }{\left(\frac{B_{\star }}{B_{\text{p}}}\right)}^{\frac{1}{2}}.$$(22)

As we consider a constant planetary magnetic field, a_cav evolves in the same way as the magnetic field at the stellar surface. In the case of our reference system, the orbital semi-major axis is always beyond the cavity formation limit (green curve), which means that the planet is able to sustain a magnetosphere throughout the whole ESPEM simulation. This will be generally the case in what follows, and we highlight the special cases when a magnetic cavity is formed and influences the secular evolution.

We now focus on the influence of the initial semi-major axis on the secular evolution of the system. To this end, we vary the initial semi-major axis of our reference case by considering a_ini = 0.019, 0.025,  0.035, and 0.045 AU. For the three highest values of a_ini, only outward migration occurs initially as they orbit outside the co-rotation radius. Remote planets migrate less efficiently because both tidal and magnetic torques decrease with higher semi-major axis (from dark to light colors in Fig. 4). However, an increase in the initial semi-major axis by a factor of 1.8 leads to a decrease in the magnetic torque by a factor of 2.5 and a drop in tidal torque by at most a factor of 30 when the dynamical tide dominates. Hence, for planets that are able to sustain a magnetosphere, the tidal torque presents a higher sensitivity to the orbital distance than the magnetic torque, which means that for remote planets the secular evolution will likely be dominated by magnetic torques for a higher fraction of the system’s lifetime. In the case of a_ini = 0.019 AU, the planet is initially located below the co-rotation radius and beyond the tidal excitation limit. The planet therefore migrates inward efficiently because of the rise in dynamical tide until it is engulfed very early on, after about 5 Myr (which explains why evolutionary tracks are so short in that case). In this case, the dynamical tide is so efficient that the addition of magnetic torques does not change the already fast evolution of the star–planet system.

In what follows we assess the sensitivity of the secular evolution of a given system to the free parameters of the ESPEM model, namely a_ini,  P_rot,ini,  M_⋆,  M_p, and B_p. To this end, we focus on characterizing the time at which the co-rotation radius exceeds the semi-major axis of the orbit. This gives us as a first idea of the sensitivity of our model to the initial conditions and physical prescriptions we chose. For instance, we have seen that the crossing of the co-rotation radius can be delayed by hundreds of millions of years when the initial semi-major axis varies from 0.025 AU to 0.045 AU. We also found that the addition of magnetic torque to the tidal torques induces a delay of the order of 10 Myr in our reference case. Let us now characterize the sensitivity of our model to stellar (Sect. 3.2) and planetary (Sect. 3.3) parameters.

Fig. 4 Secular evolution of a star–planet system formed by a fast-rotating K star (M⋆ = 0.8 M⊙,  Prot,ini = 1.4 days) orbited by a strongly magnetized hot Neptune (Mp = 0.1 MJup,  Bp = 10 G) for four different initial semi-major axes: aini = 0.019 AU, 0.025 AU,  0.035 AU,  0.045 AU (dark to light colors). The thick curves correspond to our reference case discussed in Sect. 3.1.1 with aini = 0.035 AU. Top panel: semi-major axis (solid lines), co-rotation radius of the star (black dashed lines), and limit of excitation of the dynamical tide (black dotted lines). Wind braking + tides are shown in shades of red; wind braking + tides + magnetic effects are shown in shades of blue. The gray bands on the left correspond to the disk-locking phase. The cavity formation limit, corresponding to Λp = 1, is shown in green. Bottom panel: tidal (shades of red) and magnetic (shades of blue) torques in the case of an evolution with all the combined interactions. The red circles correspond to the crossing of the dynamical tide excitation limit by the planet.

3.2 Influence of stellar parameters on planet migration

3.2.1 Influence of initial stellar rotation and stellar mass

We now investigate the influence of stellar rotation and stellar mass on the evolution of a star–planet system. The relative contribution of magnetic and tidal torques depending on the instantaneous stellar rotation is presented in Appendix B. To highlight the role of initial stellar rotation on the fate of the system, we consider models rotating initially slower than our reference case, that is P_rot,ini = 2.67 and 5 days.

The first striking effect of the initial stellar rotation period is that the planet is generally always closer to its host along the secular evolution for slower initial rotators. Indeed, our reference model shows a first phase of outward migration, followed by an inward migration after the crossing of the co-rotation radius. If a star rotates slowly initially, the first outward migration phase is very inefficient because both the tidal torque and the stellar magnetic field (and thus the magnetic torque) are small. On the contrary, the late inward migration phase is as efficient in all cases as stars converge on the same rotational tracks on the main sequence and therefore the magnetic torques that dominate the evolution here are of comparable amplitude. We note that planetary migration is negligible in the tidal case alone (in red in the upper panel) for the two slowest rotations. Indeed, as the dynamical tide is raised in this configuration, higher stellar rotation periods result in lower valuesof the tidal torque. However, in both cases, the addition of the magnetic torque affects the secular evolution. We also considered the particular case where the planet is initially situated exactly at the co-rotation orbit (P_rot,ini = 2.67 days); it weakly migrates outwards after the dissipation of the disk (see the gray shaded area in Fig. 5) as the star contracts and spins up, before the system evolves in the same way as our reference case, albeit with less efficient planetary migration. The tidal torque is indeed an order of magnitude lower, allowing the magnetic torque to dominate at all times. By taking magnetic interactions into account, a typical variation of P_rot,ini (from 1.4 to 5 days) leads to a delay in the crossing of the co-rotation radius of around 100 Myr.

To investigate the influence of stellar mass on tidal and magnetic torques, we now consider our reference case for three different stellar masses: M_⋆ = 0.5,  0.8, and 1 M_⊙. As shown inthe top panel of Fig. 6, the co-rotation radius, the tidal excitation limit and the semi-major axis have an overall similar evolution as in Sect. 3.1.1 due to the initial conditions adopted. The planet, initially beyond the co-rotation radius, undergoes an outward migration through the dynamical tide in all cases. Near the ZAMS, the secular evolution of the system is driven by magnetic torques in all cases as well (in blue in the top panel of Fig. 6).

In the tidal case alone (in red in the top panel of Fig. 6), the planet undergoes a more efficient migration around more massive stars at the beginning of the evolution. After dozens of millions of years, migration becomes negligible. Indeed, during the PMS, higher mass stars undergo a stronger tidal torque (in red in the bottom panel). Then, as the stellar structure has stabilized around the ZAMS, stellar spin-down leads to a continuous decrease in tidal dissipation towards the end of the evolution, varying weakly with stellar mass. In the presence of magnetic interactions, planetary migration is found to be more and more efficient as stars are less massive. Indeed, those have a higher relative convective mass, a smaller Rossby number (Eqs. (10)–(11)), and thus a stronger stellar magnetic field. This tends to enhance the stellar wind flow as well as star–planet magnetic interactions for low-mass stars. The magnetic torque then significantly affects the later evolution of the semi-major axis and tends to dominate the tidal torque over a longer phase (in blue in the top panel of Fig. 6). More precisely, when magnetic interactions overcome their tidal counterparts (at t ~ 30 Myr for M_⋆ = {0.8, 1} M_⊙, and immediately after the dissipation of the disk for M_⋆ = 0.5 M_⊙), the three evolutions behave similarly until the crossing of the co-rotation radius. Indeed, the magnetic torque depends weakly on stellar mass during the PMS and the beginning of the MS. As less massive stars have a lower co-rotation radius, the planet undergoes an outward migration during a longer phase. This allows it to reach larger semi-major axes for low values of M_⋆.

A typical change in stellar mass (from 0.5 to 1 M_⊙) induces a delay of around 600 Myr of the crossing of the co-rotation radius, which makes it the most sensitive parameter of our model. For low stellar masses (in particular the case M_⋆ = 0.5 M_⊙ in the top panel of Fig. 6), magnetic torques can lead to a migration delay of 100 Myr when compared to the case with only tidal effects.

In summary, the magnetic torque tends to dominate the evolution of the star–planet system during a longer fraction of its lifetime when the stellar magnetic field is stronger and when the dynamical tide is less efficient. This typically occurs for lower mass stars, as was already pointed out by Strugarek et al. (2017), and confirmed here with a fully dynamical evolution.

Fig. 5 Secular evolution of a star–planet system formed by a K star (M⋆ = 0.8 M⊙) orbited by a strongly magnetized hot Neptune (aini = 0.035 AU,  Mp = 0.1 MJup,  Bp = 10 G) for three different initial stellar rotation periods: Prot, ini = 1.4 days,  2.67 days,  5 days (dark to light colors). The thick curves correspond to our reference case discussed in Sect. 3.1.1. Top panel: semi-major axis (solid lines), co-rotation radius of the star (black dashed lines). Wind braking + tides are shown in shades of red; wind braking + tides + magnetic effects are shown in shades of blue. The gray bands on the left correspond to the disk locking phase. Bottom panel: tidal (shades of red) and magnetic (shades of blue) torques in the case of evolution with all the combined interactions. The white squares correspond to the ZAMS and the red circles to the crossing of the dynamical tide excitation limit by the planet.

Fig. 6 Secular evolution of a star–planet system formed by a fast-rotating star (Prot,ini = 1.4 days) orbited by a strongly magnetized hot Neptune (Mp = 0.1 MJup,  aini = 0.035 AU,  Bp = 10 G) for three different stellar masses: M⋆ = {0.5,  0.8,  1} M⊙ (dark to light colors). The thick curves correspond to our reference case discussed in Sect. 3.1.1. Top panel: semi-major axis (solid lines) and co-rotation radius of the star (black dashed lines). Wind braking + tides are shown in shades of red; wind braking + tides + magnetic effects are shown in shades of blue. The gray bands on the left correspond to the disk locking phase. Bottom panel: tidal (shades of red) and magnetic (shades of blue) torques in the case of evolution with all the combined interactions. The white squares correspond to the ZAMS and the red circles to the crossing of the dynamical tide excitation limit by the planet.

3.2.2 Influence of stellar magnetism on planet migration

We now aim to assess the influence of stellar magnetism on the secular evolution of star–planet systems. To this end, we first assess the dependency of Γ_mag on the stellar magnetic field heuristically. Indeed, the magnetic torque, as expressed in Sect. 2.6, presents the following dependencies for low alfvenic Mach numbers: $${\Gamma }_{\text{mag}}\propto \{\begin{array}{l}{M}_{\text{a}}^{0.44}{p}_{\text{tot}}^{0.72}a,\text{ if }{\Lambda }_{\text{p}}>1\hfill \\ M_{\text{a}}{p}_{\text{tot}}a,\text{otherwise (magnetic cavity)}.\hfill \end{array}$$(23)

By introducing the wind magnetic field at the planetary orbit, B_wind, the alfvénic Mach number M_a = |Ω_c − n|a∕v_A, where v_A is the Alfvén velocity, and scales as $$M_{\text{a}}\propto B_{\text{wind}}^{-1}.$$(24)

Moreover, as we consider close-in exoplanets, the total wind pressure is dominated by the magnetic component close to the star (Preusse et al. 2005; Réville et. al. 2015b), which leads to $$p_{\text{tot}}\propto B_{\text{wind}}^2.$$(25)

This results in the following dependency for the magnetic torque: $${\Gamma }_{\text{mag}}\propto B_{\text{wind}}.$$(26)

We now consider a multipolar topology of degree l for the magnetic field (B_wind ∝ a^−(l+2), as the star is at the center of both the wind and the planetary orbit). The magnetic torque then becomes $${\Gamma }_{\text{mag}}\propto B_{\star }{a}^{-(l+2)}\mathcal{F}(a),$$(27)

where $\mathcal{F}$ is a function of the semi-major axis, independent of the degree l at first order, and which is linked to wind acceleration. Even if those scaling laws are based on strong assumptions on the wind model, we can see that a more complex magnetic topology (corresponding to higher l values) implies a stronger dependency of the magnetic torque on the semi-major axis, which can make it more sensitive than the tidal torque itself. If the stellar magnetic field is dominated by small scales (large l) and its large-scale components are weak, the magnetic torque then also weakens efficiently and no longer affects the evolution of remote planets.

In addition, a change in the scaling law of B_⋆ (Eq. (19)) can affect the relative importance of the magnetic torque in our model. To illustrate this we consider the alternative prescription proposed by Ahuir et al. (2020), which shows the steepest Rossby number dependency: $$B_{\star }[\text{G}]=2.0{\left(\frac{Ro}{R{o}_{\odot }}\right)}^{-1.65}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{-1.04},Ro>R{o}_{\text{sat}}.$$(28)

To keep a consistent wind model, the T_c and n_c prescriptions need to be updated as (Ahuir et al. 2020): $$T_{\text{c}}[\text{MK}]=1.5{\left(\frac{Ro}{R{o}_{\odot }}\right)}^{-0.04}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{0.05},Ro>R{o}_{\text{sat}},$$(29) $$n_{\text{c}}[{\text{cm}}^{-3}]=7.25\times {10}^7{\left(\frac{Ro}{R{o}_{\odot }}\right)}^{-0.64}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{1.49},Ro>R{o}_{\text{sat}}.$$(30)

A steeper dependency on the Rossby number implies higher values of B_⋆ at young ages, and thus stronger star–planet magnetic interactions at the beginning of the evolution of the system compared to a scenario in which B_⋆ and Ro⁻¹ scale linearly (Eqs. (19), (20), and (21)). Moreover, as $Ro\propto M_{\star }{R}_{\star }^{1.2}$ during the MS in the Sadeghi Ardestani et al. (2017) formulation, the stellar magnetic field in the present scenario is also more sensitive to M_⋆. Because of the solar normalization, a change in the prescription of B_⋆ then increases the stellar magnetic field of less massive stars, and emphasizes the significance of the magnetic torque on the secular evolution of the star–planet system.

The influence of the choice of a magnetic scenario on the evolution of our reference case is illustrated in Fig. 7. Planetary migration through Eq. (28) (upper panel) is more efficient (light blue curve) than for our reference case (dark blue curve) until t ~ 700 Myr. Afterwards, both evolutions are similar. Indeed, as can be seen from the scaling laws, Eq. (28) induces a more intense magnetic torque for systems younger than the Sun. Hence, at the beginning of the evolution, the magnetic torque is increased by nearly an order of magnitude (in blue in the bottom panel of Fig. 7). The crossing of the co-rotation radius is delayed by about 10 Myr, showing a mild influence of the stellar magnetism scaling law on the system here. As the planet migrates further away from the star than in the case of a linear B_⋆ − Ro⁻¹ scaling, the tidal torque becomes less efficient during the MS (in red in the bottom panel of Fig. 7). However, this does not affect the planetary migration, as the magnetic torque already dominates the evolution in these phases.

Fig. 7 Secular evolution of a star–planet system formed by a fast-rotating K star (M⋆ = 0.8 M⊙,  Prot,ini = 1.4 days) orbited by a strongly magnetized planet (aini = 0.035 AU,  Bp = 10 G) in the case ofa linear B⋆ − Ro−1 relationship (in dark colors; see Eqs. (19), (20), and (21)) and a superlinear B⋆ − Ro−1.65 relationship (in light colors; see Eqs. (28), (29) and (30)). The thick curves correspond to our reference case discussed in Sect. 3.1.1. Top panel: semi-major axis (solid lines) and co-rotation radius of the star (black dashed lines). The gray bands on the left correspond to the disk-locking phase. Bottom panel: tidal (dark and light red) and magnetic (shades of blue) torques in the case of evolution with all the combined interactions. The red circles correspond to the crossing of the dynamical tide excitation limit by the planet.

3.3 Influence of planetary mass and magnetic field on planet migration

We now aimto unravel the influence of planetary intrinsic properties on the fate of the system. We focus in this work on the influence of the mass of the planet as well as its magnetic field. The assessment of the former has a direct influence on the mass dependency of the different torques. Indeed, from the Murray & Dermott (1999) and the Strugarek (2016) formulations we have $${\Gamma }_{\text{tide}}\propto M_{\text{p}}^2$$(31) $${\Gamma }_{\text{mag}}\propto R_{\text{p}}^2{B}_{\text{p}}^{0.56}.$$(32)

In view of the relative dependence on M_p and R_p of the two torques, the mass–radius relationship used to describe the planet will have a significant influence on the relative contributions of magnetic and tidal interactions.

We now consider different scenarios to estimate the magnetic field of the planet. First, if we assume a planetary magnetic field independent of all the other quantities of our model, by relying on Eq. (12) the magnetic torque scales as $${\Gamma }_{\text{mag}}\propto \{\begin{array}{ll}{M}_{\text{p}}^{0.56},\hfill & M_{\text{p}}<2.0M_{\oplus }(6.29\times {10}^{-3}{M}_{\text{Jup}})\hfill \\ M_{\text{p}}^{1.18},\hfill & 2.0M_{\oplus }\le M_{\text{p}}<0.4M_{\text{Jup}}\hfill \\ M_{\text{p}}^{-0.08},\hfill & M_{\text{p}}\ge 0.4M_{\text{Jup}\text{.}}\hfill \end{array}$$(33)

The tidal torque (Eq. (31)) in the scenario we consider (Eq. (33)) is therefore more sensitive to the planetary mass than the magnetic torque.

To investigate these different sensitivities to planetary mass in more detail we consider a lower and a higher planetary mass with respect to our reference case, that is: M_p = 0.01 and 1 M_Jup. As shown inFig. 8, in the tidal case alone, the migration of massive planets is more efficient (red curves in upper panel) than for their less massive counterparts, as the tidal torque is stronger. In the presence of magnetic interactions (in blue), the evolution of the super-Earths differs significantly from the tidal case alone (upper panel). Indeed, the tidal torque decreases more substantially than the magnetic torque for a decreasing planetary mass (see the bottom panel of Fig. 8). Hence, in the case of super-Earths, the magnetic torque (see the blue curves in Fig. 8) is more likely to drive the evolution of the star–planet system. By taking magnetic interactions into account, a typical variation of M_p (from 10⁻² to 1 M_Jup) leads to a delay of around 100 Myr in the crossing of the co-rotation radius. Thus, a change of two orders of magnitude in planetary mass has a similar effect on secular evolution, manifesting as an increase in the initial semi-major axis by a factor of two and an increase in the initial stellar rotation period by a factor of five.

One can assess the robustness of the previous results by considering other hypotheses for the planetary magnetic field. For instance, one can assume that the magnetic field of the planet we consider behaves in the same way as what is observed in the Solar System. Then, according to Shkolnik & Llama (2017), in the case of close-in planets with a synchronous rotation and a dipolar magnetosphere, we obtain $$B_{\text{p}}\propto {\left(\frac{M_{\text{p}}}{P_{\text{orb}}}\right)}^{1.21}{R}_{\text{p}}^{-3}.$$(34)

In such a scenario, according to the Chen & Kipping (2017) M_p − R_p relationships (Eq. (12)), B_p increases for an increasing planetary mass if M_p < 2.0 M_⊕ and M_p ≥ 0.4 M_Jup. This leads to the following mass dependencies for the magnetic torque: $${\Gamma }_{\text{mag}}\propto \{\begin{array}{ll}{M}_{\text{p}}^{0.77},\hfill & M_{\text{p}}<2.0M_{\oplus }(6.29\times {10}^{-3}{M}_{\text{Jup}})\hfill \\ M_{\text{p}}^{0.87},\hfill & 2.0M_{\oplus }\le M_{\text{p}}<0.4M_{\text{Jup}}\hfill \\ M_{\text{p}}^{0.66},\hfill & M_{\text{p}}\ge 0.4M_{\text{Jup}}.\hfill \end{array}$$(35)

This leads to a stronger mass dependency of the magnetic torque for jovian planets as well as super-Earths and a weaker M_p -dependency in the case of neptunian planets. However, as in the constant surface magnetic field scenario, the tidal torque is more sensitive to planetary mass than the SPMI torque. Hence, we find this effect to be negligible and choose to consider a constant planetary field in the remainder of this work for the sake of simplicity.

Let us now assess the sensitivity of our results to the amplitude of the planetary field. In the case of a magnetized planet, as ${\Gamma }_{\text{mag}}\propto R_{\text{p}}^2{B}_{\text{p}}^{0.56}$, a higher planetary magnetic field leads to more efficient star–planet interactions and therefore to a more significant migration. Such a trend is illustrated in Fig. 9 by varying the surface magnetic field of the planet in our reference case. In this configuration, a factor of ten in the magnetic field leads to an increase in the dipolar torque by a factor of 3.6, which has a significant influence on the semi-major axis of the planet during the main sequence. Such a variation plays a critical role if the planet is likely to be engulfed. However, a change of one order of magnitude in the planetary magnetic field only leads to a delay in exceeding the co-rotation radius of the order of 10 Myr, making B_p the least sensitive free parameter of our model.

When a magnetic cavity is formed around the planet, as the geometrical cross-section intervenes in the magnetic torque, the latter is independent of B_p. The magnetic torque then scales as ${\Gamma }_{\text{mag}}\propto R_{\text{p}}^2$. From the Chen & Kipping (2017) relations, the magnetic torque is therefore less sensitive to the planetary mass than its tidal counterpart in the magnetic cavity regime, but is always negligible compared to the dynamical tide because it decreases by one order of magnitude.

Even though planetary magnetism can affect the secular evolution of the star–planet system significantly, the magnetic field strength of extrasolar planets remains poorly constrained. For instance, tidal effects may heat the planetary core, hence affecting its dynamics. Those processes may alter the planetary dynamo in a manner that is not yet fully understood. However, we can gain initial insight into the possible values in various ways. For instance, McIntyre et al. (2019) estimated the magnetic moment of rocky planets from dynamo models, which leads to a surface magnetic field of between 1.5 × 10⁻² and 1.45 G, assuming a dipolar topology. We can also rely on what is observed in the Solar System to assess B_p. While the maximal field strength in this system is observed in Jupiter with a value of around 4 G, the Shkolnik & Llama (2017) scaling law applied to a super-Jupiter of mass M_p= 10 M_Jup and an orbital period of 0.5 days (situated near the Roche limit of a 1 M_⊙ star for example) leads to a planetary magnetic field that can reach 10 G. Hence, we adopt in the rest of this work two different values for the planetary magnetic field: B_p = 1 G, and B_p = 10 G, the latter acting as an upper bound compared to the values measured in the Solar System. Such values of B_p seem to be in agreement with the possible detections of exoplanet radio emission (e.g., τ Boötis b, we refer the reader to Turner et al. 2021). However, one has to keep in mind that it is possible for some planets to have an even stronger magnetic field. Indeed, values as high as 28 G or even hundreds of Gauss have been estimated for hot Jupiters using observed abnormal stellar activity correlated to the planet orbital motion (Cauley et al. 2015, 2019).

Fig. 8 Secular evolution of a star–planet system formed by a fast rotating K star (M⋆ = 0.8 M⊙,  Prot,ini = 1.4 days) orbited by a strongly magnetized planet (aini = 0.035 AU,  Bp = 10 G) for three different planetary masses: Mp = {0.01,  0.1,  1} MJup (dark to light colors). The thick curves correspond to our reference case discussed in Sect. 3.1.1. Top panel: semi-major axis (solid lines), co-rotation radius of the star (black dashed lines). Wind braking + tides are shown in shades of red; wind braking + tides + magnetic effects are shown in shades of blue. The gray bands on the left correspond to the disk-locking phase. Bottom panel: tidal (shades of red) and magnetic (shades of blue) torques in the case of evolution with all the combined interactions. The white squares correspond to the ZAMS and the red circles to the crossing of the dynamical tide excitation limit by the planet.

Fig. 9 Secular evolution of a star–planet system formed by a fast-rotating K star (M⋆ = 0.8 M⊙,  Prot,ini = 1.4 days) orbited by a magnetized hot Neptune (Mp = 0.1 MJup,  aini = 0.035 AU) for two different values of the planetary magnetic field: Bp = {1, 10} G (in dashed and solid line, respectively). The thick curves correspond to our reference case discussed in Sect. 3.1.1. Top panel: semi-major axis (solid lines) and co-rotation radius of the star (black dashed lines). Wind braking + tides are shown in shades of red; wind braking + tides + magnetic effects are shown in shades of blue. The gray bands on the left correspond to the disk-locking phase. Bottom panel: tidal (shades of red) and magnetic (shades of blue) torques in the case of evolution with all the combined interactions. The white square corresponds to the ZAMS and the red circles to the crossing of the dynamical tide excitation limit by the planet.

3.4 Impact on stellar rotation

We now investigate the impact of the planet on the rotation rate of its host star. Indeed, the transfer of angular momentum between the planet and the star can affect stellar rotation, especially in the case of a planet spiralling inwards, which may spin up the host star by transferring its orbital angular momentum into stellar spin (Yee et al. 2020).

We can assess the upper bound of the angular momentum transferred from the orbit to the star by considering the case of a planet engulfment. In this configuration, the initial orbital angular momentum L_orb,ini is entirely transferred to the stellar envelope over the migration timescale. Hence, a planet-hosting star with a rotation period P_rot and an angular momentum L_⋆ has a modified rotational state compared to an isolated star with a rotation period of P_{rot, al} and an angular momentum L_⋆,al. In this context, the difference in angular momenta δL = L_⋆ − L_⋆,al can be written as $$\frac{\delta L}{L_{\star ,al}}=\frac{P_{\text{rot,al}}}{P_{\text{rot}}}-1=\frac{L_{\text{orb,ini}}}{L_{\star ,al}}.$$(36)

The corresponding difference in rotation periods δP = P_rot − P_rot,al then becomes $$\frac{\delta P}{P_{\text{rot,al}}}=\frac{-L_{\text{orb,ini}}/L_{\star ,al}}{1+L_{\text{orb,ini}}/L_{\star ,al}}~-\frac{M_{\text{p}}\sqrt{G{M}_{\star }{a}_{\text{ini}}}{P}_{\text{rot,al}}}{2\pi I_{\star }},$$(37)

where I_⋆ is the total moment of inertia of the star. Hence, stellar rotation is the most impacted by the star–planet interaction when the engulfed planet is massive or when it has a high initial semi-major axis. Furthermore, all things equal, a later engulfment leads to a larger relative over-rotation of the star because older stars rotate slower. As an example, a 1 M_⊙ star engulfing a planet at t = 1 Gyr (corresponding to P_rot,al = 9.33 days) undergoes a three-times-smaller spin-up than if it had destroyed the same planet at solar age (for which P_rot,al = 28 days). More details on this process and on the influence of the different physical quantities can be found in Gallet et al. (2018) and Benbakoura et al. (2019). Here we simply illustrate the influence of the magnetic torque on stellar rotation by varying the planetary mass.

To this end, we consider a star–planet system formed by a median rotating solar twin (M_⋆ = 1 M_⊙, P_rot,ini = 5 days, a_corot,ini = 0.057 AU) orbited by a magnetized planet (a_ini = 0.02 AU,  B_p = 1 G) for three different planetary masses: M_p = {0.01,  0.1,  1} M_Jup. As visible inFig. 10, the planet can have a strong impact on stellar rotation. Indeed, as the planet migrates inwards, angular momentum is transferred from the orbit to the star, resulting in an increase of the stellar spin until planetary engulfment. More detailson the modeling of such a phenomenon in ESPEM can be found in Benbakoura et al. (2019). After the destruction of the planet, the ${\Omega }_{\star }^3$ dependency of wind braking makes the star spin down efficiently, hence causing the two rotational histories (single star and host star engulfing a planet) to converge again. Depending on its lifetime, the star may not come exactly back to an unperturbed rotational state. Furthermore, for low-mass planets, their migration is less efficient, which results in a later engulfment, which in turn leads to a later and smaller peak in stellar rotation, because of the lower initial orbital angular momentum. The magnetic effects accentuate the behavior obtained based on tidal effects only: here the planet experiences an earlier engulfment compared to the pure tidal case (for M_p= {1,  0.1} M_Jup; see the dark blue curves in Fig. 10). The destruction of these planets leads to a spin-up of their host star of around ΔP_rot = 8.9 days and 22.4 days (for M_p = {1,  0.1} M_Jup respectively). These values are slightly lower than the over-rotations obtained through tidal effects only, as planetary engulfment happens earlier. It is also worth noting that a stellar spin-up of around 20% lasts between 4 and 6 Gyr in this configuration. Thus, by relying on gyrochronology, stellar age could be underestimated by about 10%. For M_p= 0.01 M_Jup, tidal effects alone are not strong enough to disrupt the planet. The addition of the magnetic torque allows an engulfment at the end of the evolution (see the light blue curves in Fig. 10), leading to a stellar spin-up of around ΔP_rot = 18.5 days. Hence, star–planet magnetic interactions may actually lead to a late destruction of hot super-Earths.

Fig. 10 Secular evolution of a star–planet system formed by a median rotating solar twin (M⋆ = 1 M⊙,  Prot,ini = 5 days) orbited by a magnetized planet (aini = 0.02 AU,  Bp = 1 G) for three different planetary masses: Mp = {0.01,  0.1,  1} MJup (dark to light colors). Top panel: semi-major axis (solid lines). In shades of red: wind braking + tides, in shades of blue: wind braking + tides + magnetic effects, in shades of purple: Roche limit, defined as in Benbakoura et al. (2019). In dashed-dotted black: stellar radius. Middle panel: time evolution of the rotation rate of the stellar envelope (solid line) and of the stellar core (dashed line). Bottom panel: departure of the rotation period of the planet-hosting star compared to an isolated star.

Fig. 11 Evolution between 1 Myr and 10 Gyr of a planetary population designed from eight orbital periods Porb,ini, equally spaced in logarithm between 0.6 and 10 days, and four initial stellar rotation periods Prot,ini between 1 and 10 days. Those initial parameters correspond to the white circles. Panel A: sketch representing the evolution of the Steady population (S, in yellow), the Young Migrators (Y, in dark orange), and the Old Migrators (O, in light orange) in the (Porb, Prot) plane. The black arrows correspond to the possible paths of each population. Panel B: evolution of the planetary population for M⋆ = 1 M⊙ and Mp = 0.01 MJup. Panel C: evolution of the planetary population for M⋆ = 1 M⊙ and Mp = 1 MJup. Panel D: evolution of the planetary population for M⋆ = 0.6 M⊙ and Mp = 0.01 MJup. Panel E: evolution of the planetary population for M⋆ = 0.6 M⊙ and Mp = 1 MJup. Magnetic and tidal dominance are shown in blue and red, respectively. The shades of blue and red correspond to the overall migration timescale of the system, displayed in logarithmic scale (see key on the right). Regions where the overall migration timescale is greater than the age of the Universe are shown in gray. The black dashed line corresponds to co-rotation and the black dotted line to the tidal excitation limit.

4 Star–planet interactions and classification of planetary populations

4.1 ESPEM sample and migration timescales

Knowing the influence of each parameter of our model on the fate of a given star–planet system, we now adopt a more global approach by simultaneously varying them. To this end we design a sample including 7000 ESPEM simulations which browses the parameter space. The range for each free parameter we consider as well as their distribution are presented in Table 2. To highlight the behavior of star–planet systems through magnetic and tidal interactions, we first focus on a subsample comprising eight initial orbital periods P_orb,ini and four initial stellar rotation periods P_rot,ini, which are evenly spaced in logarithm, two planetary masses M_p = {0.01,  1} M_Jup, and two stellar masses M_⋆ = {0.6,  1} M_⊙.

We show in Fig. 11 the evolutive track of our representative subsample of star–planet models in the (P_orb,  P_rot) plane. The efficiency of planetary migration at each time-step is assessed through the characteristic timescale $${\tau }_{\text{mig}}=\left|\frac{P_{\text{orb}}}{{\dot{P}}_{\text{orb}}}\right|=\frac{1}{3}\left|\frac{L_{\text{orb}}}{{\dot{L}}_{\text{orb}}}\right|.$$(38)

As P_orb ∝ a^3∕2 and L_orb ∝ a^1∕2, for a given migration, the orbital period is more impacted than the semi-major axis, which is itself more sensitive than the orbital angular momentum. Hence, our expression of τ_mig gives the shortest timescale characterizing planet migration. The overall migration timescale due to the sum of tidal and magnetic torques can be split into two contributions $${\tau }_{\text{mig}}={\left(\frac{1}{{\tau }_M}+\frac{1}{{\tau }_T}\right)}^{-1},$$(39)

where τ_M and τ_T are the magnetic and tidal migration timescales respectively. From the expressions of Γ_mag and Γ_tide, those timescales can be written as $${\tau }_T[\text{Myr}]=89.31\times \left(\frac{Q^{\prime }}{{10}^6}\right){\left(\frac{M_{\star }}{M_{\odot }}\right)}^{\frac{8}{3}}\left(\frac{M_{\text{Jup}}}{M_{\text{p}}}\right){\left(\frac{R_{\odot }}{R_{\star }}\right)}^5{\left(\frac{P_{\text{orb}}}{1\text{d}}\right)}^{\frac{13}{3}}$$(40) $$\begin{array}{l}{\tau }_M={\tau }_{M,0}{\left(\frac{R_{\text{J}}}{R_{\text{p}}}\right)}^2\left(\frac{M_{\text{p}}}{M_{\text{Jup}}}\right){\left(\frac{M_{\star }}{M_{\odot }}\right)}^{\frac{5}{3}}{\left(\frac{R_{\star }}{R_{\odot }}\right)}^{-4}{\left(\frac{B_{\star }}{10\text{ G}}\right)}^{-2}{\left(\frac{P_{\text{orb}}}{1\text{d}}\right)}^{\frac{7}{3}}\text{}\text{},\hfill \end{array}$$(41)

with $${\tau }_{M,0}[\text{Myr}]=\{\begin{array}{l}405.3\times {\left(\frac{10R_{\text{p}}}{R_{\text{m}}}\right)}^{1.68}{M}_{\text{a}}^{-0.44},R_{\text{m}}\ge R_{\text{p}}\hfill \\ 2.095\times {10}^5,\text{ otherwise},\hfill \end{array}$$(42)

where $R_{\text{m}}=R_{\text{p}}{\Lambda }_P^{\frac{1}{6}}$ is the size of the planet magnetosphere if we assume a dipolar magnetic field. As we only consider subaflvenic interactions between the star and the planet, the M_a dependency only increases the migration timescale through magnetic torques. Therefore, one can neglect the influence of this quantity to provide a lower bound for τ_M.

The instantaneous overall migration timescale is shown for each system in Fig. 11. As already presented in Sect. 3 for each parameter of the model independently, we can see from Fig. 11 that tidal effects tend to dominate star–planet magnetic interactions (shown in red) for high stellar and planetary masses, long rotation periods (beyond several dozen days in the case of the equilibrium tide) or short rotation periods (below two approximately days in the case of the dynamical tide), and small semi-major axis.

4.2 Classification of planetary populations

Overall, we find that three populations arise from the action of star–planet interactions. The Steady population is composed of planets for which the mean migration timescale is greater than the age of the Universe along most of their evolution because of negligible tidal and magnetic torques. This population is situated at long orbital periods (from 1 to 10 days depending on stellar rotation period; see gray areas in Fig. 11). The Young Migrators are exoplanets that experience a significant migration during the PMS of their host stars. Situated at short orbital and rotation periods (less than two days approximately), these latter are subject to extreme star–planet interactions. Two fates are possible for those systems: If the planet is initially below the co-rotation radius, it is engulfed by the star during the early stages of its evolution. Conversely, if the planet is initially beyond the co-rotation radius, it migrates outward efficiently and may join the Steady population during the MS. However, if magnetic and tidal interactions are efficient enough later on (in the case of low stellar and planetary masses; see panel D in Fig. 11), these can still undergo a significant migration during the MS of their host star, thus becoming Old Migrators (our third migration population; see Fig. 11). Situated at high rotation periods and low orbital periods, the latter population presents the highest sensitivity to the physical parameters of our model, such as the planetary mass and the stellar mass. Their location in the (P_orb, P_rot) plane is that of the most efficient angular momentum transfer from the planet to the star, as is made more explicit below.

The evolution of the three aforementioned populations allows us to define different regions in the (P_orb, P_rot) plane. First, a depleted area is visible at low orbital and rotation periods (Teitler & Königl 2014; Benbakoura et al. 2019). Such a region becomes larger for high stellar and planetary masses, up to orbital periods of 6 days and rotation periods of 2 days in our sample (see e.g., panel C). The depopulation of this region is due to the rapid engulfment of Young Migrators within the co-rotation radius, as well as the efficient outward migration of individuals from the same population initially beyond r_corot. As the dynamical tide dominates the evolution of the systems in this region, the efficiency of this interaction has a direct influence on the extension of the depleted area.

Furthermore, Old Migrators are able to efficiently transfer their orbital angular momentum to their host star. As seen in panel C of Fig. 11, such planets, spiraling inwards, spin up the star and lead to a break in the gyrochronology (e.g., Gallet et al. 2018; Benbakoura et al. 2019), the highest rotation periods remaining unreached (Gallet & Delorme 2019). Such a deserted area at high rotation periods and low orbital periods, only visible for the most massive planets, is extended for high stellar masses (panel C). Then, star–planet systems with an orbital period of less than 3 days and a stellar rotation period of more than 20 days cannot appear, because of the spin-up of the host star. Such a process is driven by magnetic interactions (in blue) for the least massive stars of our sample and by the equilibrium tide (in red) for their massive counterparts.

In addition, the region in the (P_orb,  P_rot) plane populated by Young and Old Migrators defines an area of influence of star–planet interactions. Such a region is the most extended in the case of less massive planets orbiting K-type stars (panel D). There, planets with orbital periods up to 10 days may undergo some migration through the magnetic torque. Indeed, the enhancement of stellar magnetism and stellar wind leads to strong star–planet magnetic interactions. In addition, for giant planets orbiting K-type stars (panel E) and super-Earths orbiting G-type stars (panel B), the range in orbital periods of star–planet interactions is entirely defined by the magnetic torque, as it favors the migration of more distant planets around slower rotators (cf. Sects. 3.1.2 and 3.2.1).

Therefore, a larger number of planets from our sample undergoes slow migration, as the associated timescale is of the order of 1 Gyr. Furthermore, in the case of super-Earths orbiting slowly rotating G-type stars (panel B), the area of influence of star–planet interaction is dictated by the equilibrium tide, as the closest planets are engulfed by the star. If the planet and the star are both massive (panel C), the most distant planets undergoing star–planet interactions are Young Migrators, moving away from the star under the action of the dynamical tide. The area of influence strongly depends on the initial stellar rotation in this case. It can extend up to an orbital period of 2 days for P_rot = 10 days, and beyond 10 days for P_rot = 0.2 days. Magnetic effects then expand the area of influence of star–planet interactions to slower rotators, for which the dynamical tide is less efficient. It is therefore necessary to consider magnetic interactions to determine if a planet is likely to undergo a migration, because in most cases they affect the area of influence of star–planet interactions.

4.3 Influence on the global distributions

The existence of different planetary populations with diverse evolution paths is likely to influence the global distributions of stellar rotation periods and orbital periods in our ESPEM sample. In order to highlight the role played by the different physical parameters of our model, we represent in Fig. 12 the global distribution in P_orb for two different planetary masses (see panel A in Fig. 12) and three different planetary magnetic fields (see panel B in Fig. 12). These distributions are obtained by considering 2000 evolutionary states between the dissipation of the disk and the terminal-age main sequence (TAMS) of our solar-mass models and five initial stellar rotations.

Planetary mass may affect the distribution in orbital periods of the sample significantly. Indeed, as M_p increases, a steeper P_orb distribution is observed (in red in the panel A of Fig. 12), as planetary migration becomes more efficient with higher M_p (see Sect. 3.3). Therefore, less close-in planets are likely to be detected due to higher engulfment rate. A similar effect is observed by increasing the planetary magnetic field (panel B). Indeed, this directly increases the intensity of the magnetic torque, which reinforces the migration of close-in planets. As seen in panel B of Fig. 12, if M_⋆ = 1 M_⊙ and M_p= 0.01 M_J, an increase of one order of magnitude of B_p induces a stronger depopulation only at low orbital periods. However, under more favorable conditions, a stronger planetary magnetic field may favor the repopulation of the same region by more distant planets, as is the case for the lowest stellar masses (see panel D in Fig. 11).

Star–planet interactions may also influence the distribution of P_rot in our sample. Following the same analysis, Fig. 13 shows the percentage of stars with modified rotation as a function of stellar rotation period. In accordance with the results of the previous section, the transfer of angular momentum from the planetary orbit to the star is favored for giant planets orbiting the most massive stars. Thus, the most significant influence of star–planet interactions on the P_rot distribution is due to the spin-up of the more massive stars, which depopulates rotation periods longer than 20 days in favor of P_rot values between 8 and 20 days (in red in Fig. 13), corresponding to stellar ages ranging between 80 Myr and 2 Gyr. However, such effects affect at most 0.5% of the (1 M_⊙,  1 M_Jup) ESPEM subsample, which represents 0.07% of the whole sample. Indeed, such a configuration requires theengulfment of massive planets, thus reducing its probability of occurrence. For M_p= 0.6 M_⊙, only 0.1% of the corresponding subsample (i.e., 0.01% of the total ESPEM sample) undergoes modified rotation for P_rot ≥ 30 days (corresponding to systems older than 5 Gyr). We can therefore consider that star–planet interactions essentially impact the distribution in orbital periods, the distribution in stellar rotation periods being marginally affected here.

Fig. 12 Panel A: global distribution in orbital periods for Mp = 0.01 MJup in black, andMp = 1 MJup) in red.The planetary magnetic field is fixed at 10 G and the stellar mass at 1 M⊙. Panel B: global distribution in orbital periods for Bp = 0 G (in red), Bp = 1 G (in blue), and Bp = 10 G (in black). A star of mass 1 M⊙ and a planet of mass 0.01 MJup are considered here.

Fig. 13 Relative difference in the Prot distribution due to the presence of a planet, for (M⋆, Mp) = (1 M⊙,  0.01 MJup) in black, (M⋆, Mp) = (1 M⊙,  1 MJup) in red, (M⋆, Mp) = (0.6 M⊙,  0.01 MJup) in blue, and (M⋆, Mp) = (0.6 M⊙,  1 MJup) in violet.

Fig. 14 Observed distribution of rotation periods from MMA14 as a function of the effective temperature. The black circle correspond to the solar values.

5 Building a synthetic population of exoplanets

5.1 Observational data

Knowing the possible evolutions of planetary populations as well as the influence of our model physical parameters on the distributions of orbital and stellar rotation periods, we now aim to confront our results with the statistics of observed star–planet systems. To do so, we focus on data from the Kepler mission, and more specifically on two studies performed by McQuillan et al. (2013, hereafter MMA13) and McQuillan et al. (2014, hereafter MMA14).

The MMA14 study provides the largest homogeneous rotation dataset in the Kepler field to date, involving F-type to M-type stars; it includes 34 030 Kepler MS stars whose rotation period has been measured through an autocorrelation-based method, which represents 25.6% of the 133 030 MS Kepler targets detected at that time (excluding known eclipsing binaries and Kepler objects of interest). In this sample, only effective temperatures lower than 6500 Kare considered in order to keep only solar-type stars with convective envelopes. We reproduce the P_rot distribution of the MMA14 targets as a function of their effective temperature in Fig. 14. The MMA14 sample can be seen to present a concentration of stars with rotation periods ranging between 10 and 40 days and effective temperatures between 4000 and 6000 K. Thus, slow-rotating F-type and G-type stars are the most represented in this sample. Furthermore, the distribution of rotation periods presents an upper envelope that increases with decreasing temperature. The Sun is situated close to this upper envelope, which means that few stars older than the Sun have been detected. Indeed, more slowly rotating stars in the MS would not have been removed by exclusion processes of evolved stars. This dataset is used in the remainder of this study to design a synthetic stellar population taking into account the observational biases of the Kepler field as well as the potential influences of the stellar distribution in the galaxy.

MMA13 is one of very few studies to date that combine the orbital period of detected exoplanets and the stellar rotation of their host star, making the associated dataset a valuable tool in the study of star–planet interactions. The main highlight of the study was the observed lack of close-in planets around fast rotators (we also refer the reader to Pont 2009). The MMA13 dataset includes 737KOI with detected orbital and stellar rotation periods. The targets are all in the MS, and were selected using the same processes as in the MMA14 sample. Thus, comparing these two studies allows us to keep identical target selection processes as well as detection methods between the populations of isolated and planet-hosting stars. Moreover, for the sake of consistency, as our model does not deal with interactions between several planets, we do not take into account multiplanetary systems in the MMA13 dataset. Therefore, in the remainder of this work we aim to compare ESPEM results with the MMA13 sample which has been filtered to exclude detected multiplanetary systems. As seen in the top panel of Fig. 15, the distribution in stellar effective temperatures is marginally affected by the choice of dataset (MMA13 vs. MMA14) as well as the removal of multiplanetary systems (black line). Therefore, we can rely on the T_eff distribution from the MMA14 study to design a stellar population and compare the ESPEM results to the filtered MMA13 sample without adding significant biases. In addition, as shown in the bottom panel of Fig. 15, the initial MMA13 sample as well as its filtered analog mostly contain planets with radii between 0.6 and 8 R_⊕, corresponding to planetary masses ranging between 0.2 and 72 M_⊕ according to Chen & Kipping (2017) conversion laws (cf. Eq. (12)). Thus, low-mass planets are extensively represented in these samples.

5.2 Synthetic populations and Kepler observations: comparison of global distributions

In the remainder of this work, we aim to generate a synthetic population of star–planet systems out of our whole sample of ESPEM simulations, and then to compare the distributions we obtain with the filtered MMA13 dataset. To do so, we first design a realistic population of stars based on the MMA14 study. For the sake of consistency, as we assume a two-layer structure for the star, we consider stellar masses ranging between 0.5 and 1.1 M_⊙ with a bin size of 0.1 M_⊙. This corresponds to effective temperatures of between 3700 and 6000 K, which allows us to cover most of the MMA14 sample. The associated surface gravity logg ranges between 4.25 and 5.07, which is consistent with the selection processes used in MMA13 and MMA14 to exclude likely giants (Ciardi et al. 2011). Furthermore, in order to account for the rotational evolution of stars in open clusters (Gallet & Bouvier 2015), in our synthetic populations generated with ESPEM we consider five initial rotation periods evenly distributed between 1 and 10 days. We then rely on a common sampling of stellar age for all evolutionary tracks in order to remove biases due to the adaptative time-step used in the ESPEM integrator (for more details, we refer the reader to Benbakoura et al. 2019) as well as to account for the relative duration of the different stages of stellar evolution for each stellar mass. More precisely, stellar age is sampled between 1 Myr and 10 Gyr to obtain 2000 equally spaced instants. Then, for a given star, only stellar ages that fall in the MS are taken into account. To make sure that, in the absence of planets, our stellar sample reproduces the MMA14 stellar distribution, we introduce a coefficient C_stellar(P_rot, T_eff) to weight a given star–planet system at a given age according to the rotation period of the host star as well as its effective temperature(for more details, we refer the reader to Appendix A). This allows us to indirectly bias the stellar masses and ages of the systems according to the observations of the Kepler mission.

We then include planets by considering 40 initial semi-major axes ranging between 5 × 10⁻³ and 0.2 AU, evenly spaced in logarithm. The corresponding orbital periods are then situated between 0.12 and 46 days, which allows us to initially populate all the orbital distances for which star–planet interactions may act efficiently. Five planetary masses between 0.5 Earth masses and 5 Jupiter masses, uniformly spaced in logarithm, are considered in this planetary population. This corresponds to planetary radii of between 0.8 and 12.7 R_⊕, thus covering the majority of the planets considered by the MMA13 study, as seen in Fig. 15. Each planetary mass is weighted in our population in order to reproduce the R_p distribution of the filtered MMA13 subsample. Moreover, three configurations are investigated for a given population of planetary systems: we make it evolve through tidal effects only, or along with magnetic effects by considering B_p = 1 G and 10 G. As the initial configuration of the systems and the engulfment ratio are unknown in the MMA13 study, we only take into account the planets that survived. This allows us to define a density of probability of presence DPP_ESPEM(P_orb, P_rot). The probability $\text{d}{\mathbb{P}}_{\text{ESPEM}}$ of finding a star–planet system with an orbital period included within the interval [P_orb, P_orb + dP_orb] and a stellar rotation period in [P_rot, P_rot + dP_rot] is then defined as $$\text{d}{\mathbb{P}}_{\text{ESPEM}}(P_{\text{orb}},P_{\text{rot}})={\text{DPP}}_{\text{ESPEM}}(P_{\text{orb}},P_{\text{rot}})\text{d}{P}_{\text{orb}}\text{d}{P}_{\text{rot}}.$$(43)

To compare these results with the observed distributions, we assume that each individual from the ESPEM sample corresponds to a detected star–planet system. This way, the probability ${\mathbb{P}}_{\text{ESPEM}}$ can be interpreted as a distribution of detected systems assuming that only star–planet interactions affect the initial population. However, several biases may be involved in the MMA13 distributions (such as the probability of transit, which is equal to R_⋆ ∕a, and also biases linked to the selection of KOIs as well as the detectability of stellar rotation period). Thus, synthetic distributions may show a larger number of systems with increasing orbital period when compared with the corresponding observed planet population. Hence, to compare consistently synthetic and observed distributions, we first consider a subsample of the filtered MMA13 sample, for which orbital periods are shorter than 20 days. This way, we only focus on close-in planets, for which our study is relevant. Each system is then weighted by the inverse of the probability of transit a∕R_⋆ to remove the associated bias. In the remainder of Sect. 5.2, both synthetic and observed distributions are normalized to 1. Finally, we systematically quantify the differences between the planetary global distributions by performing a two-sample Kolmogorov-Smirnov (KS, Kolmogorov 1933; Smirnov 1948) test. The corresponding statistics of the test and p-values are presented in Table C.1.

The probability of finding a planet-hosting star at a rotation period P_rot, regardless of the position of the planet, is presented in panel A of Fig. 16. By construction, the synthetic populations produced by ESPEM for isolated stars reproduce the MMA14 distribution (in violet in panel A of Fig. 16). In the presence of a planet, the ESPEM distributions show fewer fast rotators (for which P_rot < 10 days) and more stars whose rotation period ranges between 10 and 30 days. Hence, they tend to reproduce the MMA13 distribution (shown in gray in panel A of Fig. 16). Moreover, taking into account magnetic torques for a planetary magnetic field B_p = 1,  10 G marginally affects the P_rot distributions. This implies that only tidal effects are likely to play a role in the distribution of stellar rotation periods. Furthermore, no significant difference in the KS statistic and the corresponding p-value is observed between the synthetic populations (we refer the reader to Appendix C for more details). Such a discrepancy between the P_rot distribution of isolated stars and planet-hosting stars cannot be explained by the alteration of stellar rotation through star–planet interactions,whether the planet is detected or not. Indeed, as seen in Sect. 4.3, planetary migration marginally affects the distribution in stellar rotation periods; the most favorable case of a massive star and a massive planet emphasizes the lowest periods of rotation (see the red curve in Fig. 13). We can therefore conclude that the engulfment of planets orbiting fastrotators through tidal effects makes the detection of these stars less likely, favoring slower rotators in the P_rot distribution of MMA13.

The distributions in orbital periods, regardless of stellar rotation, are shown in panel B of Fig. 16. Taking into account star–planet magnetic interactions, in particular for B_p = 10 G (see the dark blue curve in panel B), significantly affects the P_orb distributions by depopulating orbital periods shorter than 2 days. Such a modification then makes it possible to approach the MMA13 distribution for P_orb < 1 days (see the gray histogram in panel B) with a confidence level of 1σ (we refer the reader to Appendix A for more details about the calculation of the error bars), where the populations for B_p = 0 and 1 G show a more significant excess of planets at these orbital periods compared to observations (red and light blue curves in panel B). More precisely, taking into account star–planet magnetic interactions with a stronger planetary magnetic field tends to improve the goodness of fit by lowering the value of the test statistic from 0.32 to 0.27 and increasing the corresponding p-value from 0.17 to 0.34 (we refer the reader to Table C.1 for more details). For orbital periods of greater than 10 days, the ESPEM distributions become uniform, as star–planet interactions are negligible for remote planets. In any case, taking the magnetic torque for large values of B_p into accountis necessary to account for the distribution of close-in planets. Such a scenario seems to be supported by planetary dynamo considerations. Indeed, as planetary rotation is synchronized in our model, close-in planets may rotate with periods shorter than1 day. They are therefore likely to generate high magnetic fields. Moreover, Cauley et al. (2019) inferred values of B_p of between 20 G and 120 G for planets whose orbital period ranges between 2 and 4 days. Such a configuration would lead to even stronger modifications of the P_orb distribution than what is presented in this work. Furthermore, as we consider a sufficiently wide range of initial semi-major axes, corresponding to orbital periods of between 0.12 and 46 days, the initial conditions adopted to generate our synthetic populations have no influence on the planetary distributions at low orbital periods. Despite a better goodness of fit at high planetary magnetic fields, we see an excess of exoplanets for all the synthetic populations at low orbital periods compared to the observed distribution.This may highlight that other star–planet interactions not considered in this work, or the initial distribution of orbital periods resulting from planetary formation, may have played a role in shaping the observed population from MMA13.

Fig. 15 Top: distribution in effective temperature from MMA14 (in red), MMA13 (in blue), and a subsample of MMA13 excluding multiplanetary systems (in black). Bottom: distribution in planetary radii from MMA13 (in blue) and a subsample of MMA13 excluding multiplanetary systems (in black).

Fig. 16 Distribution of stellar rotation periods (left) and orbital periods (right). The Prot distribution from the MMA14 study is shown in violet. The gray bars correspond to the distributions obtained with the unbiased MMA13 sample. The ESPEM distributions with Bp = 0, 1, 10 G are shownin red, light blue and dark blue respectively. The calculation of the error bars for the synthetic distributions is presented in Appendix A.

Fig. 17 Observed distribution of stellar rotation periods from the whole MMA13 sample as a function of the orbital period of the planet. The purple dashed line corresponds to the lower edge of the distribution fitted in MMA13. The black dashed line corresponds to the synchronization between Porb and Prot.

Fig. 18 Distribution of orbital periods for super-Earths (left column) and giant planets (right column). Each row corresponds to a population of young (Region 1, for which Prot <4.7 days), middle-aged (Region 2, for which 4.7 ≤ Prot [days] < 20), and old (Region 3, Prot ≥ 20 days) star–planet systems. The gray histogram corresponds to the distributions obtained with the unbiased MMA13 sample. The ESPEM distributions with Bp = 0, 1, 10 G are shown in red, light blue and dark blue respectively. The gray, red, and purple horizontal bands at the top of each panel correspond to the range of values in each region of the co-rotation, the dynamical tide excitation limit, and the lower edge of the MMA13 distribution, respectively. The calculation of the error bars for the synthetic distributions is presented inAppendix A.

5.3 Synthetic populations and Kepler observations: comparison of young, middle-aged, and old star–planet systems

We now turn to the detailed distributions of star–planet systems at different ages. To do so, we rely on the (P_orb, P_rot) distribution of the whole MMA13 sample shown in Fig. 17, which allows us to define three regions. For stars rotating with a period shorter than around 4.7 days (the so-called Region 1 in Fig. 17), few planets are detected. Moreover, in Region 2, corresponding to stellar rotation periods ranging between 4.7 and 20 days, the detected exoplanets have an orbital period greater than around 1 day. Finally, a wide range of orbital periods, namely between 0.4 and 500 days, is observed for exoplanets orbiting around the slowest rotators (P_rot > 20 days; see the Region 3 in Fig. 17).

For each of these regions, each corresponding to a row in Fig. 18, we compare the distributions of the synthetic populations to their MMA13 counterparts in the case of super-Earths (M_p < 10 M_⊕, see the first column of Fig. 18) and giant planets (M_p ≥ 10 M_⊕, see the second column of Fig. 18). The observed star–planet systems are selected and unbiased with the same method as in Sect. 5.2. Both synthetic and observed distributions of super-Earths and giant planets for all the stellar rotation periods considered are then normalized to 1.

In the case of fast rotators (panels A and B), Young Migrators have already migrated outwards or been engulfed, depending on their initial position relative to the co-rotation orbit (see the horizontal gray bands at the top of each panel in Fig. 18). The frequency of occurrence then may be higher for orbital periods longer than 3 days for the most massive planets (panel B), and longer than around 1 day in the case of super-Earths, as star–planet interactions are less efficient (panel A). Taking into account star–planet magnetic interactions has a marginal influence on the distribution of giant planets, as tidal interactions dominate secular evolution. However, a slight shift towards higher orbital periods is observed if M_p < 10 M_⊕ (see red and dark blue curves in panel A), with more planets being engulfed. In those regions, the small number of observations and the low probability of occurrence in the case of synthetic populations prevent significant discrimination between the different distributions (see upper rows in Table C.2).

In the case of giant planets located in Region 2 (panel D), no significant migration takes place for orbital periods longer than 3 days. The associated distribution is then almost uniform. When P_orb < 3 days, Old Migrators get closer to the star through the action of the equilibrium tide and to a lesser extent star–planet magnetic interactions. This induces a depopulation of the lowest orbital periods, slightly enhanced by the adjunction of the magnetic torque. More precisely, the case B_p = 10 G leads to slightly better KS statistic compared to the other synthetic populations. In the case of the slowest rotators (Region 3 ; see panel F in Fig. 18), the frequency of occurrence decreases with decreasing orbital period. This can be attributed to planetary engulfment, which is at the origin of the deserted area presented in Sect. 4.2. Moreover, as the extension of this region decreases with stellar mass (see Sect. 4.2; we also refer the reader to Gallet & Delorme 2019), the frequency of the occurrence of massive planets decreases smoothly with their orbital period. Such a behavior results from the action of tidal effects, with magnetic torque playing a negligible role. Moreover, we see that the probability of detecting a massive planet around slow rotators in the MMA13 sample is generally lower than what can be expected from synthetic populations. As giant planets are more likely to be detected than super-Earths, such a discrepancy might not be entirely attributed to the low number of systems detected in this region. Therefore, other migration mechanisms may be responsible for the engulfment of giant planets, thus reducing their frequency of occurrence around slowly rotating stars. For instance, the excitation of tidal gravity waves in the stellar radiative zone may affect the evolution of the system (Barker & Ogilvie 2010; Guillot et al. 2014; Ahuir et al. 2021). Such a contributionwill be studied in future work. Furthermore, magnetic fields of strength greater than 10 G (Cauley et al. 2015, 2019) may also significantly affect the distributions of the giant planets by depopulating the shortest orbital periods.

Star–planet magnetic interactions significantly affect P_orb distributions in the case of low-mass planets orbiting slower rotators (panels C and E). Indeed, Older Migrators efficiently get closer to thestar at the beginning of the MS thanks to the magnetic torque. The closest planets are then likely to be engulfed through the combined action of the equilibrium tide and magnetic torque. In this configuration, we recover the characteristics of the global distribution in orbital periods. Therefore, taking into account a strong magnetic torque (for B_p = 10 G; see the darkblue curve on panels C and E of Fig. 18) significantly modifies the distributions by depopulating the orbital periods shorter than two days, and these distributions then show a better agreement with the observations of MMA13 at low orbital periods with a confidence level of 1σ. Indeed, the cases B_p = 0 G and B_p = 1 G present an excess of planets at periods of between 1 and 2 days. More precisely, taking into account star–planet magnetic interactions with a planetary magnetic field of 10 G significantly improves the KS statistic and the corresponding p-value compared to the other synthetic populations (we refer the reader to Table C.2 for more details). However, as in the global P_orb distribution, an excess of planets at orbital periods shorter than 3 days is still noticeable.

Thus, it is possible to approach the P_orb − P_rot distribution observed in Kepler systems (e.g., from the MMA13 study) by relying on a realistic stellar population (such as the MMA14 distribution) and star–planet interactions after the dissipation of the protoplanetary disk. More precisely, the magnetic torque tends to modify the distribution of super-Earths around slower rotators, which improves the agreement between synthetic populations and observations.

6 Summary and discussions

6.1 Results

In this paper, we focus on the secular evolution of star–planet systems by taking stellar magnetic braking and star–planet magnetic and tidal interactions into account simultaneously. The two latter effects act together on planetary migration and stellar rotation. Furthermore, both interactions may dominate the other throughout secular evolution depending on the initial configuration of the system and the evolutionary phase considered. More precisely, tidal effects are found to dominate star–planet magnetic interactions for high stellar and planetary masses as well as low semi-major axis. This implies in particular that super-Earths close to their host star essentially migrate through magnetic torques. Those interactions may actually lead to a late destruction of hot super-Earths. The dynamical tide governs the evolution of planets orbiting fast rotators while slower rotators evolve through magnetic interactions. However, at high rotation periods, the equilibrium tide may become the dominating contribution during planetary engulfment. Moreover, stellar rotation may be significantly impacted at any age due to the engulfment of a close-in planet. In this configuration, an alteration of the stellar rotation period of several dozens of days may last more than a few million years. However, such events are scarce, as they require high stellar and planetary masses, or sufficiently high planetary magnetic fields in the case of super-Earths, to be truly effective.

Three populations of star–planet systems emerge from the combined action of magnetic and tidal torques. Systems undergoing negligible migration define an area of influence of star–planet interactions, which can extend for instance up to orbital periods of 10 days for super-Earths orbiting K-type stars. Furthermore, for sufficiently large planetary magnetic fields, the magnetic torque determines the extension of this region. Planets outside this influence region form our first population, which we dub a “Steady” population.

The second population we identified, which we call the “Young Migrators”, is composed of planets initially close to fast rotators that migrate efficiently during the PMS. They may either be expelled away from the star or be rapidly engulfed, which engenders a depleted region at low rotation and orbital periods.

Finally, we identified a third population of “Old Migrators”, composed of planets migrating inward around slower rotators, which happens during the MS. This population is more sensitive to the physical parameters involved in our modeling. They can also lead to an efficient angular momentum transfer from the planetary orbit to stellar rotation, like Young Migrators. They could induce a break of gyrochronology for high stellar and planetary masses, as the star can efficiently spin up. This population finally creates a region at high stellar rotation periods and low orbital periods not populated by star–planet systems.

Furthermore, star–planet interactions significantly impact the global distribution in orbital periods. Indeed, for higher planetary mass and planetary magnetic fields, magnetic and tidal torques are more efficient. Low orbital periods are then more likely to be depopulated. However the global distribution in stellar rotation periods is marginally affected. Indeed, around 0.5% of G-type starsand 0.1% of K-type stars may spin up because of planetary engulfment. As a significant stellar over-rotation requires the destruction of massive planets, its probability of occurrence in a given planetary population is found to be relatively low.

Finally, we designed synthetic populations based on observed stellar distributions (such as the MMA14 distribution). We found that star–planet magnetic interactions, after the dissipation of the disk, significantly affect the distribution of super-Earths around slower rotators, which improves the agreement between synthetic populations and observations, while tidal effects shape the distribution of giant planets. Such a result is obtained without relying on specific prescriptions for the initial semi-major axis distribution of planets after the disk dissipation. However, we found that all these populations present an excess of exoplanets at low orbital periods compared to the distribution observed in Kepler systems (e.g., from the MMA13 study). This may indicate that additional star–planet interactions not taken into account in this work, such as the dynamical tide in the stellar radiative zone (e.g., Terquem et al. 1998; Goodman & Dickson 1998; Barker & Ogilvie 2010) or magnetic interactions in the unipolar regime (Laine et al. 2008; Laine & Lin 2012), are at play in shaping the observed close-in planet population in the Kepler field. Another possibility is that part of the observed distribution is actually already present immediately after the disk dissipation, and differs from the initial uniform distribution in P_orb we assume here. A combination of the two aspects could be shaping the observed distribution, and a detailed investigation is left for future work.

6.2 Discussions

The results presented in this work may have strong implications on exoplanet detection. For instance, they imply that the detection of jovian planets with orbital periods shorter than 6 days around solar twins rotating with a period shorter than 2 days is very unlikely because of the presence of a depleted area. Moreover, similar systems having an orbital period of shorter than 3 days and a stellar rotation period of longer than 20 days could hardly be detected because of potential planetary engulfment. Such a drop in the frequency of occurrence is less significant when considering weaker stellar and planetary masses. More generally, planets with orbital periods of up to 10 days orbiting stars with a rotation period shorter than 10 days are most likely to undergo migration through magnetic (and to a lesser extent tidal) interactions. In the case of slower rotators, stellar mass, initial rotation, and planetary mass have a strong impact on the possibility of planetary migration (cf. Sect. 4.2).

Moreover, regarding the influence of star–planet interactions on stellar rotation, with ESPEM we recover the results from previous studies (e.g., Zhang & Penev 2014; Gallet et al. 2018; Benbakoura et al. 2019; Gallet & Delorme 2019). Indeed, as seen in Fig. 10, the engulfment of a jovian planet by a solar twin may lead to an alteration of the stellar rotation period of more than 90%, which corresponds to an error of 45% in the estimation of stellar age if we assume the Skumanich law. It is worth noting that an over-rotation of around 20% is likely to last for a few billion years. A deviation from gyrochronology is also possible for stars hosting super-Earths as star–planet magnetic interactions may lead to their engulfment later on during the main sequence. The induced spin-up is less significant is this case. We show in this work that these effects are barely noticeable in global P_rot distributions, as at most 0.07% of a population of star–planet systems may see its stellar rotation significantly altered. A signature of planetarymigration may nevertheless become apparent at a given age, in particular in stellar open clusters (Teitler & Königl 2014; Gallet et al. 2018), and confronted with alternative scenarios that have been proposed to account for the observed P_rot distribution of low-mass stars (Brown 2014; van Saders et al. 2016).

Theoretical models predict that planetary migration within the protoplanetary disk is more efficient than the effects considered in this work (e.g., Baruteau et al. 2014; Bouvier & Cébron 2015; Heller 2019). This leads to the idea that the planetary population achieved immediately after disk dissipation is close to being steady until later phases of stellar evolution after the main sequence. Nonetheless, we show that the magnetic and tidal interactions occurring after the dissipation of the disk have a significant influence on planetary populations during the PMS and MS. As a result, our work shows that the observed populations need to be carefully studied to derive constraints on planet formations, as both disk (not taken into account here) and post-disk migration come into play for the population on short-period orbits.

In addition, studying planetary populations for orbital periods of less than 2 days is essential for understanding planetary systems andtheir evolution. Indeed, such orbital periods seem to highlight the presence of strong planetary magnetic fields (cf. Sect. 5.3). They may also highlight the presence of some star–planet interactions not taken into account in this work (see Sect. 6.3). Uncertainties remain regarding the value of the planetary magnetic field B_p. Indeed, this latter could reach values much higher than those assumed in this work (e.g., Cauley et al. 2015, 2019). In particular, from dynamo considerations, magnetic fields as high as 4000 G are expected in hot Jupiters during the PMS (Hori 2021). This would significantly affect the planetary distributions by enhancing the depopulation of star–planet systems at short orbitalperiods, whether in the case of super-Earths or even giant planets. Finally, although star–planet interactions may be necessary toaccount for planetary populations, they only concern planets close to their host star, which are thus located below the habitable zone (for an extensive study of the secular evolution of this region, we refer the reader to Gallet et al. 2017b).

6.3 Perspectives

The present work is a first attempt to study the secular evolution of star–planet systems through magnetic and tidal interaction with a fully dynamical approach. We also attempt to investigate their role in shaping the distributions of planetary populations. In order to best account for the observations, we need to take these types of interactions into account in a comprehensive manner. For instance, to get a complete picture of tidal dissipation in stars, the dynamical tide in the stellar radiative zone needs to be taken into account (e.g., Zahn 1975; Goldreich & Nicholson 1989; Goodman & Dickson 1998; Terquem et al. 1998), which is likely to compete with the dissipation of inertial waves in convective layers (Ivanov et al. 2013) and to affect secular evolution of star–planet systems (Barker & Ogilvie 2010; Barker 2011; Guillot et al. 2014; Barker 2020). Moreover, in differentially rotating stellar convective zones, tidal inertial waves may interact with mean flows at critical layers; they can therefore either deposit or extract angular momentum from or to the surrounding fluid, which leads to exchanges of angular momentum between the star and the planet (Astoul et al. 2021). Finally, tides can be affected by stellar and planetary magnetic fields (Wei 2016; Lin & Ogilvie 2018; Astoul et al. 2019).

Regarding star–planet magnetic interactions, the unipolar interaction may occur for a weakly magnetized planet with a low magnetic diffusivity (Laine et al. 2008; Laine & Lin 2012). Taking such a regime into account is likely to increase the magnetic torque by several orders of magnitude and to make the subsequent planetary migration even more efficient (Strugarek et al. 2017). New behaviors could then be expected. As an example, if the planet gets closer to its host star, the latter may significantly spin up, which may enhance the stellar magnetic field, favoring a transition from the dipolar to the unipolar regime, which has been found to be at the origin of a strong increase in magnetic torque. Planetary migration is then likely to enter a runaway regime.

Here we only considered isolated, circularized star–planet systems. More complex geometries, for example including eccentricities, inclinations (Kaula 1961), and dynamical interactions in multi-planet systems (e.g., Laskar et al. 2012; Bolmont et al. 2015), have to be implemented in a future work to fully characterize the PMS and MS evolution of multi-planetary systems.

Acknowledgements

We would like to thank the anonymous referee for helpful comments and suggestions regarding our work. The authors thank Emeline Bolmont for the help and information she has generously provided. We also thank Cillia Damiani and Noë Brucy-Ciaramella for fruitful discussions regarding our work. The authors acknowledge funding from the European Union’s Horizon-2020 research and innovation programme (Grant Agreement no. 776403 ExoplANETS-A) as well as the PLATO CNES funding at CEA/IRFU/DAp. A.S. and A.S.B. acknowledge funding by ERC WHOLESUN 810218 grant, INSU/PNST and CNES Solar Orbiter. A.S. and J.A. acknowledge funding from the Programme National de Planétologie (PNP). S.M. and J.A. acknowledge funding by the European Research Council through the ERC grant SPIRE 647383.

Appendix A Synthesis of a population of star–planet systems

A.1 Stellar population

To confront the MMA13 planetary distribution with the one obtained through the action of star–planet interactions, we rely on a synthetic population of star–planet systems generated with our ESPEM code. The goal of this section is to explain in more detail the methodsto generate this population.

Fig. A.1 Values of the Cstellar coefficient asa function of the stellar rotation period Prot and the effective temperature Teff.

We first need to design a stellar population. To do so, we consider stellar masses ranging between 0.5 and 1.1 M_⊙ with a bin of 0.1 M_⊙. Five initial rotation periods evenly distributed between 1 and 10 days are chosen to account for the rotational evolution of stars in open clusters (Gallet & Bouvier 2015). Stellar age is sampled between 1 Myr and 10 Gyr to obtain 2000 dates equally spaced, regardless of the simulation being processed. Then, for a given star, only the stellar ages falling inside the MS are taken into account. The stellar sample is then generated thanks to ESPEM by making each star evolve without a planet. This way, the stellar rotation period evolves due to changes in stellar structure, redistribution of angular momentum inside the star, and magnetic braking. We then determine the frequency of occurrence f_ESPEM,al(T_eff, P_rot) of a star at a given pair (T_eff, P_rot). By defining the same frequency of occurrence f_MMA14(T_eff, P_rot) for the MMA14 sample, we introduce the coefficient $$C_{\text{stellar}}(T_{\text{eff}},P_{\text{rot}})=\frac{f_{\text{MMA14}}(T_{\text{eff}},P_{\text{rot}})}{f_{\text{ESPEM,al}}(T_{\text{eff}},P_{\text{rot}})}.$$(A.1)

Such a ratio quantifies the difference between the P_rot distribution obtained by ESPEM for isolated stars and the one observed in the MMA14 study. This difference may be due to the limitations of our rotational evolution model, to the distribution of stars in the galaxy, to observational biases, or to uncertainties on the detection of rotation periods in MMA14 (for an assessment on the reliability of the MMA14 results, see Santos et al. 2019).

As shown in Fig. A.1, high values of C_stellar(P_rot, T_eff) indicate an under-representation of the corresponding stars in the ESPEM population compared to the MMA14 distribution. In order to reproduce the latter, it is necessary to give greater weight to stars with effective temperatures between 5000 and 6000 K and rotation periods between 1 and 30 days. This then amounts to favoring the presence of G-type stars younger than the Sun within the synthetic stellar population.

We now aim to assess in a first approach the influence of the uncertainties in effective temperature and stellar rotation period on the C_stellar coefficient. To this end, we consider that each star from the MMA14 sample has an effective temperature T_eff and a stellar rotation period P_rot. We define the associated uncertainties σ[P_rot], given in MMA14, and σ[T_eff]. For the latter quantity, we assume a typical value of 200 K (we refer the reader to Huber et al. 2014, for more details). We assume as a first approximation that each observation from MMA14 follows a 2D normal distribution $\mathcal{N}(T_{\text{eff}},\sigma [T_{\text{eff}}])\times \mathcal{N}(P_{\text{rot}},\sigma [P_{\text{rot}}])$.

We now introduce a bin in effective temperature and stellar rotation period (T_k, P_k), with dimensions (ΔT_k, ΔP_k), and p_ik the probability that the observation i falls into the bin k. Hence we have: $$\begin{array}{ll}{p}_{ik}=\hfill & \left[{\Phi }_{T_{\text{eff}},\sigma [T_{\text{eff}}]}\left(T_k+\frac{\Delta T_k}{2}\right)-{\Phi }_{T_{\text{eff}},\sigma [T_{\text{eff}}]}\left(T_k-\frac{\Delta T_k}{2}\right)\right]\hfill \\ \hfill & \times \left[{\Phi }_{P_{\text{rot}},\sigma [P_{\text{rot}}]}\left(P_k+\frac{\Delta P_k}{2}\right)-{\Phi }_{P_{\text{rot}},\sigma [P_{\text{rot}}]}\left(P_k-\frac{\Delta P_k}{2}\right)\right],\hfill \end{array}$$(A.2)

where $${\Phi }_{\mu ,\sigma }(x)=\frac{1}{2}\left[1+\text{erf}\left(\frac{x-\mu }{\sigma \sqrt{2}}\right)\right]$$(A.3)

is the cumulative distribution function of a normal distribution with a mean μ and a standard deviation σ. The number N_k of observations in the bin k can now be interpreted as a sum of Bernoulli variables B_ik, equal to 1 if the observation i is in the bin k, and 0 otherwise. By assuming that all the observations from the MMA14 sample have been performed independently, the standard deviation associated to N_k becomes $$\sigma [N_k]=\sqrt{{\sum }_i{p}_{ik}(1-p_{ik})}.$$(A.4)

By introducing the number of stars N_MMA14 in the MMA14 sample, we obtain $$\sigma [f_{\text{MMA14}}]=\frac{\sigma [N_k]}{\sqrt{N_{\text{MMA14}}}}.$$(A.5)

Next, from Eq. (A.1), one can assess the associated standard deviation σ[C_stellar] as $$\frac{\sigma [C_{\text{stellar}}]}{C_{\text{stellar}}}=\frac{\sigma [f_{\text{MMA14}}]}{f_{\text{MMA14}}}.$$(A.6)

As shown in Fig. A.2, the relative uncertainty of C_stellar is higher at low effective temperatures and low rotation periods, reaching values as high as 100%. However, it marginally affects synthetic populations,as few systems from the MMA14 sample are located in this region. In the regions favored by the MMA14 sample, σ[C_stellar]∕C_stellar is less than 1%. A more extensive study of those uncertainties is left for future work.

A.2 Star–planet population

Once the stellar parameters have been chosen and consistently biased, we include the planets by choosing 40 initial semi-major axes ranging between 5 × 10⁻³ and 0.2 AU, uniformly spaced in logarithm, and five planetary masses between 0.5 Earth masses and 5 Jupiter masses, uniformly spaced in logarithm. From these initial conditions, we can define three different star–planet populations by making them evolve through tidal effects only, or along with magnetic effects by considering a planet magnetic field equal to 1 and 10 G. Each system is weighted by the C_stellar factor. This way, the stellar distribution is consistent with the MMA14 study. Furthermore, each planetary mass is also weighted to reproduce the MMA13 distribution of planetary radii. We then compute a biased frequency of occurrence $f_{\text{ESPEM}}^{\text{biased}}(P_{\text{orb}},P_{\text{rot}})$, accounting for both the MMA13 and MMA14 distributions: $$f_{\text{ESPEM}}^{\text{biased}}=f_{\text{ESPEM}}(P_{\text{orb}},P_{\text{rot}})C_{\text{stellar}}(T_{\text{eff}},P_{\text{rot}})f_{\text{MMA13}}(M_p).$$(A.7)

The associated distribution function DF_ESPEM is defined as $$\text{d}{n}_{\text{ESPEM}}(P_{\text{orb}},P_{\text{rot}})={\text{DF}}_{\text{ESPEM}}(P_{\text{orb}},P_{\text{rot}})\text{d}{P}_{\text{orb}}\text{d}{P}_{\text{rot}},$$(A.8)

where dn_ESPEM is the fraction of systems with an orbital period included within the interval [P_orb, P_orb + dP_orb] and a stellar rotation period in [P_rot, P_rot + dP_rot]. As the rotation periods as well as the orbital periods are sampled, we calculate this distribution by considering the frequency of occurrence at each bin of surface ΔP_orbΔP_rot as $${\text{DF}}_{\text{ESPEM}}(P_{\text{rot}},P_{\text{orb}})=\frac{f_{\text{ESPEM}}^{\text{biased}}(P_{\text{rot}},P_{\text{orb}})}{\Delta P_{\text{orb}}\Delta P_{\text{rot}}}.$$(A.9)

The function obtained is then normalized such as $${\displaystyle \iint {\text{DF}}_{\text{ESPEM}}}(P_{\text{orb}},P_{\text{rot}})\text{d}{P}_{\text{orb}}\text{d}{P}_{\text{rot}}=1-\mathcal{D},$$(A.10)

with $\mathcal{D}$ the ratio of the planets from the initial population that have been engulfed. Such an approach makes it possible to compare synthetic populations with different values of B_p, as the rate of planet engulfment may change because of the variable efficiency of the magnetic torque (see Table A.1).

Fig. A.2 Values of σ[Cstellar]∕Cstellar as a functionof the stellar rotation period Prot and the effective temperature Teff.

However, in the case of planetary populations observed by the Kepler mission we have no idea of the associated initial population. In particular, the rate of planets engulfed by their host star is unknown. To compare our synthetic populations with the observations, we only take into account the planets in our sample that have survived. To do so, we define a probability density of presence DPP_ESPEM(P_orb, P_rot), obtained by normalizing the distribution function DF_ESPEM(P_orb, P_rot) so that $${\displaystyle \iint {\text{DPP}}_{\text{ESPEM}}}(P_{\text{orb}},P_{\text{rot}})\text{d}{P}_{\text{orb}}\text{d}{P}_{\text{rot}}=1.$$(A.11)

Fig. A.3 Probability of presence ${\mathbb{P}}_{\text{ESPEM}}$ (solid lines) and fraction nESPEM of systems present in the synthetic population (dashed lines) as a function of orbital period. In red: ESPEM distributions with Bp = 0 G. In light blue: ESPEM distributions with Bp = 1 G. In dark blue: ESPEM distributions with Bp = 10 G.

Thus, as shown in Fig. A.3, only the relative shape of the distributions is considered with such a normalization. More precisely, a higher concentration of planets at a given orbital period induces higher probabilitiesof presence, regardless of the number of planets destroyed in the sample. The probability of presence ${\mathbb{P}}_{\text{ESPEM}}$ (solid lines in Fig. A.3) is thus higher than the fraction n_ESPEM (dashed lines in Fig. A.3) of systems present in the sample. This is particularly true in case B_p = 10 G, because the associated planet destruction rate is the highest.

To assess the influence of uncertainties of the MMA14 sample on these distributions, we consider for the sake of simplicity the relative uncertainty of C_stellar averaged in effective temperatures. Furthermore, to decorrelate the contributions from the MMA13 and MMA14 samples, we do not consider any uncertainty in planetary radii, as the same distribution is used in both observed and synthetic populations. Hence, from Eq. (A.7), the standard deviation $\sigma \left[{\text{DPP}}_{\text{ESPEM}}\right]$ linked to the density of probability of presence can be assessed as $$\frac{\sigma [{\text{DPP}}_{\text{ESPEM}}]}{{\text{DPP}}_{\text{ESPEM}}(P_{\text{orb}},P_{\text{rot}})}=\left(\frac{\sigma [C_{\text{stellar}}]}{C_{\text{stellar}}}\right)(P_{\text{rot}}).$$(A.12)

We then rely on this prescription to build the error bars associated to the distributions presented in Sects. 5.2 and 5.3.

Appendix B Influence of instantaneous stellar rotation

To highlight the role of instantaneous stellar rotation on the fate of the system, we consider the ratio between thetidal and magnetic torques as a function of the stellar rotation period for the reference case presented inSect. 3.1 and a similar star–planet system with an initial stellar rotation period P_rot,ini = 9 days and an initial semi-major axis a_ini = 0.025 AU. In both cases, the planet is able to sustain a magnetosphere. As seen in Fig. B.1, when inertial waves are excited in the stellar envelope, the associated torque tends to dominate for lower rotation periods. Thus, the dynamical tide generally dominates the evolution during most of the PMS. If the inertial waves cannot be excited by the tidal potential, only the equilibrium tide contributes to the evolution of the system. In this case the ratio between the tidal and magnetic torques increases for higher rotation periods, as the tidal torque is less sensitive to stellar rotation than the magnetic torque. Hence, at the end of the MS (seethe highest rotation periods in Fig. B.1) the tidal torque can be of the same order of magnitude as the magnetic torque or even dominate the latter, depending on the initial configuration of the system.

Fig. B.1 Ratio between the tidal and magnetic torques as a function of the stellar rotation period for our reference case (solid line) and for a similar star–planet system with Prot,ini = 9 days and aini = 0.025 AU (dotted line). In gray: regions where the overall migration timescale is greater than the age of the universe. Blue background:magnetic dominance. Red background: tidal dominance. The white circles correspond to the beginning of the evolution. The white squares correspond to the ZAMS. The orange circles correspond to solar age.

Appendix C Results of the Kolmogorov–Smirnov test

We perform a two-sample Kolmogorov–Smirnov test to assess whether the observed and synthetic samples show the same underlying probability distribution in orbital and stellar rotation periods, which constitutes the null hypothesis. To perform the KS test, the probability distributions are first normalized to 1. We then compute the KS statistic ${\mathcal{D}}_{\text{KS}}$, defined as $${\mathcal{D}}_{\text{KS}}={\text{sup}}_{P_{\text{orb}}}\left|F_{\text{ESPEM}}(P_{\text{orb}})-F_{\text{MMA13,unbiased}}(P_{\text{orb}})\right|,$$(C.1)

where F corresponds to the cumulative distribution function. We then rely on the Kolmogorov distribution K to estimate the p-value as follows: $$\mathbb{P}\left(K\ge \sqrt{\frac{N}{2}}{\mathcal{D}}_{\text{KS}}\right)=2{\sum }_{k=1}^{+\infty }{(-1)}^{k-1}{e}^{-2k^2{\left(\sqrt{\frac{N}{2}}{\mathcal{D}}_{\text{KS}}\right)}^2},$$(C.2)

where N is the size of both observed and synthetic samples. We present the global cumulative distribution functions in P_orb and P_rot, for both the observed and synthetic samples, in Fig. C.1. The corresponding KS statistics and p-values, mentioned in Sect. 5.2, are presented in Table C.1. We also present the cumulative distribution functions for super-Earths and giant planets, for both the observed and synthetic samples, in Fig. C.2. The corresponding KS statistics and p-values, mentioned in Sect. 5.3, are presented in Table C.2.

Fig. C.1 Cumulative distributions of stellar rotation periods (left) and orbitals periods (right). In gray: unbiased MMA13 sample. In red: ESPEM distributions with Bp = 0 G. In light blue: ESPEM distributions with Bp = 1 G. In dark blue: ESPEM distributions with Bp = 10 G.

Fig. C.2 Cumulative distributions of orbital periods for super-Earths (left column) and giant planets (right column). Each row corresponds to a population of young (Region 1, for which Prot <4.7 days), middle-aged (Region 2, for which 4.7 ≤ Prot (days) < 20), and old (Region 3, Prot ≥ 20 d) star-planet systems. In gray: unbiased MMA13 sample. In red: ESPEM distributions with Bp = 0 G. In light blue: ESPEM distributions with Bp = 1 G. In dark blue: ESPEM distributions with Bp = 10 G.

References

All Tables

All Figures

Fig. 1 Sketch of the main features and locations of interest involved in a star–planet system undergoing tidal and magnetic interactions. The planet (in blue) orbits the star (in orange) with an orbital angular velocity n. As a result of the presence of a planet, the star presents a bulge misaligned with the line joining the centers of the two celestial bodies, which induces a lag angle δ. This angle has been greatly exaggerated here for visualization purposes (it is indeed much smaller than 1 degree in almost all cases). If the planetary mean motion is greater than twice the stellar rotation rate, inertial waves cannot be excited in the stellar convective zone and no dynamical tide is raised within the convective envelope of the star (see the purple area). The relative motion between the planet and the ambient wind (represented with orange arrows) leads to the formation of Alfvén wings (orange lobes around the planet) if the planet is below the Alfvén radius (in blue). Beyond thisdistance, no Alfvén wings can connect the star and the planet (see the gray area). For both tidal and magnetic interactions, a planet situated below the co-rotation radius (see the black dashed line) undergoes an inward migration, and a planet situated beyond this distance migrates outwards. The relative position of the different orbits of interest mayvary depending on the initial configuration of the system considered.

In the text

Fig. 2 Schematic view of the system and its interactions (adapted from Benbakoura et al. 2019). The radiative core (in yellow) and the convective envelope (in orange) exchange angular momentum (green arrows). Stellar wind carries away angular momentum from the envelope and spins the star down (purple arrows). Stellar rotation and planetary orbit are coupled through tidal (red arrows) and magnetic effects (blue arrows).

In the text

Fig. 3 Secular evolution of the rotation rate of the convective envelope (solid lines) and of the radiative core (dashed lines) for stars with a mass M⋆ = {1,  0.8,  0.5} M⊙ (from top to bottom). Slow (orange), median (light blue), and fast (dark blue) rotators are considered. The solar rotation rate at solar age is represented by a black circle in each panel. The dots with error bars correspond to the 25th, 50th, and 90th percentiles of rotational distributions of observed stellar clusters published by Gallet & Bouvier (2015). In black: Skumanich law normalized to the Sun.

In the text

Fig. 4 Secular evolution of a star–planet system formed by a fast-rotating K star (M⋆ = 0.8 M⊙,  Prot,ini = 1.4 days) orbited by a strongly magnetized hot Neptune (Mp = 0.1 MJup,  Bp = 10 G) for four different initial semi-major axes: aini = 0.019 AU, 0.025 AU,  0.035 AU,  0.045 AU (dark to light colors). The thick curves correspond to our reference case discussed in Sect. 3.1.1 with aini = 0.035 AU. Top panel: semi-major axis (solid lines), co-rotation radius of the star (black dashed lines), and limit of excitation of the dynamical tide (black dotted lines). Wind braking + tides are shown in shades of red; wind braking + tides + magnetic effects are shown in shades of blue. The gray bands on the left correspond to the disk-locking phase. The cavity formation limit, corresponding to Λp = 1, is shown in green. Bottom panel: tidal (shades of red) and magnetic (shades of blue) torques in the case of an evolution with all the combined interactions. The red circles correspond to the crossing of the dynamical tide excitation limit by the planet.

In the text

Fig. 5 Secular evolution of a star–planet system formed by a K star (M⋆ = 0.8 M⊙) orbited by a strongly magnetized hot Neptune (aini = 0.035 AU,  Mp = 0.1 MJup,  Bp = 10 G) for three different initial stellar rotation periods: Prot, ini = 1.4 days,  2.67 days,  5 days (dark to light colors). The thick curves correspond to our reference case discussed in Sect. 3.1.1. Top panel: semi-major axis (solid lines), co-rotation radius of the star (black dashed lines). Wind braking + tides are shown in shades of red; wind braking + tides + magnetic effects are shown in shades of blue. The gray bands on the left correspond to the disk locking phase. Bottom panel: tidal (shades of red) and magnetic (shades of blue) torques in the case of evolution with all the combined interactions. The white squares correspond to the ZAMS and the red circles to the crossing of the dynamical tide excitation limit by the planet.

In the text

Fig. 6 Secular evolution of a star–planet system formed by a fast-rotating star (Prot,ini = 1.4 days) orbited by a strongly magnetized hot Neptune (Mp = 0.1 MJup,  aini = 0.035 AU,  Bp = 10 G) for three different stellar masses: M⋆ = {0.5,  0.8,  1} M⊙ (dark to light colors). The thick curves correspond to our reference case discussed in Sect. 3.1.1. Top panel: semi-major axis (solid lines) and co-rotation radius of the star (black dashed lines). Wind braking + tides are shown in shades of red; wind braking + tides + magnetic effects are shown in shades of blue. The gray bands on the left correspond to the disk locking phase. Bottom panel: tidal (shades of red) and magnetic (shades of blue) torques in the case of evolution with all the combined interactions. The white squares correspond to the ZAMS and the red circles to the crossing of the dynamical tide excitation limit by the planet.

In the text

Fig. 7 Secular evolution of a star–planet system formed by a fast-rotating K star (M⋆ = 0.8 M⊙,  Prot,ini = 1.4 days) orbited by a strongly magnetized planet (aini = 0.035 AU,  Bp = 10 G) in the case ofa linear B⋆ − Ro−1 relationship (in dark colors; see Eqs. (19), (20), and (21)) and a superlinear B⋆ − Ro−1.65 relationship (in light colors; see Eqs. (28), (29) and (30)). The thick curves correspond to our reference case discussed in Sect. 3.1.1. Top panel: semi-major axis (solid lines) and co-rotation radius of the star (black dashed lines). The gray bands on the left correspond to the disk-locking phase. Bottom panel: tidal (dark and light red) and magnetic (shades of blue) torques in the case of evolution with all the combined interactions. The red circles correspond to the crossing of the dynamical tide excitation limit by the planet.

In the text

Fig. 8 Secular evolution of a star–planet system formed by a fast rotating K star (M⋆ = 0.8 M⊙,  Prot,ini = 1.4 days) orbited by a strongly magnetized planet (aini = 0.035 AU,  Bp = 10 G) for three different planetary masses: Mp = {0.01,  0.1,  1} MJup (dark to light colors). The thick curves correspond to our reference case discussed in Sect. 3.1.1. Top panel: semi-major axis (solid lines), co-rotation radius of the star (black dashed lines). Wind braking + tides are shown in shades of red; wind braking + tides + magnetic effects are shown in shades of blue. The gray bands on the left correspond to the disk-locking phase. Bottom panel: tidal (shades of red) and magnetic (shades of blue) torques in the case of evolution with all the combined interactions. The white squares correspond to the ZAMS and the red circles to the crossing of the dynamical tide excitation limit by the planet.

In the text

Fig. 9 Secular evolution of a star–planet system formed by a fast-rotating K star (M⋆ = 0.8 M⊙,  Prot,ini = 1.4 days) orbited by a magnetized hot Neptune (Mp = 0.1 MJup,  aini = 0.035 AU) for two different values of the planetary magnetic field: Bp = {1, 10} G (in dashed and solid line, respectively). The thick curves correspond to our reference case discussed in Sect. 3.1.1. Top panel: semi-major axis (solid lines) and co-rotation radius of the star (black dashed lines). Wind braking + tides are shown in shades of red; wind braking + tides + magnetic effects are shown in shades of blue. The gray bands on the left correspond to the disk-locking phase. Bottom panel: tidal (shades of red) and magnetic (shades of blue) torques in the case of evolution with all the combined interactions. The white square corresponds to the ZAMS and the red circles to the crossing of the dynamical tide excitation limit by the planet.

In the text

Fig. 10 Secular evolution of a star–planet system formed by a median rotating solar twin (M⋆ = 1 M⊙,  Prot,ini = 5 days) orbited by a magnetized planet (aini = 0.02 AU,  Bp = 1 G) for three different planetary masses: Mp = {0.01,  0.1,  1} MJup (dark to light colors). Top panel: semi-major axis (solid lines). In shades of red: wind braking + tides, in shades of blue: wind braking + tides + magnetic effects, in shades of purple: Roche limit, defined as in Benbakoura et al. (2019). In dashed-dotted black: stellar radius. Middle panel: time evolution of the rotation rate of the stellar envelope (solid line) and of the stellar core (dashed line). Bottom panel: departure of the rotation period of the planet-hosting star compared to an isolated star.

In the text

Fig. 11 Evolution between 1 Myr and 10 Gyr of a planetary population designed from eight orbital periods Porb,ini, equally spaced in logarithm between 0.6 and 10 days, and four initial stellar rotation periods Prot,ini between 1 and 10 days. Those initial parameters correspond to the white circles. Panel A: sketch representing the evolution of the Steady population (S, in yellow), the Young Migrators (Y, in dark orange), and the Old Migrators (O, in light orange) in the (Porb, Prot) plane. The black arrows correspond to the possible paths of each population. Panel B: evolution of the planetary population for M⋆ = 1 M⊙ and Mp = 0.01 MJup. Panel C: evolution of the planetary population for M⋆ = 1 M⊙ and Mp = 1 MJup. Panel D: evolution of the planetary population for M⋆ = 0.6 M⊙ and Mp = 0.01 MJup. Panel E: evolution of the planetary population for M⋆ = 0.6 M⊙ and Mp = 1 MJup. Magnetic and tidal dominance are shown in blue and red, respectively. The shades of blue and red correspond to the overall migration timescale of the system, displayed in logarithmic scale (see key on the right). Regions where the overall migration timescale is greater than the age of the Universe are shown in gray. The black dashed line corresponds to co-rotation and the black dotted line to the tidal excitation limit.

In the text

	Fig. 12 Panel A: global distribution in orbital periods for Mp = 0.01 MJup in black, andMp = 1 MJup) in red.The planetary magnetic field is fixed at 10 G and the stellar mass at 1 M⊙. Panel B: global distribution in orbital periods for Bp = 0 G (in red), Bp = 1 G (in blue), and Bp = 10 G (in black). A star of mass 1 M⊙ and a planet of mass 0.01 MJup are considered here.
In the text

	Fig. 13 Relative difference in the Prot distribution due to the presence of a planet, for (M⋆, Mp) = (1 M⊙,  0.01 MJup) in black, (M⋆, Mp) = (1 M⊙,  1 MJup) in red, (M⋆, Mp) = (0.6 M⊙,  0.01 MJup) in blue, and (M⋆, Mp) = (0.6 M⊙,  1 MJup) in violet.
In the text

	Fig. 14 Observed distribution of rotation periods from MMA14 as a function of the effective temperature. The black circle correspond to the solar values.
In the text

	Fig. 15 Top: distribution in effective temperature from MMA14 (in red), MMA13 (in blue), and a subsample of MMA13 excluding multiplanetary systems (in black). Bottom: distribution in planetary radii from MMA13 (in blue) and a subsample of MMA13 excluding multiplanetary systems (in black).
In the text

Fig. 16 Distribution of stellar rotation periods (left) and orbital periods (right). The Prot distribution from the MMA14 study is shown in violet. The gray bars correspond to the distributions obtained with the unbiased MMA13 sample. The ESPEM distributions with Bp = 0, 1, 10 G are shownin red, light blue and dark blue respectively. The calculation of the error bars for the synthetic distributions is presented in Appendix A.

In the text

	Fig. 17 Observed distribution of stellar rotation periods from the whole MMA13 sample as a function of the orbital period of the planet. The purple dashed line corresponds to the lower edge of the distribution fitted in MMA13. The black dashed line corresponds to the synchronization between Porb and Prot.
In the text

Fig. 18 Distribution of orbital periods for super-Earths (left column) and giant planets (right column). Each row corresponds to a population of young (Region 1, for which Prot <4.7 days), middle-aged (Region 2, for which 4.7 ≤ Prot [days] < 20), and old (Region 3, Prot ≥ 20 days) star–planet systems. The gray histogram corresponds to the distributions obtained with the unbiased MMA13 sample. The ESPEM distributions with Bp = 0, 1, 10 G are shown in red, light blue and dark blue respectively. The gray, red, and purple horizontal bands at the top of each panel correspond to the range of values in each region of the co-rotation, the dynamical tide excitation limit, and the lower edge of the MMA13 distribution, respectively. The calculation of the error bars for the synthetic distributions is presented inAppendix A.

In the text

	Fig. A.1 Values of the Cstellar coefficient asa function of the stellar rotation period Prot and the effective temperature Teff.
In the text

	Fig. A.2 Values of σ[Cstellar]∕Cstellar as a functionof the stellar rotation period Prot and the effective temperature Teff.
In the text

	Fig. A.3 Probability of presence ${\mathbb{P}}_{\text{ESPEM}}$ (solid lines) and fraction nESPEM of systems present in the synthetic population (dashed lines) as a function of orbital period. In red: ESPEM distributions with Bp = 0 G. In light blue: ESPEM distributions with Bp = 1 G. In dark blue: ESPEM distributions with Bp = 10 G.
In the text

Fig. B.1 Ratio between the tidal and magnetic torques as a function of the stellar rotation period for our reference case (solid line) and for a similar star–planet system with Prot,ini = 9 days and aini = 0.025 AU (dotted line). In gray: regions where the overall migration timescale is greater than the age of the universe. Blue background:magnetic dominance. Red background: tidal dominance. The white circles correspond to the beginning of the evolution. The white squares correspond to the ZAMS. The orange circles correspond to solar age.

In the text

	Fig. C.1 Cumulative distributions of stellar rotation periods (left) and orbitals periods (right). In gray: unbiased MMA13 sample. In red: ESPEM distributions with Bp = 0 G. In light blue: ESPEM distributions with Bp = 1 G. In dark blue: ESPEM distributions with Bp = 10 G.
In the text

Fig. C.2 Cumulative distributions of orbital periods for super-Earths (left column) and giant planets (right column). Each row corresponds to a population of young (Region 1, for which Prot <4.7 days), middle-aged (Region 2, for which 4.7 ≤ Prot (days) < 20), and old (Region 3, Prot ≥ 20 d) star-planet systems. In gray: unbiased MMA13 sample. In red: ESPEM distributions with Bp = 0 G. In light blue: ESPEM distributions with Bp = 1 G. In dark blue: ESPEM distributions with Bp = 10 G.

In the text