Tidal evolution of close-in exoplanets and host stars

ABSTRACT

The evolution of exoplanetary systems with a close-in planet is ruled by the tides mutually raised on the two bodies and by the magnetic braking of the host star. This paper deals with consequences of this evolution and some features that can be observed in the distribution of the systems' two main periods: the orbital period and the stars rotational period. The results of the simulations are compared to plots showing both periods as determined from the light curves of a large number of Kepler objects of interest. These plots show important irregularities as a dearth of systems in some regions and accumulations of hot Jupiters in others. It is shown that the accumulation of short-period hot Jupiters around stars with rotation periods close to 25 d results from the evolution of the systems under the joint action of tides and braking, and requires a relaxation factor for Solar-type stars of around |$10 \, \mathrm{ s}^{-1}$|⁠.

1 INTRODUCTION

One difficulty faced by tidal evolution theories is the extremely slow time-scale of these phenomena and the almost impossibility of observing them at work. So far, tidal infall rates have only been estimated for two exoplanets: WASP-12b (29 ± 2 ms yr⁻¹, cf. Yee et al. 2020) and Kepler-1658b (131 ± 20 ms yr⁻¹, cf. Vissapragada et al. 2022), both obtained from analysing variations in the orbital period of the planet. For all other systems, tidal evolution must be deduced indirectly from present-day dynamical structures that may have been generated or affected by tidal interactions.

Together with magnetic braking, the transfer of angular momentum from the companion orbit to the rotation of the host star plays an important role in defining the rotational and orbital evolution of the system (e.g. Bolmont et al. 2012; Teitler & Königl 2014; Ferraz-Mello et al. 2015). Key elements in this case are the physical parameters of the star and of the planetary companion, and the age of the system – one important datum of difficult determination. We may also use our knowledge of the distribution of some physical parameters to infer some properties of the tidal models. One example is the analysis done by Hansen (2010) using the distribution of eccentricities, periods, and masses of the known planets, which allowed him to constrain the dissipation values of stars and planets.

The top left-hand frame of Fig. 1 is an excerpt of the full figure published by McQuillan, Mazeh & Aigrain (2013) and based on the three first years of public Kepler data. It shows the distribution of stellar rotational periods P_rot = 2π/Ω₀ and planetary orbital periods P_orb = 2π/n; hereafter the P–P diagram. As shown in the inset, the size of each circle is proportional to the planetary radius, while its colour is determined by the effective temperature of the host star. We note an almost absence of systems with fast rotating host stars and T_eff < 5800 K. This last characteristic is just due to the loss of angular momentum of the colder stars due to stellar winds, the so-called magnetic braking (see Bouvier, Forestini & Allain 1997). The hottest stars do not have the outer convective layer responsible for the winds and so do not migrate upwards in this figure as the colder stars.

As indicated by McQuillan et al. (2013), the large majority of close-in planets (P_orb < 5 d) are located above the continuous line drawn in all panels of Fig. 1, while a few also appear to lie close to a synchronous state (P_orb = P_rot) highlighted with a dashed line. The region between both these lines appears surprisingly empty.

The same feature is found when restricting the data to smaller planets, as seen in the middle left-hand frame, where only bodies with R < 5R_⊕ are shown. Several explanations were proposed for the lack of planets below the continuous line, including the combined effects of tidal evolution, braking, and magnetic star–planet interactions (e.g. Teitler & Königl 2014; Ahuir, Mathis & Amard 2021) as well as tidal capture of near-parabolic bodies in multiplanet systems (Lanza & Shkolnik 2014). None of these works, however, focus on the distribution of large planets, seen here in the lower left-hand plot of Fig. 1. Not counting the two Kepler objects of interest (KOIs) close to the synchronous line, most bodies seem to show a significant correlation, with larger stellar rotation associated to larger orbital periods and lower values of P_rot for hot Jupiters closer to the star.

The right-hand frames of Fig. 1 show the same distributions in the P–P diagram, but employing more recent data. Stellar rotations for main sequence M and K stars are taken from Santos et al. (2019) while similar information for G and F stars are found in Santos et al. (2021). These were then compared with the latest catalogue of Kepler candidates, while the same data base also allowed to identify and remove possible eclipsing binaries. The total number of data points increased from 1079 to 1698. Even so, it is important to keep in mind that many correspond to Kepler candidates (KOI) and therefore do not constitute confirmed planets. This explains why some data points in the left-hand column are absent in the latest catalogues.

While some characteristics of the distribution in the P–P diagram remain, others have undergone noticeable changes. The diagonal continuous line now appears less decisive as a lower bound for orbital periods of small close-in planets, and the middle right-hand plot shows several new systems for a wider of range of stellar rotations. In fact, more recent data from both Kepler and TESS (Messias et al. 2022) suggest that the dearth of close-in small planets around rapidly rotating stars could be due to a lack of data and thus not statistically significant.

Conversely, the correlation between P_orb and P_rot in hot Jupiters (lower right-hand plot) appears even more pronounced and better defined, and practically all bodies are located in a moraine-like accumulation. The only exception is KOI 554, found very close to the synchronous line. This system appears as unconfirmed in NASA’s lists and is a probable false positive.

In this paper we present a simple dynamical model to analyse the orbital and rotational evolution of a single planet around a star, under the combined effects of tidal interactions and magnetic braking. We will show that the observed distribution of small planets as well as the accumulation observed for giant planets can be well reproduced, and, at least for Solar-type stars, allow for a reliable estimation for the stellar relaxation factor γ.

2 THE DYNAMICAL MODEL

We computed the evolution tracks of close-in planets taking into account tidal effects and the magnetic braking of the star. Braking of the stellar rotation was computed following Bouvier et al. (1997) while the tidal evolution of both the planet and host star were calculated using the creep tide theory (Ferraz-Mello 2013, 2015). Equations and details of the dynamical model are given in Appendix  A.

The codes adopted for this paper differ slightly from those used in the study of the interplay of tidal evolution and stellar wind braking in the rotation of stars hosting massive close-in companions (Ferraz-Mello et al. 2015). The main improvements include the introduction of the actual fluid Love number k_f so that the equations correspond to the case of a differentiated body whose layers are aligned homogeneous ellipsoidal shells, and the introduction of the effects due to the shortening of the polar axis by the tidal potential (Ferraz-Mello 2015).

In all the simulations, the physical set-up was assumed frozen. Certainly, some parameters involved in the calculations are expected to vary during the whole simulated story. The mass of the star is expected to vary a small amount during the system lifetime. This variation could be included in the simulations, but would not be able to change significantly the results. Large variations may be expected in the planetary masses if the planet is too close to the star because of evaporation of its outer layers. However, the influence of these variations in the simulation results is negligible. The variation in the orbital parameters when masses are not constant is proportional to the derivative of the sum of the two masses (Hadjidemetriou 1963); the planet has nothing but a tiny fraction of the total mass of the system and so, the sum of the masses will be almost unaffected by variations in the mass of the planet even if they are large.

Other parameters showing variation with the evolution of the star are the fluid Love number and the radius of the star. They are related to the density profile of the star. The radius of the star may have a significant variation, but it appears in the equations only through k_fR²; models of the evolution of the internal structure of the stars (Claret 2019) show that the product k₂R² and the moment of inertia have just a small variation during the time in which the stars remain in the main sequence and variations can be neglected. More complex models taking into account the variations in the internal structure of the star would introduce new non-universal unknown parameters without introducing significant changes in the results.

Table 1 summarizes the main physical parameters adopted for each type of planet, including nominal values for the relaxation factor γ. We denote by α the multiplicative factor in the expression of the body’s moment of inertia, i.e. C = αmR². The adopted values of the planetary k_f are based on those calculated for similar Solar system planets by Gavrilov & Zharkov (1977). Different choices for k_f affect the estimations of the relaxation factor γ since these two quantities are entangled in the tidal model. In a first approximation, for gaseous bodies, the tidal variation of the elements is proportional to the ratio k_f/γ (see Appendix  A). However, in the present case, the contribution of the tides on the planet to the variation of the orbital period is insignificant and the planetary γ and k_f are given only to complete the information on the parameters used in the simulations.

Fig. 2 shows the results of some preliminary simulations involving a Jupiter-like planet orbiting a Solar-type star. While the planetary parameters were taken from Table 1, for the central mass we adopted a moment of inertia, |$C_0=0.07 m_0 R_0^2$|⁠, similar to the present-day Sun. Its fluid Love number was chosen equal to |$k_{\mathrm{ f}_0}=0.05$|⁠, as obtained using the equations derived by Folonier, Ferraz-Mello & Kholshevnikov (2015), and assuming a density profile equal to the standard solar model. The stellar relaxation factor was fixed at γ = 10 s⁻¹, in the middle of the range indicated by previous studies (see Ferraz-Mello 2022), but corrected for the smaller value of |$k_{\mathrm{ f}_0}$|⁠.

The left-hand and center plots show, in blue, the time evolution of the planet’s orbital period P_orb for four different initial values: |$2 \, , \, 2.5 \, , \, 3$|⁠, and 4  d. For each case we considered two different initial eccentricities: e = 0 and e = 0.2, while the initial rotation frequency of the planet was taken equal to 2π/Ω₁ = 10 d. This value is arbitrary, but we found no significant difference as long as the body was initially subsynchronous (P_orb < P_rot). Finally, both graphs differ in the initial rotational period of the star, as indicated by the text on top. Fast rotators are considered on the left while initial slow rotators are considered for the center plot. The tidal evolution was followed according to the equations described in the Appendix A, considering both the star and planet as extended bodies and including Cayley functions E_0,k and E_2,k up to order k = 7.

The pink lines show the evolution of the planetary rotational period. The orbital circularization, together with the synchronization of the planetary spin (Ω₁ = n), occurs early in the evolution of the system, as indicated by the superposition of the pink and blue lines. The longest time-scale is found for P_orb(t = 0) = 4 d, where the synchronization occurs after 1 Gyr. However, the subsequent orbital evolution is almost negligible even after T = 10 Gyr. The stellar rotation period, however, fueled by magnetic braking, increases monotonically, reaching values of the order of 30–40 d at the end of the simulation.

Planets closer to the star tell a different story. Synchronization occurs very early in the system’s history and most of the changes in P_orb occur in a scenario dominated by stellar tides and where the planetary counterparts are negligible. Tidal evolution for these initial conditions are much stronger than those for P_orb(t = 0) = 4 d, and the bodies are engulfed by the star in time-scales between 1 and 10 Gyr. Conservation of angular momentum implies that the orbital infall of the planet also causes a decrease in the rotational period of the star. Given the large planet mass adopted for these runs, this effect is able to counteract the magnetic braking and at some point during the system’s evolution P_rot peaks and starts to decrease in value. Although this behaviour is well known (e.g. Ferraz-Mello et al. 2015), its effect on the distribution of hot Jupiters in the P–P plane has yet to be explored.

A comparison between the left and center plots seem to indicate that the evolution of the system is only weakly dependent on the star’s initial rotational period. The orbital evolution of the planet appears virtually identical in both cases, except for a change in the time-scale of the order of 0.5 Gyr, The values of P_rot are of course different, but both tend to very similar values after ∼5 Gyr. However, a more complex dependence is noted with respect to the initial eccentricity.

The thin continuous lines in the right-hand side plot of Fig. 2 shows the results of two simulations with the same initial values for P_orb and P_rot (indicated on top) but different initial eccentricities. While the general trend is similar in all cases, the time-scale and the maximum value attained by P_rot varies significantly, even though the planet reaches a synchronous state with almost circular orbit before 0.2 Gyr. This result is interesting since it indicates very different outcomes and time-scales even though most of the system’s evolution occurs in circular orbit and dominated by stellar tides.

An explanation may be found precisely during the road towards synchronization. Since the planet begins with P_orb < P_rot, the early orbital decay rate is much stronger than that expected for Ω₁ ∼ n. This may be observed in the behaviour of P_orb during the first stage of the system’s evolution. Once synchronization is reached, the orbital decay levels out and the subsequent decay is much shallower. The broad light-coloured lines show the evolution of ‘equivalent’ systems, characterized by e = 0, Ω₁ = n and initial semimajor axis a_equiv chosen such that both the planet’s orbital and star’s rotational evolution mimics the original system after synchronization. For initial rotation rates such that Ω₁ ≪ n we found |$a_{\rm equiv} \simeq a_{\rm ini} \left(1-e_{\rm ini}^2 \right)$|⁠, compatible with the conservation of the orbital angular momentum.

The overlap between the evolution of eccentric systems and their equivalent counterparts has far-reaching consequences. The most obvious is that we can simulate the tidal/braking interaction assuming circular orbits and synchronized planets as long as we adopt a_equiv as the initial semimajor axis. The tidal equations are thus simpler, the number of differential equations are reduced and the numerical integrations run much faster.

The second consequence is more relevant to our study and, perhaps, more debatable. Basically, we may say that when analysing the origin of the observed distribution of exoplanets in the P–P diagram, we need not be concerned about the primordial distribution of semimajor axes and eccentricities, but solely the distribution of a_equiv. We can thus approach our problem using the simplified tidal equations as long as we keep in mind that the initial separation between star and planet must be considered as representative and not equal to the true primordial value.

3 EVOLUTION OF SYSTEMS WITH A CLOSE-IN GASEOUS GIANT

We can now proceed to analyse whether our simple dynamical model can explain the distribution of planets in the P–P diagram, and see what tidal parameters are most suited for such a process. We begin studying hot Jupiters around Solar-type stars. As shown in Fig. 1, the distribution of large (R > 10R_⊕) KOIs taken from Santos et al. (2019, 2021), defines a moraine-like accumulation with a positive correlation between P_orb and P_rot, at least for close-in planets (or planetary candidates) with orbital periods up to 5 d. The only exception is KOI 554 which, as discussed previously, is probably a false positive. The same overall distribution is also found in older data, such as McQuillan et al. (2013) although perhaps less streamlined.

There are currently four different mechanisms proposed to explain the origin of hot Jupiters. Tidal capture following eccentricity excitation from planet–planet scattering (Beaugé & Nesvorný 2012) or perturbations from a stellar companion (Naoz, Farr & Rasio 2012), disc-induced planetary migration (Crida & Batygin 2014), and secular chaos (Wu & Lithwick 2011). The distribution of misalignment angles and planet multiplicity seem to indicate that no single process acted alone and at least part of the known population is a consequence of disc–planet interactions while the rest may be traced to orbital circularization after tidal capture (Dawson & Johnson 2018). Among the different dynamical processes leading to tidal capture, the distribution of misalignment angles seems to favor planet–planet scattering (Martí & Beaugé 2015) which may have occurred shortly after the dissipation of the gaseous disc (Lega, Morbidelli & Nesvorný 2013; Izidoro et al. 2021).

Both disc–planet interactions and post-dissipation instabilities point towards an early accumulation of hot Jupiters in their current orbital distance and, consequently, we can assume that the central star was still a rapid rotator at that time. Although the slow process of secular chaos may not be ruled out, it is currently difficult to evaluate how much it may have contributed to the observed distribution. Thus, even though we cannot rule out that some hot Jupiters may have reached their current orbital distance after magnetic braking drove the star to a slow rotation rate, there is little evidence to indicate that their number was substantial.

The three frames of Fig. 3 show, for different values of the stellar relaxation factor γ₀, the evolution of 4000 initial conditions in the P–P diagram, all consisting of a Jupiter-size planet in circular orbit and initial orbital period in the range P_orb ∈ [0.5, 5] d. The central star was assumed Solar-type with an initial stellar rotation P_rot = 1 d. The orange lines show the evolutionary tracks of four characteristic initial conditions; their starting orbital periods were chosen to be 2, 2.5, 3, and 4 d, respectively. Since we are assuming initial circular orbits and synchronous planets, the initial values of P_orb correspond to the equivalent semimajor axes, as defined in the previous section.

The arrows indicate the flow of the system. The characteristic evolution time-scale for stellar rotation, driven by magnetic braking, is much shorter than the evolution time-scale for the orbital period of the planet that is driven by tides. Therefore, even for an initial stellar rotation period shorter than the initial orbital periods, the evolutionary tracks begin as almost vertical straight lines. This stage, however, is temporary and after 10⁸−10⁹ yr the system reaches the domain of the accumulation highlighted in Fig. 1. For very small semimajor axes the rotation of the star is then accelerated due to the transfer of orbital angular momentum to the rotation of the star by means of the tides raised by the planet on the star. The planet falls on the star in a few Gyrs.

For wider systems, the tidal effects are not strong enough to cause the infall of the planet, but the evolution almost stops after ∼6 Gyr and the systems no longer shows a significant evolution in the P–P diagram. These hot Jupiters accumulate in a moraine-like structure. Like in the geological process with this name, big planets are carried along with a flow to the domain where the flow becomes weaker and accumulate there.

Since the evolutionary tracks of different initial conditions do not cross, we may colour-code the P–P diagram indicating the age at which the system reached a given spot. Lighter tones correspond to early stages in the evolution, while the dark brown region highlights the position attained by all initial conditions at times between 1 and 10 Gyr, the estimated range of ages for these systems (Lanza & Shkolnik 2014). To further aid in following the time evolution of the systems, the black lines show the isochrones for log₁₀(T) = 7, 8, 9, 10, where the time is given in years.

These results show that our simple dynamical model (tidal evolution + magnetic braking) lead to a distribution of Jupiter-size planets in the P–P plane that strongly resembles the observed distribution of hot Jupiters. The moraine shape is a natural consequence of the interplay between both phenomena on the stellar rotation, and may be used to estimate the value of the stellar relaxation factor γ₀ that appears to lead to a better correlation. As shown in the top frame of Fig. 3, a value of γ₀ = 2 s⁻¹ generates an excessively efficient tidal decay and a significant portion of the observed hot Jupiters lie above the isochrone associated to log₁₀(T) = 10. The larger relaxation factor considered in the middle frame leads to a much better fit, with practically all the population embedded in the region between T = 10⁹ yr and T = 10¹⁰yr. A similar conclusion may be drawn from the lower plot where an even larger value of γ₀ is considered.

We may thus deduce from these graphs that the observed distribution of hot Jupiters is consistent with a tidal evolution dominated by stellar tides with a relaxation factor between 10 and 50 s⁻¹. As the orbital distance of the planet decreases, so does the stellar rotation, leading to a moraine-type accumulation in the P–P diagram. Consequently, for massive close-in planets the stellar rotational period does not necessarily grow monotonically with time and care must be taken when using P_rot as a direct proxy for stellar age.

As we mentioned at the beginning of this analysis, most of the hot Jupiters are expected to have reached their present location before the stellar rotation slowed significantly. Even if this was not the case, we have found that the initial rotational period of the star exerts almost no influence on the evolution. As shown in the simulations carried out in Fig. 2, even a relatively slow initial rotation of 10 d led to the same evolutionary tracks in the P–P plane after ∼1 Gyr. Even the final system ages for a given value of P_orb or P_rot only varied by ∼0.5 Gyr at the end of the simulations. We therefore believe that the above analysis should be fairly robust and independent of the initial conditions.

4 EVOLUTION OF SYSTEMS WITH A CLOSE-IN SUB-JOVIAN PLANET

Regardless of which stellar relaxation factor better fits the data, it seems that the moraine is (at least partially) caused by a funneling of the dark-toned region as P_orb → 0. In turn, this effect is associated to a change in sign of the time derivative of the star rotation P_rot. Regardless of the age of the system when this occurs, as soon as the tidally induced speed-up of the stellar rotation surpasses the slow-down caused by magnetic braking, the evolutionary tracks converge and lead to an orbital infall of the planets in a tight formation.

To understand when this phenomena occurs and under what system parameters, we look for solutions of the equation:

where the first term is the derivative of the stellar spin Ω₀ due to stellar tides and the second is the contribution from magnetic braking. Full expressions for both are given in Appendix  A. Assuming |ν|/γ₀ = 2|n − Ω₀|/γ₀ ≪ 1 and n ≫ Ω₀, we can approximate the tidal term by

Similarly, since we expect the moraine to occur for slow rotators where the magnetic braking is not so efficient, for this region we can approximate the second term by:

Introducing both expressions into (equation 1) we obtain that the maximum stellar rotational period |$P^{(\rm max)}_{\rm rot}$|⁠, before tidal effects begin to dominate, is approximately given by:

Relating a and P_orb through Kepler’s third law, and taking base 10 log on both sides, we finally obtain

where Λ is a constant that depends only on the physical properties of the system, and is approximately given by

We thus obtain a linear relation between both periods in a log–log scale; moreover for any given value of P_orb the critical value of the stellar rotation period should increase for smaller planetary masses and larger stellar relaxation factors.

Fig. 4 shows the critical value of P_rot as a function of the orbital period, for four different planet types; each is identified by the text alongside the colour line. Masses were taken from Table 1 while the star was assumed Solar-type with α₀ = 0.07 and |$k_{\mathrm{ f}_0} = 0.05$|⁠. Continuous lines correspond to γ₀ = 10 s⁻¹, while their dashed counterparts were obtained assuming γ₀ = 50 s⁻¹.

Jupiter-type planets, characterized by large masses, generate strong tidal interactions with the star, thus leading to relatively small values of |$P^{(\rm max)}_{\rm rot}$| even for large relaxation factors. This is consistent with the existence of the moraine-like distribution of hot Jupiters in the P–P diagram, especially noticeable for orbital periods below P_orb ∼ 3 d.

The case of smaller gaseous/icy planets (Saturn and Neptune) is less clear. Although the results of Fig. 4 seem to predict a moraine-type distribution, at least for Saturn-size bodies and for very small orbital periods, the observability also depends on the distribution of planets. This is analysed in Fig. 5 for two different values of γ₀. As was done in Fig. 3 for hot Jupiters, the orange curves show the evolutionary curves of six initial conditions with different P_orb and the same stellar rotation period P_rot = 1 d. Orbits are assumed circular and with n = Ω₁.

The distribution of confirmed planets and Kepler candidates is again shown with filled circles whose colours are indicative of the effective temperature of the star. The paucity of bodies in this size range with orbital period less than 2–3 d is well known (Beaugé & Nesvorný 2013) and usually referred to as the sub-Jovian desert. It is interesting to note that most of the candidates at this orbital distance belong to cool stars with T_eff ≲ 5500 K. Consequently, in our simulations of tidal + braking evolutions we assumed a central star with m₀ = 0.8m_⊙, R₀ = 0.77R_⊙, and |$\alpha _0 = k_{\mathrm{ f}_0} = 0.05$| (Claret 2019).

Given the lack of sub-Jovian planets very close to the star, it is difficult to correlate their distribution with our dynamical model, let alone discuss which tidal relaxation factor better fits the data. However, except for a single system with P_orb ≃ 1 d and orbiting a fast rotator, the rest of the population does not show any inconsistency with the expected evolutionary tracks, nor does it seem necessary to include additional phenomena into the model.

Finally, Fig. 6 shows analogous results, but now focusing on small planets (R ≤ 3R_⊕). Since the known population is large, we restricted the analysis to Solar-type stars. We verified that there is no significant difference with the results obtained from bodies around cooler or hotter stars. Again we assume that the close-in planets reached this region when the star was still a fast rotator. As with the hot Jupiters described previously, the main scenarios for the origin of these systems is planetary migration and planet-scattering following disc dispersal (e.g. Izidoro et al. 2021). Both processes are believed to have occurred before magnetic braking had the opportunity to slow the stellar rotation significantly.

As expected from the values of |$P^{(\rm max)}_{\rm rot}$| predicted from Fig. 4, no moraine-type structure is observed; moreover the distribution in the P–P diagram appears almost flat for time-scales between 10⁹ and 10¹⁰ yr.

As discussed recently by Messias et al. (2022), we also find no relevant evidence of the lower bound for orbital periods, at least in the form suggested by McQuillan et al. (2013). The rapid slow-down of P_rot due to magnetic braking and the tidal decay for very short period planets seems to account for the observed distribution quite well. Both values of γ₀ lead to similar outcomes, although the shape of the tracer orange curves seems like a better fit for γ₀ = 10 s⁻¹. For orbital periods P_orb ≳ 3 d, tidal effects are negligible and the rotation of the corresponding stars evolve following only the rotational braking due to stellar winds. The time-scale in this case follows closely the Skumanich law (Skumanich 1972).

5 THE STAR’S RELAXATION FACTOR

In our simulations for Jupiter-size planets, we considered three different values for the stellar relaxation factor γ₀: 2, 10, and 50 s⁻¹ (see Fig. 3). They show that, at least for Solar-type stars, a value of γ₀ between ∼10 and ∼50 s⁻¹ appears to better represent the possible evolution of these systems to reach the moraine-like domain. Smaller values lead to evolutionary tracks in which the hot Jupiter falls on the star before the system reaches the moraine-like accumulation.

This result may be compared to those obtained by Ferraz-Mello et al. (2015) for the rotation of the host stars in several systems with large companions (exoplanets or brown dwarfs). There, the result was rather centered on γ₀ = 50 s⁻¹ due to a value for the fluid Love number |$k_{\mathrm{ f}_0}$| some five times larger than the value adopted here.

While the distribution of smaller planets did not allow to further constrain γ₀, the observed structures are consistent with the values derived from hot Jupiters. Similarly, the distribution of sub-Jovian planets around cooler stars do not appear to show any significant difference with respect to that expected for Solar-type stars. We thus found no evidence that the relaxation factor of stars could be a strong function of the stellar type.

It is interesting to see how the above results are translated in terms of the stellar quality factor Q₀ (or its variant |$Q^{\prime }_0=1.5 Q_0/k_{\mathrm{ f}_0}$|⁠) used in some current tidal friction theories. Using the relation given in Appendix A2 and supposing |ν| ≪ γ₀, we may write

We obtain the curves shown in Fig. 7. We note the big variation of Q₀ along the path in all solutions making evident that the choice of Q₀ to parametrize the dissipation in evolving systems is not a good one. The definition of Q₀ mixes a property of the body, the relaxation, with two frequencies of the system, Ω₀ and n. Finally, we omitted in Fig. 7 the parts of the paths in which Ω₀ < n and the neighbourhood of the synchronization where the definition of Q₀ becomes singular.

6 CONCLUSION

This paper shows that the distribution of the points representing the short period Kepler systems and KOIs (P < 5 d), for which the period of rotation of the star and the orbital period could be determined, presents some features that may be explained by the joint action of the tidal evolution of the system and the magnetic braking of the star. The main features in the two periods (P–P) diagram are: (1) The existence of a moraine-like accumulation of systems hosting a hot Jupiter with orbital period in the range 1.5−5 d, all in an inclined zone around the stellar rotation period 25 d. (2) The absence of any correlation between orbital period and stellar rotation in the case of small (R ≲ 3R_⊕) planets.

The creep tide theory allowed us to calculate the evolution of systems with one exoplanet in orbit around a Sun-like star showing that systems hosting a hot Jupiter with orbital period shorter than 5 d have a fast evolution upwards in the P–P diagram and stops evolving exactly in the moraine-like accumulation seen in the diagram or, if the initial orbital period is much shorter, fall on the star. In the case of systems hosting small planets, tidal evolution quickly pulls down planets with orbital periods initially smaller than 1.5–2 d; those with periods slightly larger first evolve upwards thanks to the magnetic braking and then slowly fall towards the star. The distribution of these systems in the P–P diagram is consistent with a flow ruled by tidal evolution and magnetic braking.

The boundaries of distributions seen in the P–P diagram are determined by the intensities of these two agents. In the case of the tidal evolution, we have found that the most probable relaxation factor of the Sun-like star lies between 10 and 50 s⁻¹. If the dissipation is much larger or much smaller than these values, the accumulation would not be located in the place where they are observed in the P–P diagram constructed with either the results of McQuillan et al. (2013) or those of Santos et al. (2019) and Santos et al. (2021). This result is in agreement with the values obtained for the stars relaxation factor by Ferraz-Mello et al. (2015), if we take into account that the fluid Love number used there (0.26) was some five times larger than the value used in this paper (0.05) and that, for gaseous bodies, the variation of the orbital elements is proportional to k_f/γ. In this paper, we used a determination of the Sun’s k_f obtained with the density profile of the standard solar model and the formulas given by Folonier et al. (2015) to compute the fluid Love number of layered non-homogeneous bodies.

The use of so-called equivalent initial semimajor axis a_equiv in our simulations allowed us to avoid the lack of information regarding the primordial eccentricities of the planets as they arrived to the region close to the star. Finally, the lack of a moraine-like structure or any significant correlation between P_rot and P_orb for low-mass planets is also in accordance with our model, and seems additional and strong evidence that tidal interactions, together with magnetic braking, have probably been the driving forces behind many of the observed dynamical features of close-in planets and their host stars.

ACKNOWLEDGEMENTS

The authors wish to thank the insightful analysis of an anonymous referee that helped improve this work, both in content and presentation. SFM acknowledges the support of CNPq (Proc. 303540/2020-6) and FAPESP (Proc. 2016/13750-6 ref. Mission PLATO). CB is grateful to SECYT/UNC for financial support.

DATA AVAILABILITY

All data analysed in this paper is publicly available, but a copy may be requested from the authors.

References

APPENDIX: THE EQUATIONS OF MOTION

A1 Tidal variation of the elements

The equations used to calculate the tidal variations of the orbital elements: semimajor axis and eccentricity, consider the density profile of the bodies, but imposes some constraints in order to avoid a huge number of free parameters. They assume that the body is formed by co-rotating homogeneous layers (Ferraz-Mello, Folonier & Gomes 2022). They are

and

where |$\overline{\epsilon }_\rho$| is the mean flattening of the equivalent Jeans homogeneous spheroid:

In this expression m is the mass of the tidally deformed body, M is the mass of the companion whose attraction is creating the tidal potential, R is the radius of the deformed body, and a is the semimajor axis of the relative orbit of the two bodies. No hypothesis is done on the relative size of the two masses. In fact, in the general case, we have to consider the tides raised in the star by the planet and also the tides raised in the planet by the star, so m ≪ M and M ≪ m. In both cases we use the same equations and we have added their contributions to the variations of the orbital elements to get the total effect.

In the applications described in this paper, both were considered, but the contribution of the tides on the exoplanet are too small to be significant. The role of the tides on the exoplanet is to quickly drive the rotation of the planet to the stationary almost synchronous state where it remains trapped for the rest of the time (see Fig. 2).

The tidal equations involve some known functions:

• The Cayley functions of the orbital eccentricity e:

  $$\begin{eqnarray} E_{q,p}(e) = \frac{1}{2\pi }\int _0^{2\pi }\left(\frac{a}{r}\right)^3\cos {(qv+(p-q)\ell)\ \mathrm{ d}\ell }, \end{eqnarray}$$

  (A4)

  where r is the modulus of the position vector, ℓ is the mean anomaly, and v is the true anomaly. For small eccentricities, we may use the polynomial approximations published by Cayley (1861) or, equivalently, the Hansen coefficients (Hughes 1981).
• The other functions are

  $$\begin{eqnarray} \sin {2\sigma _{k}} &=& \frac{2\gamma (\nu +kn)}{{\gamma ^2+(\nu +kn)^2}}; \\ \sin {2\sigma ^{\prime \prime }_{k}} &=& \frac{2\gamma kn}{{\gamma ^2+k^2 n^2}}; \end{eqnarray}$$

  (A5)

  where γ is the relaxation factor and ν = 2Ω − − − 2n is the mean semidiurnal frequency (Ω is the rotational velocity of the deformed body and n is the orbital mean motion).

The other parameter appearing in the variation equations is the fluid Love number, k_f. This parameter is related to the mass concentration of the body and may be determined from a model of the density profile of the body (Folonier et al. 2015). It is equal to twice the apsidal motion constant introduced in the study of close binary stars (see Kopal 1953). In the case of the present-day Sun, the standard solar model leads to a value close to 0.05. In previous applications, a rule using the relationship between the fluid Love number and the moment of inertia (k_f = 15C/4mR²) was used; however, this relation, valid for homogeneous bodies, overestimates the value of k_f in the cases under study in this paper.

The variation of the angular velocity of rotation of the deformed body is given by

where C is the moment of inertia of the deformed body and G the gravitation constant.

A2 Tidal dissipation

The tidal energy released inside the star is negligible when compared to the thermal energy produced by the hydrogen burning near the center of the star. However, many authors use the quality factor Q to parameterize the tidal effects and the equation giving the tidal dissipation is necessary to allow the comparison of the relaxation factor of the creep tide theory to the quality factor used in several current tidal friction theories.

In the creep tide theory, the variation of the mechanical energy is

or, in a first approximation,

(N.B. E_{2, 0} ≃ 1) (see Ferraz-Mello et al. 2022). This equation is the same used in other theories for motions far from the synchronization when we introduce the quality factor Q through

as discussed, for example, in Efroimsky & Lainey (2007).

A3 Braking of the star rotation due to stellar activity

The magnetic braking of the star rotation was computed using the results of Bouvier et al. (1997) for stars with masses in the range 0.5M_⊙ < M_s < 1.1M_⊙:

where ω_sat is the value at which the angular momentum loss saturates (fixed at |$\omega _{\rm sat} = 3, 8, 14 \, \Omega _\odot$| for 0.5, 0.8, and 1.0 M_⊙ stars, respectively) and B_W is a factor depending on the star moment of inertia, mass, and radius through the relation

(Bouvier et al. 1997). The subscript s is used to stress the fact that the star parameters are being considered.

The above form of the law is valid after the star has completed its contraction (the stellar moment of inertia C_s no longer changes significantly) and is fully decoupled from its primeval disc. F-type stars are not expected to be affected by the magnetic braking.