© ESO 2021

1 Introduction

Planets with very short orbital periods have been discovered by observing their transits across the disks of their host stars. Considering those with an orbital period P≲1 day, Sanchis-Ojeda et al. (2014) found an occurrence rate ranging from 0.15 ± 0.05 percent for F-type dwarf stars to 1.1 ± 0.4 percent for M-type dwarfs with a prevalence of planets with radii between 0.8 and 1.25 Earth radii. For those planets, the occurence rate appears to be similar up to P≲2 days (Fressin et al. 2013; Dressing & Charbonneau 2015).

These planets are amenable to a measurement of their mass through radial velocity techniques because they induce wobbles of their host stars with amplitudes of a few m s⁻¹ that are detectable with current spectrographs, while a correction of the stellar jitter due to magnetic activity can usually be performed (e.g., Malavolta et al. 2018). In some cases, they are also detectable in reflected light thanks to their closeness to their hosts, thus allowing a study of their atmospheric thermal emission and distribution of their surface brightness (see, e.g., Essack et al. 2020, for details). Such studies are not presently possible for telluric planets orbiting at larger separations from their host stars.

Owing to the small distance from their host stars, such planets interact strongly with them through the gravitational, the radiation, the magnetic, and the stellar wind fields. Gravitation is responsible for very strong tidal interactions affecting the orbit and the rotation of these planets (e.g., Mathis 2018), while the intense irradiation produces surface temperatures of thousands of kelvins on the illuminated hemisphere as well as a strong atmospheric evaporation powered by the stellar flux in the extreme ultraviolet (Owen 2019). Magnetic interactions are also possible as recently discussed by, for example, Strugarek et al. (2019) in the case of Kepler-78b, a planet with an orbital period of only 8.5 h.

The close separation between these planets and their host stars leads to a synchronization of the planet rotation and a damping of the orbit eccentricity on timescales shorter than 1–10 Myr owing to the strong tidal interactions (see Sect. 2). If the orbit eccentricity is not excited by other planets in the system, the internal heating of these planets is dominated by the slow release of heat stored in the planet core after their formation and by the decay of radioactive elements (e.g., Nimmo et al. 2020). However, another internal heating mechanism can operate during short phases of the planet evolution and this work is dedicated to introduce it.

Young solar-like stars are usually fast rotators and show an intense magnetic activity sustained by a vigorous hydromagnetic dynamo. In these stars, strong magnetic fields with intensities up to several tens of teslas can be amplified and stored in the overshoot layers beneath their convection zones where such strong fields are probably organized in slender flux tubes that are stable against magnetic buoyancy thanks to the subadiabatic stratification of the layers. When their magnetic field intensity reaches a critical value, they become unstable to undulatory instability and emerge forming bipolar active regions in the photosphere characterized by the appearance of starspots (e.g., Caligari et al. 1995; Granzer et al. 2000, and references therein). These strong-field flux tubes have a slightly lower density than the surrounding medium and are not distributed symmetrically in longitude over the star. Therefore, they produce a slight deviation of the stellar mass distribution from an axisymmetric configuration that in turn leads to a non-axisymmetric quadrupole component in the outer gravitational field of the rotating star (see Sect. 2).

When the stellar rotation period is twice the orbital period of the planet, such a quadrupole moment produces a resonant coupling between the stellar rotation and the orbital motion leading to an oscillation in the star-planet separation. As a consequence, a remarkable tidal dissipation is induced into the planet leading to an internal heating that can profoundly affect its internal structure (see Sects. 2 and 3). Although the spin-orbit resonance has a limited duration because of the evolution of the stellar rotation, the energy released inside the planet can produce a complete melting of its interior that may enhance its differentiation and affect the chemical composition of its surface and its atmosphere,thus leading to potentially permanent changes that could be observable (see Sect. 4).

2 Model

In this section, we introduce the different ingredients of our model starting from the effects of strong magnetic flux tubes on the internal density stratification and the outer gravitational potential of the host star (Sect. 2.1); discuss the tidal deformation of the star and the planet (Sect. 2.2); write the Lagrangian function of the star-planet system from which we derived the equations of motion for its components (Sects. 2.3 and 2.4); and include the tidal damping of the orbital eccentricity (Sect. 2.5). Our equations are then solved focusing on the radial orbital motion (Sects. 2.6) and the libration of the planet (Sects. 2.7 and 2.8). Finally, a simple model for the evolution of the stellar rotation and the orbital semimajor axis is introduced (Sect. 2.9) in order to show how a resonance between the stellar rotation and the radial motion of the planet can happen along the evolution of the system. During such a resonance, a remarkable tidal dissipation occurs inside the planet as we shall see in some illustrative examples in Sect. 3.

2.1 The outer gravitational field of a magnetically active star

Late-type main-sequence stars have an internal structure characterized by an outer convection zone and a radiative interior, similar to the Sun. Immediately beneath their convection zone, there is a layer characterized by undershooting convection that is called the overshoot layer by analogy with a similar layer present in upper main-sequence stars (Zahn 1991). Such a layer is characterized by a subadiabatic stratification and a strong radial shear in the Sun (Schou et al. 1998) that allow to store and amplify intense magnetic fields. They reach intensities between 10 and 200 T, depending on the depth of the convection zone and the rotation period of the star before becoming unstable and emerge through the convection zone. Owing to the conservation of the angular momentum, the flux tubes do not emerge keeping their initial circular symmetry around the rotation axis (see, e.g., Caligari et al. 1995, for a discussion of the force balance and the geometry of flux tubes initially in equilibrium), but become non-axisymmetric with the azimuthal mode m = 1 usually being the first to become unstable in rapidly rotating stars as differential rotation amplifies the field strength (van Ballegooijen 1982a; Moreno-Insertis et al. 1992; Ferriz-Mas & Schüssler 1994; Granzer et al. 2000). After emergence, these flux tubes are characterized by a strong buoyancy that makes them vertical with their roots anchored into the overshoot layer (van Ballegooijen 1982b; Schüssler et al. 1994; Caligari et al. 1995).

The density inside these vertical flux tubes is lower than in the surrounding medium because the external pressure is balanced in part by the interior magnetic pressure. The equation governing the transversal mechanical equilibrium of a slender flux tube is (e.g., van Ballegooijen 1982b, and references therein) $$p_{\text{i}}(r)+\frac{B^2(r)}{2\mu }=p_{\text{e}}(r),$$(1)

where p_i(r) and p_e (r) are the plasma pressures inside and outside the flux tube, respectively, while B(r) is the magnetic field strength at the radialdistance r from the center of the star and μ the magnetic permeability of the plasma. By differentiating Eq. (1) with respect to the radius r and using the equation for the radial pressure gradient dp_i, e∕dr = −g(r) ρ_i, e(r), we find $${\rho }_{\text{i}}(r)-{\rho }_{\text{e}}(r)=\frac{1}{2\mu g(r)}\frac{\text{d}{B}^2(r)}{\text{d}r},$$(2)

where g(r) is the acceleration of gravity and ρ_i, e(r) the density of the plasma inside and outside the magnetic flux tube, respectively. Equation (2) allows to compute the density deficit inside a slender magnetic flux tube when a model for the radial variation of the field inside the tube itself is specified. The simplest approach is to compute B(r) by assuming that the field is mainly vertical due to the strong magnetic buoyancy and that ∇ ⋅ B = 0. Such an approximation fails in the outermost layers of a star, but this has a very limited impact on the evaluation of the stellar quadrupole moment because those layers have a much lower density than the interior of the convection zone. We refer the reader to Sect. 2.3 of Lanza (2020a) for details on the calculations. A qualitative sketch of the effect of the density deficit in a vertical magnetic flux tube is given in Fig. 1 where the quadrupole gravitational potential associated with such a density deficit is illustrated.

The light curves of magnetically active stars and the Doppler imaging of their photospheres reveal that their magnetic flux tubes are unevenly distributed versus the stellar longitude (see Sect. 2.1 of Lanza 2020a, for details). In addition to observations, hydromagnetic dynamo models also show the predominance of non-axisymmetric modes with an azimuthal wavenumber m = 1 in stars rotating several times faster than the Sun (e.g., Käpylä et al. 2013; Viviani et al. 2018; Käpylä 2021). Therefore, the net effect of several magnetic flux tubes simultaneously present in the convection zone of an active star is that of producing a non-axisymmetric gravitational quadrupole moment when they are preferentially concentrated in one hemisphere as indicated by the observations and suggested by the above dynamo models. Sometimes, configurations with two starspots of similar sizes on opposite hemispheres are observed (e.g., Collier Cameron et al. 2009), that point to an instability of magnetic flux tubes dominated by the m = 2 mode as suggested by models in some specific regimes of stellar rotation and depth of the convection zone (e.g., Granzer et al. 2000). In such a case, the gravitational quadrupole moment of the star is reduced, but the octupole moment could be significant and produce a similar resonance effect when the first harmonic of the rotation period corresponds to twice the orbital period of the planet. For simplicity, we shall not explore this possibility in the present work and limit ourselves to a model dominated by a non-axisymmetric density perturbation associated with the azimuthal m = 1 mode. The single flux tube geometry depicted in Fig. 1 is adopted as the simplest configuration to model this case. The effects of several flux tubes simultaneously present in the stellar convection zone can be computed by summing the contributions of the individual tubes because their effects on the quadrupole moment tensor are simply additive (see Eqs. (4) and (5) below).

The outer gravitational potential Φ of the active star at a given point P can be expressed as: $$\Phi =-\frac{G{m}_{\text{s}}}{r}-\frac{3G}{2r^3}{\sum }_{i,k}\frac{Q_{ik}{x}_i{x}_k}{r^2},$$(3)

where G is the constant of gravitation, m_s the mass of the star, r the distance of P from the barycenter O of the star, Q_ik the quadrupole moment tensor of thestar, and x_i the coordinates of P in a Cartesian reference frame of origin O, while the indexes i, k = 1, 2, 3 specify the Cartesian coordinates. The quadrupole moment tensor can be expressed in terms of the inertia tensor of the mass distribution of the star as $$Q_{ik}=I_{ik}-\frac{1}{3}{\delta }_{ik}\text{Tr}I,$$(4)

where δ_ik is the Kronecker δ tensor and TrI = I_xx + I_yy + I_zz is the trace of the inertia tensor I, that is, the sum of its diagonal components. The components of the inertia tensor are given by $$I_{ik}={\displaystyle {\int }_V\rho }(x)x_i{x}_k\text{d}V,$$(5)

where ρ is the density, x the position vector, and V the volume of the star over which the integration is extended.

Equation (3) can be simplified by adopting a Cartesian reference frame whose axes are directed along the principal axes of inertia of the star so that only the diagonal components of the inertia tensor are different from zero. The principal axes of inertia are generally coincident with the axes of symmetry of the star. We assume that the z-axis is directed along its rotation axis, while the x- and y-axes are in the equatorial plane of the star. Introducing a set of spherical polar coordinates (r, θ, α), where r is the distance from the barycenter O of the star, θ the colatitude measured from the z-axis, and α the longitude, we have: $$\Phi (r,\theta ,\alpha )=-\frac{G{m}_{\text{s}}}{r}-\frac{3G}{2r^3}\left(-2Q+3Q{\sin }^2\theta +\frac{1}{2}T{\sin }^2\theta \cos 2\alpha \right),$$(6)

where the two scalars Q ≡ (Q_xx + Q_yy)∕2 = −Q_zz∕2 and T ≡ Q_xx − Q_yy express the axisymmetric and non-axisymmetric components of the quadrupole tensor Q_ik referred to its principal axes, respectively. In that reference frame, only two scalars are required to specify the tensor because it is traceless by definition (cf. Eq. (4)).

Assuming that the magnetic flux tube in Fig. 1 lies in the equatorial plane of the star along the y-axis, the principal axes of inertia will be the y-axis itself and the two perpendicular axes x in the equatorial plane and z orthogonal to the equatorial plane. By applying Eq. (5), we obtain $I_{xx}=2m_{\text{A}}{x}_{\text{A}}^2$, I_yy = I_zz = 0, and Tr I = I_xx, where m_A and x_A are the mass and the distance from the barycenter of each of the two point masses A and A′ that simulate the effect of the density deficit inside the magnetic flux tube. In this way, we find $Q=(1/3)m_{\text{A}}{x}_{\text{A}}^2$ and $T=2m_{\text{A}}{x}_{\text{A}}^2$. Substituting these expressions into Eq. (6), we obtain the outer gravitational potential of the star in this simple approximation.

The quadrupole term coincides with the expression obtained by developing the potentials of the individual point masses m_A at the outer point P in series of the small ratio (|x_A|∕r) ≪ 1 up to the second order (cf. Murray & Dermott 1999, Sect. 5.3). It shows that the principal axis of inertia x coincides with the line joining the two point masses that is also the axis of symmetry of the bulge of the density distribution of the star. We note that the moment of inertia for rotation about the x-axis, that is, I_yy + I_zz, is minimal.

Going beyond this simplified model, the quadrupole terms produced by a vertical flux tube in the equatorial plane have been computed by Lanza (2020a) following a method based on Eqs. (1) and (2) (see Sects. 2.3 and 2.4 of Lanza 2020a, for details). Considering the configuration depicted in Fig. 1 with the axis of the flux tube directed along the y-axis of the reference frame and the angle α measured from the x-axis, they are $$Q=\frac{\pi }{3}\frac{B_0^2{r}_{\text{b}}^4}{\mu }J{\sin }^2{\theta }_0\cos {\theta }_0,$$(7) $$T=2\pi \frac{B_0^2{r}_{\text{b}}^4}{\mu }J{\sin }^2{\theta }_0\cos {\theta }_0,$$(8)

where B₀ is the intensity of the magnetic field at the base of the flux tube in the overshoot layer that is at a radial distance r_b from the center of the star, θ₀ the angular radius of the flux tube, and $J$ is the integral $$J={\displaystyle {\int }_{r_{\text{b}}}^{r_{\text{L}}}\frac{r^{\prime }}{GM(r^{\prime })}}\text{d}{r}^{\prime },$$(9)

where M(r) is the mass of the star inside the radius r and r_L the limit radius beyond which the magnetic field of the flux tube deviates significantly from the vertical because the pressure of the plasma outside the flux tube becomes too small to confine the field. The exact value of r_L is not critical because the density in the outer layers of the star where the field is no longer confined by the plasma pressure is so small that the effect on the quadrupole moment is negligible (cf. Lanza 2020a). If the fraction of the stellar surface covered by the section of the magnetic flux tube is f_s, the angle θ₀ is given by θ₀ = arccos(1 − 2f_s) or ${\theta }_0\simeq 2\sqrt{f_{\text{s}}}$ for f_s ≪ 1.

In Lanza (2020a), the axis of the flux tube was directed along the x-axis of the reference frame in Fig. 1, thus there was a minus sign in the r.h.s. of Eqs. (7) and (8). Choosing to have the tube axis along the y-axis instead of along the x-axis is equivalent to add π∕2 to the angle α in Eq. (6) which compensatesfor the sign change of Q and T with respect to the expressions in Lanza (2020a). The present choice appears to be more motivated from a physical point of view because the density deficit inside the flux tube corresponds to adding two excess masses m_A along the x-axis perpendicular to the tube axis. With the simple two-mass model, we obtain the above quadrupole moment components by assuming that the effective moment of inertia $m_{\text{A}}{x}_{\text{A}}^2$ is equal to $\pi (B_0^2{r}_{\text{b}}^4/\mu )J{\sin }^2{\theta }_0\cos {\theta }_0$.

Fig. 1 Effect of the density perturbation inside a vertical magnetic flux tube FT (indicated by the purple stub) in a stellar convection zone. For the purpose of evaluating the outer stellar gravitational potential, the mass deficit inside the tube is equivalent to adding two point masses mA at the pointsA and A′ on a plane perpendicular to the tube axis here taken to be coincident with the y axis. The x-axis of our reference frame is taken along the line joining the two masses that are symmetrically situated with respect to the barycenter of the star O at the abscissas ± xA. These two masses produce a gravitational potential outside the star that is no longer spherically symmetric and can be described by a quadrupole term. It depends on the distance OP = r of a given point P from the barycenter O and the angle α between the directions OP and OA of the line joining the two point masses (see the text).

2.2 Origin of the stellar and planet quadrupole moments

The effect of the tidal potentials acting on the star and the planet is that of producing a deviation of their density distributions from axial symmetry. The outer gravitational potential of each body can be developed in a series of spherical harmonics the angular dependence of which is given by the Legendre polynomials Y_lm (e.g., Ogilvie 2014). The term of the tidal potential of degree l and azimuthal order m produces a density perturbation in each body that is responsible for a component of the outer gravitational potential with the same degree and azimuthal order, but with a different amplitude that depends on the internal stratification and the rigidity of the body.

The ratio of the amplitude of the outer gravitational component of degree l and order m to the amplitude of the corresponding component of the tidal potential is indicated with k_lm and is called the Love number of degree l and azimuthal order m of the body. It is in general a complex number the imaginary part of which expresses the tidal lag and depends on the frequency of oscillation $\widehat{\omega }$ of the tidal potential. The real part of k_lm expresses the hydrostatic deformation of the body under the action of the tidal potential and is generally much larger than its imaginary part because the tidal lag is usually a very small angle, especially in the case of fluid bodies such as a star. Since the most important component of the tidal potential corresponds to l = m = 2, we shall focus on the Love number k₂₂ = k₂, that is the second-order Love number. It can be used to express the outer quadrupole potential of a body produced by the equilibrium tide. It corresponds to an ellipsoidal deformation whose axis of symmetry is always directed along the line joining the barycenters of the two bodies (the perturbed one and the companion producing the tidal potential), if we neglect the very small tidal lag angle.

In the case of a fluid body, like a star or a giant planet, k₂ is entirely defined by the function ρ(r) that specifies the unperturbed density of the body as a function of the distance r from its barycenter and that is spherically symmetric (e.g., Kopal 1959). However, in the case of a telluric planet the value of k₂ depends also on the rigidity of the body that takes into account the deviation of its response from the purely hydrostatic case. Therefore, the quadrupole moment of a planet of mass m_p and with a radius R_p at a separation D from a star of mass m_s can be expressed as (cf. Remus et al. 2012) $$T=k_2\left(\frac{m_{\text{s}}}{m_{\text{p}}}\right){\left(\frac{R_{\text{p}}}{D}\right)}^3{m}_{\text{p}}{R}_{\text{p}}^2;$$(10)

for the Earth k₂ = 0.295 (Lainey 2016), while for a fluid body of uniform density k₂ = 3∕2 (Ogilvie 2014).

In the case of the host star in a planetary system, the hydrostatic deformation due to the tides raised by the planet produces an hydrostatic quadrupole moment that can be neglected for our purposes because the tidal bulge is always directed along the line joining the centers of the two bodies. The tidal lag angle α_T, that measures the deviation of the axis of symmetry of the bulge from the line joining the barycenters of the two bodies, is negligible. More precisely, it is given by ${\alpha }_{\text{T}}\simeq 1/Q_{\text{s}}^{\text{T}}$, where $Q_{\text{s}}^{\text{T}}$ is the tidal quality factor of the star (Ogilvie 2014)¹, that is on the order of at least 10⁴−10⁶ in the case of late-type stars (Ogilvie & Lin 2007, see also below Sect. 2.9).

The situation is different in the case of the non-axisymmetric quadrupole produced by the density perturbation inside the magnetic flux tubes in the stellar convection zone that we considered in Sect. 2.1. In that case, the bulge of the density perturbation is not constrained to stay close to the line joining the barycenters of the two bodies, so the angle α in Fig. 1 can assume any value.

The magnetic flux tubes we are considering are not permanent structures because their magnetic fields are subject to turbulent diffusion, thus they have a finite lifetime in stellar convection zones. In the case of the Sun, the longest lived magnetic flux tubes are associated with recurrent active regions that last a few solar rotations, that is, a few months (e.g., Solanki 2003; Hathaway & Choudhary 2008). Nevertheless, in the case of young active stars, active regions can persist for several years at the same active longitude because they are much larger and the diffusion timescale is proportional to the area of the region. Lifetimes of hundreds of days are observed in Sun-like stars with rotation periods P_rot shorter than 5−7 days (e.g., Giles et al. 2017), while active longitudes persist for decades in rapidly rotating stars (P_rot≲2−4 days) with deep convection zones showing a very strong dynamo action (e.g., Lanza et al. 2006). Therefore, we can regard these active stars as endowed with non-axisymmetric quadrupole moments that persist for timescales ranging from hundreds of days to several years. Only the dynamical effect of this component of the non-axisymmetric quadrupole moment will be considered in our model because the value of the angle α can change in time, while we shall neglect the quadrupole associated with the equilibrium tides because α_T is extremely small and thus gives a virtually constant contribution to the potential that can be neglected for our purposes (cf. Sect. 2.3).

Finally, we consider that a telluric planet may have a permanent quadrupole moment associated with its solid layers in addition to the tidal (hydrostatic) quadrupole deformation. In the case of the Earth, the hydrostatic deformation is prevalent and $Q_{\text{p}}^{\text{T}}~280$ considering the Earth-Moon tidal interaction (Lainey 2016); this implies that α_T is always small (see Sect. 2.5 for more information on $Q_{\text{p}}^{\text{T}}$ in the case of telluric planets). Conversely, in the case of the Moon or Mercury, the permanent quadrupole moment associated with the persistent deviation of their bodies from axisymmetry is prevalent and is responsible for the spin-orbit coupling experienced by those bodies (Goldreich & Peale 1968; Murray & Dermott 1999, Ch. 5). Therefore, we shall include the contributionof a possible permanent quadrupole moment of the planet in our model because the corresponding lag angle α is not a priori constrained to remain small.

2.3 Lagrangian function and equations of motion of the star-planet system

To study the dynamics of our star-planet system, we apply the Lagrangian formalism (e.g., Goldstein 1950). The Lagrangian $L$ of our system is defined as $$L=T-\Psi ,$$(11)

where $T$ is the kinetic energy of the system and Ψ its potentialenergy expressed as functions of the coordinates and their time derivatives in an inertial reference frame. We consider a reference frame having its origin at the barycenter Z of the star-planet system and the xy plane of which is the orbital plane of the system. In this reference frame, the kinetic energy of the system can be written as the sum of the kinetic energy of the relative orbital motion and of the kinetic energy of rotation of each of the bodies around an axis passing through its barycenter, respectively. For the sake of simplicity, we assume that the spin of the star and of the planet are perpendicular to the orbital plane. Tides inside the planet induce the alignment and the synchronization of its spin with the orbital motion on a timescale on the order of ≲10⁵ yr, much shorter than the lifetime of the system (Leconte et al. 2010). Therefore, the only assumption that may not be always verified is that the star spin has zero obliquity. Nevertheless, we make this assumption to keep our model as simple as possible, noting that there are some counterexamples to it (see Becker et al. 2020, for more details and possible interpretations).

The kinetic energy of the orbital motion of our system can be expressed as the energy of the relative motion of a body having the reduced mass of the system m = m_sm_p∕(m_s + m_p), where m_s is the mass of the star and m_p that of the planet, around the barycenter Z of the system. The kinetic energy of rotation of each body can be expressed in terms of the moment of inertia around their spin axes and the angular velocity of rotation around those axes. Specifically, we write: $$T=\frac{1}{2}m\left({\dot{r}}^2+r^2{\dot{f}}^2\right)+\frac{1}{2}{I}_{\text{s}}{\dot{\phi }}^2+\frac{1}{2}{I}_{\text{p}}{\dot{\psi }}^2,$$(12)

where r is the relative distance between the two bodies, f the true anomaly of their relative orbital motion, I_s the moment of inertia of the star around its rotation axis, φ the rotational coordinate of the star, I_p the moment of inertia of the planet around its rotation axis, and ψ the rotational coordinate of the planet; the dot over a variable indicates its time derivative. The angles f, φ, and ψ are illustrated in Fig. 2, where the barycenters of the star and the planet are O and O′, respectively, and the xy and x′ y′ planes coincide with the orbital plane. The Cartesian reference frames xy and x′y′ have their origins at O and O′, respectively,and their axes have fixed directions in an inertial space. Specifically, the angles φ and ψ are measured from the x and x′ axes to the principal axes of inertia corresponding to the quadrupole bulges of the star and the planet, respectively (see Sect. 2.1). This definition simplifies the expression of the gravitational potential energy Ψ as $$\begin{array}{lll}\Psi \hfill & =\hfill & -\frac{G{m}_{\text{s}}{m}_{\text{p}}}{r}-\frac{3G{m}_{\text{p}}}{2r^3}\left[Q_{\text{s}}+\frac{1}{2}{T}_{\text{s}}\cos 2(f-\phi )\right]\hfill \\ \hfill & \hfill & -\frac{3G{m}_{\text{s}}}{2r^3}\left[Q_{\text{p}}+\frac{1}{2}{T}_{\text{p}}\cos 2(f+\pi -\psi )\right],\hfill \end{array}$$(13)

where Q_s and T_s are the axisymmetric and non-axisymmetric quadrupole moments of the star as defined in Sect. 2.1, while Q_p and T_p are the same quantities for the planet (see below).

The potential energy of the system in Eq. (13) includes the first-order term of the gravitational energy that assumes that the masses of the two bodies are concentrated into their barycenters and the second-order term associated with their gravitational quadrupole moments with the stellar quadrupole acting on the planet and the planet quadrupole acting on the star, respectively. The quadrupole terms of the potential energy follow from the expression (6) of the quadrupole potential written for the case of the equatorial plane of the star and the planet (θ = π∕2) and considering that the angle α = f − φ in the case of the star or α = f + π − ψ in the case of the planet (cf. Fig. 2).

The bulges due to the equilibrium tides are almost aligned with the line joining the barycenters O and O′ of the two bodies and their magnitudes are almost constant for a nearly circular orbit. Therefore, they can be assumed to contribute to Q_{s, p}, respectively, while only the persistent non-axisymmetric quadrupole of magnetic origin is included into T_s and only the permanent deformation of the planet gives rise to T_p in our model (see the discussion in Sect. 2.2).

The equations of motion of our system can be derived from its Lagrangian as $$\frac{\text{d}}{\text{d}t}\left(\frac{\partial L}{\partial \dot{q}}\right)-\frac{\partial L}{\partial q}=0,$$(14)

where t is the time and q = r, f, φ, or ψ is any of the coordinates adopted to describe the motion of the system. In this way, we obtain the following equations of motion: $$\begin{array}{lll}\hfill & \hfill & m\ddot{r}-mr{\dot{f}}^2+\frac{G{m}_{\text{s}}{m}_{\text{p}}}{r^2}+\frac{9G}{2r^4}[Q_{\text{s}}{m}_{\text{p}}+\frac{1}{2}{T}_{\text{s}}{m}_{\text{p}}\cos 2(f-\phi )\hfill \\ \hfill & \hfill & +Q_{\text{p}}{m}_{\text{s}}+\frac{1}{2}{T}_{\text{p}}{m}_{\text{s}}\cos 2(f-\psi )]=0,\hfill \\ \hfill & \hfill & m\frac{\text{d}}{\text{d}t}\left(r^2\dot{f}\right)+\frac{3G}{2r^3}\left[T_{\text{s}}{m}_{\text{p}}\sin 2(f-\phi )+T_{\text{p}}{m}_{\text{s}}\sin 2(f-\psi )\right]=0,\hfill \\ \hfill & \hfill & \hfill \\ \hfill & \hfill & I_{\text{s}}\ddot{\phi }-\frac{3G}{2r^3}{T}_{\text{s}}{m}_{\text{p}}\sin 2(f-\phi )=0,\hfill \\ \hfill & \hfill & I_{\text{p}}\ddot{\psi }-\frac{3G}{2r^3}{T}_{\text{p}}{m}_{\text{s}}\sin 2(f-\psi )=0.\hfill \end{array}$$

Summing Eqs. (16)–(18) and integrating with respect to the time, we obtain an equation expressing the conservation of the total angular momentum J of the system: $$m{r}^2\dot{f}+I_{\text{s}}\dot{\phi }+I_{\text{p}}\dot{\psi }=J.$$(19)

The spin angular momentum of the planet is usually negligible because its rotation is synchronized by tides ($\dot{f}~\dot{\psi }$) and I_p ≪ mr² for m_p ≪ m_s. The angle f − ψ is always very small because it makes damped oscillations about the equilibrium zero value as we shall see from Eq. (59) below. Therefore, we can simplify Eq. (15) by considering that when the planet spin is synchronized and the orbit is close to circular, cos2(f − ψ) ≃ 1 to the first order, thus giving $$\ddot{r}-r{\dot{f}}^2+\frac{G{m}_{\text{T}}}{r^2}+\frac{9G{m}_{\text{T}}}{2m_{\text{s}}{r}^4}\left[{\tilde{Q}}_{\text{s}}+\frac{1}{2}{T}_{\text{s}}\cos 2(f-\phi )\right]=0,$$(20)

where m_T ≡ m_s + m_p is the total mass of the system and we have defined $${\tilde{Q}}_{\text{s}}\equiv Q_{\text{s}}+\left(\frac{m_{\text{s}}}{m_{\text{p}}}\right)\left(Q_{\text{p}}+\frac{1}{2}{T}_{\text{p}}\right).$$(21)

Fig. 2 Reference frames adopted to define the true anomaly f and the rotational coordinates of the star φ and the planet ψ to express the Lagrangian of our star-planet system (see text for details). The relative separation of the barycenters of the two bodies isr = OO′. The rotational coordinate of the star is defined as the angle between the x-axis, the directionof which is fixed in an inertial space, and the principal axis of inertia of the star with respect to which the angle α in Eq. (6) is measured. A similar definition is adopted for the rotational coordinate of the planet considering ψ as the angle between the axis x′ parallel to x and the principal axis of inertia of the planet passing through its permanent quadrupole bulge (see text).

2.4 Radial motion

The eccentricity of the orbit can be assumed to be small because it is rapidly damped by the tides inside the planet (see Sect. 2.5). Therefore, the terms involving the eccentricity will be of the first order as well as all the terms involving the quadrupoles, while their products will be of the second order, so we can neglect them. Since the angle f − ψ ~ 0, we neglect the contribution of the planet spin to the variation of the orbital angular momentum and re-write Eq. (16) as $$m\frac{\text{d}}{\text{d}t}\left(r^2\dot{f}\right)\simeq -\frac{3G}{2a^3}{T}_{\text{s}}{m}_{\text{p}}\sin [2(n-\Omega )t],$$(22)

where a is the orbit semimajor axis; n the orbital mean motion, n ≡ 2π∕P_orb, with P_orb the orbital period; and Ω ≡ 2π∕P_rot the spin angular velocity of the star with P_rot its rotation period. Equation (22) is a first-order equation valid for e ≪ 1 and can be immediately integrated with respect to the time to yield for n≠Ω $$\dot{f}=n{\left(\frac{a}{r}\right)}^2\left\{1+ϵ\cos \left[2(n-\Omega )t\right]\right\},$$(23)

where $$ϵ\equiv \frac{3}{4}\left(\frac{T_{\text{s}}}{m_{\text{s}}{a}^2}\right)\left(\frac{n}{n-\Omega }\right),$$(24)

is a very small quantity (ϵ ≪ 1) because thenon-axisymmetric component of the stellar quadrupole moment T_s is very small in comparison with the moment of inertia of the orbit ma² (see Sect. 3). When the orbit is nearly circular, we can define the small non-dimensional quantity x(t) such that |x(t)|≪ 1 from $$r\equiv a(1+x),$$(25)

and substitute it into Eq. (20) taking only the first-order terms. In this way, making use also of Eq. (23), we obtain an equation for x for n≠Ω: $$\ddot{x}+n^2\left\{1+\eta +\zeta \cos \left[2(n-\Omega )t\right]\right\}x=\xi +f(t),$$(26)

where η, ζ, and ξ are small quantities defined as $$\begin{array}{lll}\hfill & \hfill & \eta \equiv -18\frac{{\tilde{Q}}_{\text{s}}}{m_{\text{s}}{a}^2},\hfill \\ \hfill & \hfill & \zeta \equiv 9\left(\frac{T_{\text{s}}}{m_{\text{s}}{a}^2}\right)\left(\frac{2\Omega -n}{2n-2\Omega }\right),\text{and}\hfill \\ \hfill & \hfill & \xi \equiv -\frac{9}{2}{n}^2\left(\frac{{\tilde{Q}}_{\text{s}}}{m_{\text{s}}{a}^2}\right)=\frac{1}{4}{n}^2\eta ,\hfill \end{array}$$

while the forcing term $$f(t)\equiv A_{\text{x}}\left[\frac{T_{\text{s}}(t)}{I_{\text{s}}}\right]\cos \left[2(n-\Omega )t\right]$$(30)

has an amplitude A_x given by $$A_{\text{x}}\equiv \frac{3}{2}{n}^2\left(\frac{I_{\text{s}}}{m_{\text{s}}{a}^2}\right)\left(\frac{3\Omega -n}{2n-2\Omega }\right).$$(31)

In Eq. (30), we indicated explicitly the dependence of the stellar quadrupole term T_s on the time because of the finite lifetime of the magnetic flux tubes that are responsible for the term itself and normalized it to the stellar moment of inertia I_s instead of the moment of inertia m_sa² that explains the additional factor appearing in Eq. (31).

2.5 Tidal damping of the eccentricity of the planetary orbit

Equation (26) implies that x can make free oscillations along the orbit with the same frequency n of the orbital motion because η, ζ ≪ 1. Such oscillations correspond to the motion along an instantaneously eccentric orbit, therefore dissipative tides are excited within the planet and the star. They lead to the damping of the oscillations themselves with a timescale τ_e ≡−e∕(de∕dt) that is given by (cf. Jackson et al. 2008): $${\tau }_{\text{e}}^{-1}=\left[\frac{63R_{\text{p}}^5\sqrt{G{m}_{\text{s}}^3}}{4Q_{\text{p}}^{\prime \text{T}}{m}_{\text{p}}}+\frac{225R_{\text{s}}^5{m}_{\text{p}}\sqrt{G/m_{\text{s}}}}{16Q_{\text{s}}^{\prime \text{T}}}\right]a^{-13/2},$$(32)

where R_s and R_p are the radii of the star and the planet, m_s and m_p their masses, respectively, while $Q_{\text{s}}^{\prime \text{T}}$ and $Q_{\text{p}}^{\prime \text{T}}$ are the modified tidal quality factors of the star and the planet defined as $Q_{\text{s},\text{p}}^{\prime \text{T}}\equiv (3/2)(Q_{\text{s},\text{p}}^{\text{T}}/k_{2\text{s},\text{p}})$ (see Ogilvie 2014). The smaller the tidal quality factor, the larger the dissipation of the tidal kinetic energy and the damping of the eccentricity.

The tidal dissipation inside a telluric planet depends on the viscoelastic properties of its interior that in turn depend on its structure, composition, and the frequency of the tidal forcing. In our case, the tidal frequency corresponds to the orbital mean motion because we consider the eccentricity tide excited by the varying separation of the planet from its host star along its orbit. The simplest treatment, based on the Maxwell rheological model, severely underestimates the dissipation for tidal frequencies much higher than the inverse of the so-called Maxwell time, which is on the order of ~ 10⁴ yr for the Earth. Conversely, the Andrade rheological model gives a better description of the dissipation at the high frequencies characteristic of our case (Tobie et al. 2019; Bolmont et al. 2020). The multi-layered structure of a telluric planet also plays a relevant role, especially if it has an inner liquid core, a differentiated mantle, or a water ice layer, thus making simple models based on an homogeneous structure usually overestimating the dissipation (cf. Fig. 6 in Bolmont et al. 2020).

The values of $Q_{\text{p}}^{\text{T}}$ and the Love numbers k₂ computed by Tobie et al. (2019) for telluric planets consisting of a liquid metallic core and a silicate mantle, divided into an upper and a lower layer to account for a mineralogical transition occurring at pressures around 25 GPa, can be adopted for our investigation. Their model reproduces the observed tidal quality factors and Love numbers of the Earth at different tidal frequencies. Considering a forcing period of 1 day, their computed $Q_{\text{p}}^{\text{T}}$ and k₂ can be applied to the case of our telluric planets with orbital periods between 4 and 24 h because their dependencies on the tidal frequency in this range are smaller than the effect of the unknown interior structure and composition (e.g., Fig. 6 of Tobie et al. 2019). The value of $Q_{\text{p}}^{\prime \text{T}}=(3/2)(Q_{\text{p}}^{\text{T}}/k_2)$ increases from ~ 10³ for a planet with the mass of the Earth, up to ~2 × 10³ for a super-Earth of 8−10 Earth masses. Comparable results have been obtained by Bolmont et al. (2020)². Therefore, we can adopt the measured value of $Q_{\text{p}}^{\prime \text{T}}=1425$, relative to the semi-diurnal tide of our Earth with a forcing period of 12.4 h (Lainey 2016), as representative also of our telluric planets, although they are not expected to have surface oceans as our Earth. In other words, the above $Q_{\text{p}}^{\prime \text{T}}$ value is assumed to refer to the mantle only, while in the present Earth the $Q_{\text{p}}^{\prime \text{T}}$ value associated with the oceanic tides is smaller and the tidal dissipation is actually dominated by the friction applied by the topography on tidal gravito-inertial waves that propagate in the oceans (e.g., Egbert & Ray 2000, 2003; Tyler 2021).

The intense tidal heating experienced by a telluric planet as we shall find in our model, can deeply modify its rheological properties. For example, the melting of its silicate mantle can produce a decrease of its viscosity by several orders of magnitude that in turn can increase tidal dissipation (e.g., Walterová & Běhounková 2020). As we shall see, these complications as well as the ignorance of the interior structure and mineralogical composition of ultra-short period planets, will not affect our conclusions because the total dissipated energy inside the planet will not depend critically on the value of $Q_{\text{p}}^{\prime \text{T}}$. Therefore, we shall assume the above Earth value for reference in our investigation.

In our case, the tidal dissipation inside a telluric planet dominates that inside the star because $Q_{\text{p}}^{\prime \text{T}}\ll Q_{\text{s}}^{\prime \text{T}}$, given that $Q_{\text{s}}^{\prime \text{T}}\text{ ga}{10}^5-{10}^6$ as we shall see in Sect. 2.9. As a consequence, we shall assume that tidal dissipation is dominated by tides inside the planet and that all the dissipated mechanical energy is converted into heat inside the planet.

In conclusion, we can recast Eq. (26) including the damping effects of tides as $$\ddot{x}+2b\dot{x}+n^2\left\{1+\eta +\zeta \cos \left[2(n-\Omega )t\right]\right\}x=\xi +f(t),$$(33)

where $2b\equiv {\tau }_{\text{e}}^{-1}$. By multiplying both sides of Eq. (33) by ma²ẋ, we obtain $$\frac{\text{d}{E}_{\text{R}}}{\text{d}t}\simeq -2bm{a}^2{\dot{x}}^2+m{a}^2\dot{x}\left[f(t)+\xi \right],$$(34)

where $$E_{\text{R}}=\frac{1}{2}m{a}^2{\dot{x}}^2+m{a}^2{n}^2{x}^2$$(35)

is the total mechanical energy of the radial oscillations and we have neglected the small terms η, ζ ≪ 1 in the l.h.s.The first term in the r.h.s. of Eq. (34), that is, − 2bm a²ẋ², is the power dissipated by the tides inside the planet that produces a damping of the radial oscillations, while the term ma² ẋf(t) is the power that excites the radial oscillations. The small constant force ξ cannot excite any oscillation, but slightly shifts the position of equilibrium of the radial oscillator, so we shall neglect it in the following.

Equation (33) is a Mathieu equation with damping and a forcing term. A first-order parametric resonance cannot occur because it would imply that the frequency of variation of the pulsation should be equal to 2n which is not possible because it would imply Ω = 0. Parametric resonances of the second or higher orders are similarly not possible because of the dissipation term (Landau & Lifshitz 1960). Therefore, we shall consider the stationary solution of Eq. (33) under the effect of the forcing term in the r.h.s. For the sake of simplicity, we neglect the small terms η, ζ ≪ 1 in the l.h.s. of Eq. (33), thus assuming that the frequency of the free oscillations of x is equal to the unperturbed orbital mean motion n.

2.6 Solution of the radial motion equation and power dissipated inside the planet

The quadrupole moment T_s is not constant because the magnetic flux tubes in the convection zone of the host star have a finite lifetime (see Sect. 2.2). We quantify this effect by considering the autocorrelation function R_f (τ) of the forcing term on the r.h.s. of Eq. (33), $$R_f(\tau )=\langle f(t+\tau )f(t)\rangle ,$$(36)

where τ is the time lagand the angular brackets indicate an ensemble average over many realizations of the forcing term f(t). It oscillates with the frequency γ ≡ 2(n − Ω), that is largein comparison with the inverse of the mean lifetime of the magnetic flux tubes producing the quadrupole term T_s (t) in Eq. (30). Therefore, R_f(τ) can be written as the product of the autocorrelation of the quadrupole term T_s by that of the rapidly oscillating term cosγt. The autocorrelation of cosγt is (1∕2)cosγτ, while that of T_s is assumed to be proportional to a decaying exponential with a typical timescale ${(D/2)}^{-1}$, that we take equal to the mean lifetime of the magnetic flux tubes. By combining these two contributions we obtain $$R_f(\tau )=\langle f(t+\tau )f(t)\rangle =\frac{1}{2}{A}_{\text{s}}^2\exp \left(-\frac{D}{2}\left|\tau \right|\right)\cos (\gamma \tau ).$$(37)

where the amplitude A_s is $$A_{\text{s}}={\left(\frac{\langle T_{\text{s}}^2\rangle }{I_{\text{s}}^2}\right)}^{1/2}{A}_{\text{x}},$$(38)

with A_x given by Eq. (31). The quadratic mean relative amplitude of the stellar quadrupole $\langle T_{\text{s}}^2\rangle /I_{\text{s}}^2$, can be obtained from the mean filling factor f_s and intensity of the magnetic field B₀ in the overshoot layer as discussed in Sect. 2.1. An estimate of the timescale ${(D/2)}^{-1}$ can be obtained by the rate of decrease of the amplitudes of successive peaks in the autocorrelation function of the light curve of an active star, while the lag between consecutive peaks gives a measure of the stellar rotation period (cf. Lanza et al. 2014; Giles et al. 2017).

The solution of Eq. (33) can be obtained by considering the Fourier transforms of the functions x(t) and f(t). We define the Fourier transform $\widehat{x}(\nu )$ of x(t) as $$\widehat{x}(\nu )=\frac{1}{\sqrt{2\pi }}{\displaystyle {\int }_{-\infty }^{\infty }x}(t)\exp (-i\nu t)\text{d}t,$$(39)

where ν is the frequency and $i=\sqrt{-1}$ the imaginary unit. The inverse Fourier transform allows us to express x(t) in terms of $\widehat{x}(\nu )$ as $$x(t)=\frac{1}{\sqrt{2\pi }}{\displaystyle {\int }_{-\infty }^{\infty }\widehat{x}}(\nu )\exp (i\nu t)\text{d}\nu .$$(40)

The power spectrum of x can be defined as $S_x(\nu )=\widehat{x}(\nu ){\widehat{x}}^{*}(\nu )$, where the asterisk indicates complex conjugation. It is equal to the Fourier transform of the autocorrelation function of x, that is, R_x (τ) ≡⟨x(t + τ)x(t)⟩ $$S_x(\nu )=\frac{1}{\sqrt{2\pi }}{\displaystyle {\int }_{-\infty }^{\infty }{R}_x}(\tau )\exp (-i\nu \tau )\text{d}\tau .$$(41)

The expectation value of the square of x, that is, $\langle x{(t)}^2\rangle$, when the mean of x is zero (⟨x⟩ = 0), can be obtained by applying the inverse Fourier transform for τ = 0 because by definition $\langle x{(t)}^2\rangle =R_x(0)$ $$\langle x{(t)}^2\rangle =\frac{1}{\sqrt{2\pi }}{\displaystyle {\int }_{-\infty }^{\infty }{S}_x}(\nu )\text{d}\nu .$$(42)

By taking the Fourier transforms of both sides of Eq. (33), neglecting the term ξ and the small terms η, ζ ≪ 1 because parametric oscillations are not excited, and finally taking the complex conjugate of the transforms themselves, we find the power spectrum of x as $$S_x(\nu )\equiv \widehat{x}(\nu ){\widehat{x}}^{*}(\nu )=\frac{S_f(\nu )}{{(n^2-{\nu }^2)}^2+4b^2{\nu }^2},$$(43)

where $S_f(\nu )\equiv \widehat{f}(\nu ){\widehat{f}}^{*}(\nu )$ is the power spectrum of the forcing f(t).

The power spectrum of f(t) is the Fourier transform of its autocorrelation function as given by Eq. (37). After application of the definition (Eq. (39)) and some algebraic rearrangements, we obtain (cf. the analogous system studied by Zhou et al. 1997) $$S_f(\nu ;\gamma ,D)=\frac{A_{\text{s}}^{2}D}{2\sqrt{2\pi }}\frac{\frac{1}{4}{D}^2+{\gamma }^2+{\nu }^2}{{\left(\frac{1}{4}{D}^2+{\gamma }^2-{\nu }^2\right)}^2+D^2{\nu }^2},$$(44)

where we have explicitly indicated that S_f(ν) depends also on the parameters γ and D that enter into the expression of the autocorrelation of the forcing function f (cf. Eq. (37)).

The mean power P_T dissipated by the tides inside the planet is given by (cf. Eq. (34)) $$P_{\text{T}}\equiv 2bm{a}^2\langle \dot{x}{(t)}^2\rangle ,$$(45)

where ⟨ẋ²(t)⟩ can be obtained from its power spectrum (cf. Eq. (42)) $$\langle {\dot{x}}^2\rangle =\frac{1}{\sqrt{2\pi }}{\displaystyle {\int }_{-\infty }^{\infty }{S}_{\dot{x}}}(\nu )\text{d}\nu =\frac{1}{\sqrt{2\pi }}{\displaystyle {\int }_{-\infty }^{\infty }{\nu }^2}{S}_x(\nu )\text{d}\nu ,$$(46)

because the Fourier transform of the time derivative is $\widehat{\dot{x}}=i\nu \widehat{x}$ that implies S_ẋ(ν) = ν²S_x(ν). By substituting Eq. (43), we obtain $$\begin{array}{lll}{P}_{\text{T}}\hfill & =\hfill & 2bm{a}^2\frac{1}{\sqrt{2\pi }}{\displaystyle {\int }_{-\infty }^{\infty }\frac{{\nu }^2{S}_f(\nu ;\gamma ,D)}{{(n^2-{\nu }^2)}^2+4b^2{\nu }^2}}\text{d}\nu \hfill \\ \hfill & \simeq \hfill & 2bm{a}^2\frac{1}{\sqrt{2\pi }}{S}_f(n;\gamma ,D){\displaystyle {\int }_{-\infty }^{\infty }\frac{{\nu }^2\text{d}\nu }{{(n^2-{\nu }^2)}^2+4b^2{\nu }^2}}\hfill \\ \hfill & =\hfill & 2bm{a}^2\frac{1}{\sqrt{2\pi }}\frac{\pi }{2b}{S}_f(n;\gamma ,D)=\sqrt{\frac{\pi }{2}}m{a}^2{S}_f(n;\gamma ,D),\hfill \end{array}$$

where we made use of the fact that the denominator in the integrand of Eq. (47) becomes very small when ν is very close to n because b ≪ n. Therefore, we can take the factor S_f(ν;γ, D) out of the integral by computing it for ν = n and integrate the remaining factor by applying the theorem of the residues (see Appendix A) to obtain Eq. (49). Equation (49) shows that the mean dissipated power P_T is independent of the dissipation rate b, which is a very useful result in view of our ignorance of the internal structure and rheology of the planet determining its modified tidal qualityfactor (cf. Sect. 2.5).

The value of S_f(n;γ, D) is given by Eq. (44) with ν = n. It reaches its maximum when $\gamma \equiv {\gamma }_{\text{res}}=\sqrt{n^2-D^2/4}$, that is, 2(n − Ω) ≃ n or Ω ≃ n∕2 because D ≪ n in very active stars in which active longitudes persist for timescales ranging from hundreds of days up to a decade and the orbital period of the planet is less than a few days. Therefore, the maximum of the dissipated power at the resonance is sharp and much greater thanthe power dissipated far from the resonance. By combining Eqs. (44) and (49), the maximum dissipated power is $$P_{\text{T}\max }=\frac{1}{2}\frac{m{a}^2{A}_{\text{s}}^2}{D},$$(50)

where we made use of the fact that D ≪ n. The power decreases rapidly from its maximum when γ varies around the resonance value γ_res. A power of one half of P_{T max} is obtained for γ ~ γ_res[1 ± D∕(2n)], considering Eq. (44) and that D ≪ n.

The angular velocity Ω will change during the evolution of the star in its pre-main-sequence and main-sequence phases owing to the changes in the stellar moment of inertia and the loss of angular momentum carried away by the stellar magnetized wind. Therefore, the value of P_T will vary as Ω evolves and crosses the resonance at Ω ≃ n∕2. As we shall see in Sect. 3, the crossing of the resonance occurs over a time interval much shorter than the internal timescale of heat redistribution inside the planet, so that all the power dissipated during the resonance is stored into the planet and released on a much longer timescale. Since the resonance is sharp because D ≪ n (see above), the total energy E_D dissipated inside the planet during the crossing of the resonance can be evaluated as $$E_{\text{D}}\equiv {\displaystyle {\int }_{-\infty }^{\infty }{P}_{\text{T}}}(t)\text{d}t\simeq \sqrt{\frac{\pi }{2}}m{a}^2{\left(\frac{\text{d}t}{\text{d}\gamma }\right)}_{\text{res}}{\displaystyle {\int }_{-\infty }^{\infty }{S}_f}(n;\gamma ,D)\text{d}\gamma ,$$(51)

where dt∕dγ is the inverse of the rate of change of the frequency of the forcing term f that is produced by the evolution of the stellar rotation and of the orbit semimajor axis under the action of the tides inside the star (cf. Sect. 2.9). Such an inverse rate of change is evaluated at the resonance, that is, for $\gamma =\sqrt{n^2-D^2/4}$, and it is taken out of the integral thanks to the sharpness of the resonance. The integral can be evaluated with the theorem of the residues (see Appendix A, Eq. (A.15)) giving $$E_{\text{D}}=\frac{\pi }{2}m{a}^2{A}_{\text{s}}^2{\left(\frac{\text{d}t}{\text{d}\gamma }\right)}_{\text{res}}.$$(52)

It is remarkable that E_D is independent of D and is a function only of A_s, that is, of the non-axisymmetric component of the stellar quadrupole moment T_s, the orbital parameters, and the inverse rate of change of γ at the resonance, that is, of the evolution of the angular momentum of the system. The amount of dissipated energy decreases rapidly with the increase of the orbit semimajor axis, scaling as E_D ∝ a⁻⁸ (cf. Eqs. (31), (38), and (52)), which implies that this heating mechanism is mostly limited to operate in planets with ultra-short-period orbits (P_orb≲1 days) around young active late-type stars.

2.7 Libration of the planet

Now we focus on the Eqs. (16)–(18) for the angular variables and consider, for simplicity, the case of a circular orbit. As we shall see, such an assumption is justified because the libration of the planet that we want to investigate occurs far from the resonance that excites the radial motion. In such a regime, x(t) is damped on a timescale much shorter than the lifetime of the system (see Sect. 3) justifying the assumption. For a general discussion of the effects of librations in the case of telluric planets, see, for example, Le Bars (2016) and references therein.

Indicating with a the constant orbital radius r, Eqs. (16)–(18) read $$\begin{array}{l}\ddot{f}+\frac{3}{2}{n}^2\left[\left(\frac{T_{\text{s}}}{m_{\text{s}}{a}^2}\right)\sin 2\alpha +\left(\frac{T_{\text{p}}}{m_{\text{p}}{a}^2}\right)\sin 2\beta \right]=0,\\ \ddot{\phi }-\frac{3}{2}{n}^2\left(\frac{m_{\text{p}}}{m_{\text{T}}}\right)\left(\frac{T_{\text{s}}}{I_{\text{s}}}\right)\sin 2\alpha =0,\\ \ddot{\psi }-\frac{3}{2}{n}^2\left(\frac{m_{\text{s}}}{m_{\text{T}}}\right)\left(\frac{T_{\text{p}}}{I_{\text{p}}}\right)\sin 2\beta =0,\end{array}$$

where we made use of Kepler III law (n²a³ = Gm_T) and introduced the angle α ≡ f − φ = (n − Ω)t and the libration angle of the planet β = f − ψ. By subtracting Eqs. (55) from (53), we obtain an equation for β as $$\ddot{\beta }+\frac{1}{2}{\omega }_{\text{p}}^2\sin 2\beta =-\frac{3}{2}{n}^2\left(\frac{T_{\text{s}}}{m_{\text{s}}{a}^2}\right)\sin \gamma t,$$(56)

where γ = 2(n − Ω) is the frequency of the forcing introduced in Sect. 2.6. This is the equation of a forced simple pendulum with pulsation ω_p given by $${\omega }_{\text{p}}^2=3n^2\left[\frac{1}{m_{\text{p}}{a}^2}+\left(\frac{m_{\text{s}}}{m_{\text{T}}}\right)\frac{1}{I_{\text{p}}}\right]T_{\text{p}}\simeq 3n^2\left(\frac{T_{\text{p}}}{I_{\text{p}}}\right),$$(57)

because the moment of inertia of the planet I_p is much smaller than m_pa² and m_s ≃ m_T.

The effect of tides is not included into Eq. (56) as it tends to synchronize the rotation of the planet with the orbital motion. In other words, tides produce a damping of $\dot{\beta }$ over a timescale τ_rot given by (cf.Eq. (9) in Gu et al. 2003) $${\tau }_{\text{rot}}^{-1}=\left|\frac{\ddot{\psi }}{n-\dot{\psi }}\right|=\left|\frac{\ddot{\beta }}{\dot{\beta }}\right|=\frac{9n}{2h_{\text{p}}{Q}_{\text{p}}^{\prime \text{T}}}\left(\frac{m_{\text{s}}}{m_{\text{p}}}\right){\left(\frac{R_{\text{p}}}{a}\right)}^3,$$(58)

where we considered that $\ddot{\beta }=-\ddot{\psi }$ because theorbital motion is not appreciably perturbed given that I_p ≪ ma² and defined $h_{\text{p}}\equiv I_{\text{p}}/(m_{\text{p}}{R}_{\text{p}}^2)$, while the other symbols have already been introduced above. For the Earth, h_p = 0.3307.

The equation for the libration angle β, including tidal effects and in the limit of a small amplitude of the libration, that is, β ≪ 1, becomes $$\ddot{\beta }+2b_{\mathcal{l}}\dot{\beta }+{\omega }_{\text{p}}^2\beta =-\frac{3}{2}{n}^2\left(\frac{T_{\text{s}}}{m_{\text{s}}{a}^2}\right)\sin \gamma t,$$(59)

where $2b_{\mathcal{l}}\equiv {\tau }_{\text{rot}}^{-1}$ is the damping constant.

Equation (59) shows that the libration of the planet can be resonantly excited when γ = 2(n − Ω) is close to ω_p with the width of the resonance controlled by the tidal dissipation through the parameter b_ℓ. Considering Eq. (57) and that T_p ~ (10⁻³−10⁻⁴)I_p (cf. Eq. (10)), we have that ω_p ≪ n, thus the resonance occurs only when Ω ~ n, that is, the angular velocity of rotation of the star is close to the mean orbital motion of the planet. On the other hand, we saw in Sect. 2.4 that the radial motion is excited when Ω ~ n∕2. Therefore, the libration resonance is well separated from the resonance that excites the radial motion provided that the damping parameters b and b_ℓ are much smaller than n and ω_p, respectively.

2.8 Power dissipated by librations inside the planet

The power P_libr dissipated by the tides inside a librating planet can be written as (Ogilvie 2014, Sect. 2.2) $$P_{\text{libr}}=\widehat{\omega }{T}_{\text{t}},$$(60)

where $\widehat{\omega }=2(\dot{f}-\dot{\psi })=2\dot{\beta }$ is the frequency of the semidiurnal tide acting on the planet and T_t is the tidal torque acting on the planet spin that we write as $T_{\text{t}}=-2b_{\mathcal{l}}{I}_{\text{p}}\dot{\beta }$. Therefore, the mean dissipated power becomes $$P_{\text{libr}}=4b_{\mathcal{l}}{I}_{\text{p}}\langle {\dot{\beta }}^2\rangle .$$(61)

Following the same procedure illustrated in Sect. 2.6, we compute the power spectrum of $\dot{\beta }$ from the power spectrum of the forcing term and assume that b_ℓ ≪ ω_p to find $$P_{\text{libr}}=\sqrt{2\pi }{I}_{\text{p}}{S}_{f_{\beta }}({\omega }_{\text{p}};\gamma ,D),$$(62)

where the power spectrum of the forcing term $$f_{\beta }\equiv -\frac{3}{2}{n}^2\left[\frac{T_{\text{s}}(t)}{m_{\text{s}}{a}^2}\right]\sin \gamma t,$$(63)

is similar to that of f in Eq. (44), the only difference between f and f_β being the factor (3Ω − n)∕γ appearing in the amplitude of f(t) that changes the amplitude of the power spectrum accordingly (cf. Eqs. (30) and (31)). In Eq. (63), we made the dependence of T_s on the time explicit to allow an easier comparison with Eq. (30).

The total energy E_β dissipated in the planetwhen the frequency γ passes across the resonance that excites the librations is obtained by integrating the power P_libr over the time. Following the same approach as in Sect. 2.6, we find $$E_{\beta }=\pi I_{\text{p}}{A}_{\beta }^2{\left(\frac{\text{d}t}{\text{d}\gamma }\right)}_{\beta \text{res}},$$(64)

where A_β is the amplitude of the power spectrum in Eq. (44) in the case of the forcing (63) and the derivative d t∕dγ is computed at the resonance, that is, for $\gamma ={\gamma }_{\beta \text{res}}\equiv \sqrt{{\omega }_{\text{p}}^2-D^2/4}$.

By comparing Eqs. (52) and (64) and noting that the amplitudes A_s and A_β as well as the derivatives (dt∕dγ) at the respective resonances are comparable with each other, we find that the ratio of the total dissipated energies is E_β ∕E_D ~ 2I_p∕(ma²) ~ 10⁻⁵−10⁻⁴ ≪ 1. Therefore, we shall neglect the energy dissipated by the librations inside the planet in comparison with the energy dissipated by the radial oscillations discussed in Sect. 2.6.

2.9 Long-term evolution of the stellar rotation and orbital semimajor axis

The quadrupole moments of the star and the planet produce a spin-orbit coupling that exchanges angular momentum between the orbital motion and the rotation of the bodies. In the most extreme cases, this effect produces a modulation of the orbital period with a relative amplitude ΔP_orb∕P_orb≲10⁻⁶. These angular momentum exchanges are cyclic, therefore, they neither affect the long-term evolution of the orbit nor of the stellar rotation. More details can be found in Lanza (2020a,b), but they do not contribute a significant tidal dissipation inside a solid planet, so we shall neglect them.

The long-term evolution of the system is driven by the loss of angular momentum in the stellar magnetized wind and by the star-planet tidal interactions. Here the angular momentum loss is parametrized following the simple model by Bouvier et al. (1997) and assuming that the star is rotating as a rigid body after settling on the zero-age main sequence (ZAMS). The internal structure of the star is assumed to be steady, thus its moment of inertia I_s is constant during the evolution of its rotation. Such approximations are useful for our illustrative applications in Sect. 3 because we shall focus on the first ~100 Myr of evolution of the stellar rotation after the ZAMS when they are verified rather well. More sophisticated models can be applied for a detailed comparison with the observations (e.g., Amard et al. 2016; Gallet et al. 2018; Benbakoura et al. 2019), but this is postponed to later works because here we focus on the physical mechanism producing the internal heating of the planet.

The angular momentum loss rate dL_w∕dt due to the magnetized wind is described by $$\frac{\text{d}{L}_{\text{w}}}{\text{d}t}=\{\begin{array}{r}\hfill -K_{\text{w}}{\Omega }^3{\left(\frac{R_{\text{s}}}{R_{\odot }}\right)}^{1/2}{\left(\frac{m_{\text{s}}}{M_{\odot }}\right)}^{-1/2}\text{for}\Omega <{\Omega }_{\text{sat}}(m_{\text{s}}),\\ \hfill -K_{\text{w}}{\Omega }_{\text{sat}}^2\Omega {\left(\frac{R_{\text{s}}}{R_{\odot }}\right)}^{1/2}{\left(\frac{m_{\text{s}}}{M_{\odot }}\right)}^{-1/2}\text{for}\Omega \ge {\Omega }_{\text{sat}}(m_{\text{s}}),\end{array}$$(65)

where K_w = 2.7 × 10⁴⁰ kg m² s⁻¹, and Ω_sat is the angular velocity corresponding to the saturation of the wind angular momentum loss rate that is a function of the star mass m_s. We shall assume Ω_sat = 14, 8, 3 Ω_⊙ for a stellar mass m_s = 1.0, 0.8, 0.5 M_⊙, respectively, and linearly interpolate for masses in between.

The tidal torque acting on the stellar rotation Γ_t can be expressed in terms of the parameters of the system following, e.g., Mardling & Lin (2002). After some re-arrangement of the terms and taking into account also the torque produced by the stellar wind, we can write the equation for the evolution of the stellar rotation as $$I_{\text{s}}\frac{\text{d}\Omega }{\text{d}t}={\Gamma }_t+\frac{\text{d}{L}_{\text{w}}}{\text{d}t}=9\frac{G{m}_{\text{p}}^2}{a}{\left(\frac{R_{\text{s}}}{a}\right)}^5\frac{1}{Q_{\text{s}}^{\prime \text{T}}}\left(\frac{n-\Omega }{n}\right)+\frac{\text{d}{L}_{\text{w}}}{\text{d}t},$$(66)

where all the symbols have already been defined.

The stellar modified tidal quality factor $Q_{\text{s}}^{\prime \text{T}}$ is remarkably uncertain because it depends on the processes that dissipate the kinetic energy of the tides raised inside the star by the orbiting planet.

Traditionally, the tidal response of a star is decomposed into a large-scale nearly-hydrostatic deformation, called the equilibrium tide, and the various kinds of waves that are excited by the time-dependent tidal potential in a reference frame rotating with the star, the so-called dynamical tide. In the case of late-type stars with close-by planets, the dynamical tide consists of the internal gravity waves excited inside their radiative zones and the inertial waves excited into their convection zones when the tidal frequency $\left|\widehat{\omega }\right|\equiv |2(n-\Omega )|\le 2\Omega$ (see Barker 2020). The excitation and dissipation of dynamical tides are associated with a remarkably larger tidal dissipation than in the case of the equilibrium tide and correspond to a decrease of the modified tidal quality factor $Q_{\text{s}}^{\prime \text{T}}$ by 2− 4 orders of magnitude.

Internal gravity waves become relevant when their amplitude increases beyond the threshold required for their efficient dissipation inside the radiative zone of the star, that requires a planet mass remarkably greater than the mass of Jupiter in the case of late-type stars of mass 0.7−1.1 solar masses, ages younger than a few Gyr, and orbital periods on the order of 1 day or shorter as in our case; see, e.g., Sect. 6 of Barker (2020) or Ahuir et al. (2021). Therefore, their dissipation does not appreciably contribute to the evolution of the systems we consider here.

On the other hand, inertial waves in the convection zone, the restoring force of which is the Coriolis force, play a relevant role in our case in the range of tidal frequencies where they are excited. The associated $Q_{\text{s}}^{\prime \text{T}}$ shows a highly erratic dependence on the tidal frequency for given stellar parameters, therefore, it is useful to adopt an average value to compute theevolution of star-planet systems as in, e.g., Mathis (2015) or Gallet et al. (2017). Such an average value is generally 2− 3 orders of magnitude smaller than that associated with the equilibrium tide and is proportional to the inverse of the square of the stellar angular velocity because the Coriolis acceleration is proportional to Ω, viz., $Q_{\text{s}}^{\prime \text{T}}\propto {\Omega }^{-2}$ (Ogilvie & Lin 2007; Ogilvie 2013; Barker 2020).

In view of these results of the tidal theory, we assume $$Q_{\text{s}}^{\prime \text{T}}=Q_{\text{s0}}^{\prime \text{T}}{({\Omega }_0/{\Omega }_{\text{ref}})}^{-2}$$(67)

in the regime of tidal frequency where inertial waves are excited with $Q_{\text{s0}}^{\prime \text{T}}=4\times {10}^7$ for m_s ≥ 0.9 M_⊙ or $Q_{\text{s0}}^{\prime \text{T}}=2\times {10}^7$ for m_s < 0.9 M_⊙, Ω₀ being the angular velocity of the star on its arrival on the ZAMS, and Ω_ref the angular velocity when the modified tidal quality factor is equal to $Q_{\text{s0}}^{\prime \text{T}}$; we assume it to correspond to a rotation period of 10 days, viz. Ω_ref = 7.27 × 10⁻⁶ s⁻¹. The values of $Q_{\text{s}}^{\prime \text{T}}$ given by Eq. (67) are similar within a factor of ≈ 1−3 to those computed by Mathis (2015) and Gallet et al. (2017) with an approximate two-layer model for late-type stars on the main sequence, while they are about a factor of 2−5 greater than the values computed by Barker (2020) with stratified models. In view of the significant uncertainties in the present tidal theory, a closer agreement between their estimates and ours is not pursued here.

When inertial waves are not excited, we increase the modified tidal quality factor $Q_{\text{s}}^{\prime \text{T}}$ by a factor on the order of 10³, thus producing a much smaller tidal coupling between the stellar spin and the orbit (cf. Ogilvie & Lin 2007; Barker 2020). More precisely, given the present uncertainties of the equilibrium tide theory to be applied in that regime (see, e.g., Sect. 2 of Barker 2020), we take the liberty of adjusting the precise value of $Q_{\text{s}}^{\prime \text{T}}$ to reproduce the current orbital period and rotation period of the star in the modeled systems starting from their adopted initial conditions (see Sect. 3 for details).

The evolution of the orbit semimajor axis a can be computed from the tidal torque Γ_t considering the conservation of the total angular momentum of the system. After some algebra, we find $$\frac{\text{d}a}{\text{d}t}=\frac{18}{Q_{\text{s}}^{\prime \text{T}}}{\left(\frac{R_{\text{s}}}{a}\right)}^5\left(\frac{m_{\text{p}}}{m_{\text{s}}}\right)\left(\frac{\Omega -n}{n}\right)na.$$(68)

In the next section, we shall make use of Eqs. (66) and (68) to compute some illustrative models of the evolution of the stellar rotation and the semimajor axis for some planetary systems in order to show how the conditions for the excitation of the radial oscillations of the planet can be realized along their evolution.

We shall focus on the dissipation inside the planet considering the evolution from the ZAMS. The angular velocity of the star on the ZAMS depends on its past evolution, in particular on its contraction during the pre-main-sequence phase, and the particular moment when the mechanical coupling with its circumstellar disk came to an end leaving the contracting protostar free to spin up. Therefore, the initial rotation period 2π∕Ω₀ on the ZAMS is strongly dependent on the protostar mass and the time when the decoupling from the disk occurred and can range from a few hours up to several days in the case of a long-lasting disk locking (e.g., Bouvier et al. 1997; Gallet & Bouvier 2015). In view of our ignorance on the pre-main-sequence evolution of our stars, we shall regard Ω₀ as a free parameter and adjust it to obtain the resonant excitation of the radial motion at some time after the ZAMS.

3 Applications

We consider three MESA stellar models (see Fields et al. 2018; Paxton et al. 2019, and references therein) to evaluate the non-axisymmetric quadrupole moment for three late-type active stars with masses typical of those hosting ultra-short-period planets. The non-axisymmetric quadrupole moment T_s is computed according to Eq. (8) considering θ₀ = 30° that corresponds to a filling factor f_s ≃ 0.07 of the whole stellar surface. The corresponding relative amplitude of the optical light modulation of the star is ≈ 25 percent that is typical of late-type stars with a rotation period close to ≈ 1 day and younger than 100 Myr as we shall consider in our applications (e.g., Messina 2021).

The integral $J$ in Eq. (9) is evaluated by considering the internal structure model of the host star as computed by means of the MESA Web interface³. We assume a metallicity Z = 0.02, a ratio of the mixing length to the pressure scale height equal to 2.0, and standard values for the other parameters and options. We list the stellar parameters in the cases of models of mass 0.7, 0.8 and 0.9 M_⊙ in Table 1 at the indicated ages that correspond to models during the main-sequence phase of their evolution. Given their slow evolution along the main sequence, these models can be regarded as representative of the stellar internal structures all along this phase for our purposes.

The magnetic field B₀ in the overshoot layer immediately below the base of the convection zone at r = r_b is chosen to have T∕I_s = 10⁻⁶, where I_s is the stellar moment of inertia. Those values of B₀ are comparable with the magnetic field instability thresholds for the emergence of the flux tubes stored in the overshoot layers of late-type, rapidly rotating stars according to the model of Granzer et al. (2000). The upper limit r_L to which the integration in Eq. (9) is extended is also listed in Table 1 as well as the magnetic energy E_mag of the vertical flux tube forf_s = 0.07 according tothe model of Lanza (2020a).

Three illustrative applications to ultra-short-period planets are introduced in this section, specifically to CoRoT-7b, Kepler-78b, and K2-141b, respectively. They are not intended as detailed and accurate models of the evolution of these systems, but only as examples to show the typical powers, total energies, and timescales involved in the heating mechanism introduced in this paper. Therefore, we shall not compare the results of our models with the observations in detail. In the case of CoRoT-7 and K2-141, another planet has been detected in the system with an orbital period longer than that of the considered innermost planet.

The presence of outer companions may be a common property of ultra-short-period planets and could have played a decisive role in their migration close to the host star (cf. Dai et al. 2018; Petrovich et al. 2019). For simplicity, we shall not consider the effects due to such companions and model the evolution assuming that these systems consist only of the host star and their closest planet.

The parameters assumed for the considered systems are listed in Table 2. To model their evolution starting from the ZAMS corresponding to the time t = 0, we adopt initial values of the orbit semimajor axis a(t = 0), orbital period P_orb(t = 0), and stellar rotation period P_rot(t = 0) = 2π∕Ω₀ as listed in Table 3. We start from rapidly rotating stars on the ZAMS and do not consider the evolution of the radius of the host star along its main-sequence lifetime. On the other hand, had the star settled on the main sequence as a slow rotator with an initial P_rot > 2P_orb, our mechanism could not operate because the spin-orbit resonance responsible for the excitation of the radial oscillations of the planetary orbit and the consequent dissipation could not occur.

Another assumption of our model is that the planet was on a very close orbit since the star started its main-sequence evolution as could happen in the case of a migration occurred inside the protoplanetary disk or an early inward scattering in multiplanet systems. Other mechanisms that could bring the planet very close to its host star after a much longer time should not produce a heating of the planet by our mechanism because the rotation of the star would have already been braked below the value required to excite the radial motion resonance (e.g., Lanza & Shkolnik 2014).

The angular velocity Ω_sat corresponding to the saturation of the angular momentum loss in the stellar wind is computed by a linear interpolation considering its dependence on the mass of the host star (cf. Sect. 2.9). We consider a fixed amplitude of the non-axisymmetric quadrupole moment T_s∕I_s = 10⁻⁶ because it is relevant only close to the resonance that excites the radial orbital oscillations of the planet. The lifetime ${(D/2)}^{-1}$ of the magnetic flux tubes responsible for the non-axisymmetric quadrupole is assumed to be 500 days in all the cases (cf. Sect. 2.2).

The stellar modified tidal quality factors are listed in Table 3 for the case when inertial waves are excited (IW) and for the case when only the equilibrium tide (ET) produces dissipation inside the star as discussed in Sect. 2.9. The value of $Q_{\text{s}}^{\prime \text{T}}$ in the case of the equilibrium tide is adjusted to reproduce both the observed current orbital period and stellar rotation period as indicated by the vertical dotted lines in Figs. 3, 5, and 7. The system ages estimated with our model at those epochs are listed in Table 3. They can differ from those in the literature as reported in Table 2 because of the large uncertainty (up to 30− 40 percent) in those estimates, mainly based on gyrochronology, as well as the simplifying assumptions adopted in our illustrative models.

For Kepler-78 and K2-141, the initial rotation period is shorter than the orbital period which leads to an increase of the orbital period and semimajor axis because of the angular momentum transferred by tides from the stellar rotation to the orbital motion during the initial phase of their evolution (cf. Figs. 6 and 8). Conversely, in CoRoT-7, P_rot is initially longer than P_orb, therefore tides start immediately to decrease the semimajor axis and the orbital period of the planet (cf. Fig. 4).

A slope change occurs in the plots of the orbital period vs. the time when the stellar rotation period is twice the orbital period. This corresponds to the transition from the strong tidal coupling maintained by the excitation of inertial waves inside the star to the regime with only the equilibrium tide. This remarkably reduces the tidal torque inside the star and leads to a much slower decay of the planetary orbit. The excitation of the radial orbital oscillations of the planet occurs close to the transition between these two tidal regimes because the resonance frequency ${\gamma }_{\text{res}}=\sqrt{n^2-D^2/4}$ is very close to the orbital mean motion n, given that D ≪ n.

The relative amplitude of the radial oscillations at resonance ${\langle x^2\rangle }^{1/2}$ is also reported in Table 3. It is computed from the power spectrum of the forcing as illustrated in Sect. 2.6 (cf. Eq. (42)). Its value is very small, thus verifying our assumption that |x| ≪ 1 in Eq. (25) in Sect. 2.4, and is extremely difficult to detect because of the short duration of the excitation in comparison to the main-sequence lifetime of the systems (see below). Any residual eccentricity would be damped by the tides inside the planet on a timescale τ_e = 1∕(2b) that is shorter than or comparable with a Myr when considering a value of $Q_{\text{p}}^{\prime \text{T}}=1425$ as in the case of our Earth (Lainey 2016).

The maximum dissipated power P_{T max} is listed in Table 3 together with the total energy dissipated inside the planet E_D and the time needed to cross the resonance, that is, the time interval during which the instantaneously dissipated power is greater than or equal to one half of the maximum power. The resonance crossing time is very short in comparison with the evolutionary timescales and the timescale needed by sub-solidus convection to transport the heat dissipated inside the mantle of a telluric planet, that is, on the order of 10⁸ yr in the case of our Earth. For this reason, we can assume that all the heat dissipated by the proposed mechanism is stored into the planetary interior. The amount of heat available per unit mass of the planet E_D ∕m_p is also listed in Table 3 and it shall be used to infer the effects of the heating in the next section.

The heat dissipation produced by our mechanism is restricted to a short time interval on the order of 10⁴ yr occurring after ≈40−50 Myr since our host stars have settled on the main-sequence. Another resonance could have occurred during the pre-main-sequence evolution of those systems when the contracting and spinning up protostar could have matched the condition for the resonance P_rot ≃ 2P_orb while evolving toward the ZAMS. For simplicity, we do not consider this possibility and limit ourselves to model only the resonance occurring after the ZAMS. Finally, we note that when P_rot = 3P_orb, the amplitude of the forcing term A_s in Eq. (38) becomes zero, thus producing a drop of the dissipated power inside the planet clearly evident in the bottom panels of Figs. 4, 6, and 8.

Fig. 3 Upperpanel: rotation period of CoRoT-7 vs. the time as obtained from our model of the stellar angular momentum and tidal evolution introduced in Sect. 2.9. The vertical dotted line marks the time corresponding to the presently measured mean rotation period of the star that should correspond to its age in the hypothesis that its rotational evolution followed our model. Lower panel: orbital period of CoRoT-7b vs. the time as computed with our tidal model in Sect. 2.9. The vertical dotted line marks the current value of the orbital period.

Fig. 4 Upper panel: enlargement of the upper plot of Fig. 3 showing the rotation period of CoRoT-7 vs. the time as obtained from our model of the stellar angular momentum and tidal evolution introduced in Sect. 2.9. Middle panel: enlargement of the lower plot of Fig. 3, showing the orbital period of CoRoT-7b vs. the time. Lower panel: power dissipated inside CoRoT-7b vs. the time. The change in the slope of the curve showing the orbital period evolution occurs when Prot = 2Porb and corresponds to the limit beyond which inertial waves are no longer excited inside the star leading to a remarkable decrease of the tidal dissipation. The same condition corresponds to the resonant excitation of the radial oscillations and the intense dissipation inside the planet.

Fig. 5 Same as Fig. 3, but for Kepler-78.

Fig. 6 Same as Fig. 4, but for Kepler-78.

4 Discussion and conclusions

We found that a resonance between the stellar rotation and the orbital motion of an ultra-short-period planet (P_orb≲1.0 days) can occur during the first phase of the evolution of the star soon after (≈ 40−50 Myr) it reached the ZAMS. The resonance happens when the stellar rotation period is very close to twice the orbital period of the planet. This implies that the star must reach the ZAMS as a fast rotator with P_rot≲1−2 days for the resonance to occur. As a consequence, a large amount of heat (≈ 10³¹−10³² J) is dissipated in the planetary interior during a short time interval ((1−5) × 10⁴ yr), much shorter than the timescale for its transport by conduction or convection in the planetary interior (cf. Jaupart et al. 2007).

Immediately after their formation, telluric planets very likely experience a short phase of complete melting owing to the heat released in their accretion and core formation processes (Jaupart et al. 2007; Chao et al. 2021). The duration of this phase is limited to ≈ 10 Myr for an Earth-mass planet because of the rapid cooling of the interior magma ocean (e.g., Abe 1997; Jaupart et al. 2007). Therefore, we can assume that the mantles of our ultra-short-period planets have had enough time to become nearly solid when they experience the rapid and intense heating produced by the process introduced above, except for their outermost layers where stellar irradiation or magnetic induction can maintain a temperature higher than ~1300−1400 K sufficient to keep silicates liquid at low pressures (Chao et al. 2021).

The intense and rapid heating produces a huge increase of the temperature in the interior of the planet considering that the specific heat for an Earth-like mantle is ~ 1250 J kg K⁻¹. From the energy per unit mass listed in Table 3, we estimate a temperature increase of ~ 10³−10⁴ K, that is sufficient in several cases to completely melt the planetary interior even including a latent heat of mantle melting of ~ 3.2 × 10⁵ J kg⁻¹ (e.g., Abe 1997; Driscoll & Barnes 2015), thus resuming the state experienced by the planet immediately after its formation and core differentiation. The amount of energy E_D dissipated into the planet does not depend on its tidal quality factor $Q_{\text{p}}^{\prime \text{T}}$, therefore this conclusion appears to be robust.

The peak power reached by the dissipation mechanism we introduced is on the order of 10¹⁹ W and is huge in comparison with the present internal heat flow in our Earth that is ~ 47 ± 2 TW (Davies & Davies 2010). The internal heat flux coming from the decay of radioactive nuclids is presently estimated to be ~ 22 TW (O’Neill et al. 2020), although the recent detection of geoneutrinos would suggest a larger value of ~ 38 ± 14 TW (Agostini et al. 2020). An extrapolation to an age of about ~30−40 Myr, indicates a contribution of the radioactive heating between 100 and 200 TW, thus much smaller than the heat power released by our mechanism. The internal heat released during the process of the planet formation and the core differentiation is assumed to be dissipated for the most part during the first 10− 30 Myr of the planet evolution, therefore, we may assume that the mean temperature of the mantle of an average telluric planet at an age of 40− 50 Myr is only about 300−400 K higher that the present Earth value (Jaupart et al. 2007) implying a cooling power on the order of 30− 50 TW at those young ages. We note that the present power due to the cooling of the interior of our Earth is ≈ 25 TW, including ≈ 8 TW coming from the slow solidification of the inner metallic core of our planet (e.g., Table 11 of Jaupart et al. 2007).

The magnetic interaction between a telluric planet endowed with an intrinsic magnetic field and its young host star can lead to a release of power up to 10¹⁹ W through magnetic reconnection (e.g., Lanza 2012) or up to 10²¹ W in the case of the energy stored by Maxwell stresses in an interconnecting magnetic loop (e.g., Lanza 2013). However, such powers are released outside the planetary body at the boundary of its magnetosphere or close to its surface, respectively, because the intrinsic planetary field and the high electric conductivity of the ionized planetary atmosphere shield the planetary interior preventing the penetration of the stellar field inside the planet. Therefore, such kinds of dipolar magnetic interactions cannot heat the planetary interior and affect its evolution. The situation is different in the so-called unipolar inductor regime, when the planet has no significant intrinsic magnetic field and has a low, although finite, electric conductivity in its mantle, because in that case the stellar field can penetrate into its interior. In such a regime, the time variation in the stellar magnetic flux across the section of the planet can dissipate powers up to ≈ 10²²−10²³ W in the case of a young rapidly rotating star (e.g., Laine & Lin 2012), which can exceed the peak power dissipated by the mechanism proposed in the present paper.

The dissipation predicted by our mechanism mainly occurs in the planetary layers most affected by the tidal deformation, which implies that most of the heating is concentrated in the outer half of the planet. Therefore, it can profoundly affect the cooling of the planetary core, for instance, inhibiting its convective motions and shutting off the operation of the planet hydromagnetic dynamo until the previous temperature gradient is re-establish which may imply timescales on the order of at least 10⁷ −10⁸ yr. This may have a remarkable effect on the planetary atmosphere that can be completely eroded by the stellar high-energy radiation and wind, given that the protection by the planetary magnetic field has ceased. Moreover, the shutting off of the planetary dynamo may change the regime of the star-planet magnetic interaction to the unipolar inductor. Such a regime, characterized by a huge internal dissipation that persists for several tens of Myr until the stellar field weakens owing to the braking of stellar rotation, can deeply affect the evolution of the planetary interior and of its atmosphere (cf. Laine & Lin 2012).

From a more general point of view, a phase of very intense heating as predicted by our model (or the unipolar inductor model) can lead to changes in the chemistry of the interior with consequences for the properties of the surface and/or of the secondary planetary atmosphere that could be detectable, at least in principle. We do not explore these consequences here because they are outside the scope of this work that is devoted to introduce the physical basis of this new heating mechanism. The evolution of a global magma ocean and its interaction with a planetary atmosphere leading to alterations of a planet’s observable properties are discussed in, e.g., Chao et al. (2021).

Given that not all late-type dwarf stars begin their ZAMS evolution as fast rotators and that ultra-short-period planets may not be already in place during the early phases of their main-sequence evolution, not all these planets may experience this phase of early and highly intense heating. A future analysis of the surface properties and of the secondary atmospheres of ultra-short-period planets may reveal the signatures of that phase, thus confirming the prediction about such a new internal heating mechanism. In any case, the main result of the present work suggests that some ultra-short-period planets may have had an evolution remarkably different from that of their siblings on longer period orbits or that migrated close to their host stars after the first ≈ 100 Myr of their main-sequence lifetime.

Fig. 7 Same as Fig. 3, but for K2-141.

Fig. 8 Same as Fig. 4, but for K2-141.

Acknowledgements

The author gratefully acknowledges an anonymous referee for a careful reading of the manuscript and several valuable comments that helped to improve this work. This investigation has been supported by the PRIN-INAF 2019 “Planetary systems at young ages (PLATEA)”, the PI of which is Dr. S. Desidera. The use of the MESA-Web interface in December 2020 to compute illustrative stellar models is also acknowledged.

Appendix A Computation of some integrals with the method of the residues

To compute the integrals appearing in Eqs. (48) and (51), we made use of the theorem of the residues (Smirnov 1964). We shall consider a complex rational function f(z) of a complex argument $z\in ℂ$ defined as $$f(z)=\frac{\varphi (z)}{\psi (z)},$$(A.1)

where ϕ(z) and ψ(z) are complex polynomials and the degree of ψ(z) is greater thanthe degree of ϕ(z) by at least two, while ψ(z) has no zeroson the real axis.

We consider the integral of f(z) along the closed curve Σ consisting of the union of the semicircle C in the positive imaginary semiplane having its center at the origin O and radius R with its diameter, that is, the segment of the real axis [−R, R] (see Fig. A.1). The theorem of the residues states that (Smirnov 1964, Sect. 58) $${\displaystyle {\int }_{\Sigma }f}(z)\text{d}z=2\pi i{\sum }_k{r}_k,$$(A.2)

where $i=\sqrt{-1}$ and the r_k are the residues of f(z) in the poles enclosed by the closed curve Σ. We can write the integral along Σ as the sum of the integral along the segment on the real axis plus the integral along C as $${\displaystyle {\int }_{\Sigma }f}(z)\text{d}z={\displaystyle {\int }_{-R}^{+R}f}(x)\text{d}x+{\displaystyle {\int }_{C}f}(z)\text{d}z,$$(A.3)

where f(x) is the restriction of the complex function f(z) to the real axis. If we make R →∞, the integral along C tends to zero becausethe function f(z) tends to zero at least as z⁻² for z →∞ and we can apply Eq. (A.2) to find: $${\displaystyle {\int }_{-\infty }^{\infty }f}(x)\text{d}x=2\pi i{\sum }_{\Im z_k>0}{r}_{z_k},$$(A.4)

where the summation is extended to all the residues of f(z) in the upper complex semiplane, that is, to all and only the residues $r_{z_k}$ in the poles z_k of f(z) whose imaginary parts are positive. The improper integral in Eq. (A.4) exists because f(x) goes to zero as x⁻² when x → ± ∞ and its denominator has no zeros on the real axis.

Since f(z) is the ratio of two polynomials, its poles are the zeros of the denominator ψ(z). In all the cases considered in this Appendix, these are simple poles (that is, roots of the first order of the polynomial ψ(z)), therefore, the residual in a pole z_k is given by $$r_k={\left[f(z)(z-z_k)\right]}_{z=z_k}=\frac{\varphi (z_k)}{{\psi }^{\prime }(z_k)},$$(A.5)

where ψ′(z) is the derivative of the polynomial ψ(z) (Smirnov 1964, Sect. 21).

A.1 Computation of the integral in Eq. (48)

In this case, ϕ(z) = z² and $\psi (z)={(n^2-z^2)}^2+4b^2{z}^2$, where z is a complex argument, while n and b are real numbers. By differencing the latter, we find ψ′(z) = 4(z² + 2b² − n²)z. The poles are given by ψ(z) = 0, that is: $${\nu }_{1,2}=ib\pm \sqrt{n^2-b^2},$$(A.6) $${\nu }_{3,4}={\nu }_{1,2}^{*}=-ib\pm \sqrt{n^2-b^2}.$$(A.7)

From these equations, it immediately follows that ${\nu }_2=-{\nu }_1^{*}$ and ${\nu }_4=-{\nu }_3^{*}$, where the asterisk indicates complex conjugation.

Fig. A.1 Complex plane with the semicircle C (in red) and its diameter, that is, the segment of the real axis [−R, +R] (in green). Their union forms the closed contour Σ along which the complex function f(z) is integrated in Eq. (A.2). For the sake of illustration, we show two poles of the complex integrand function P1 and P2 (in blue encircled in red), inside the contour Σ and two poles outside the contour Σ. Only the residues in P1 and P2 are to be considered in the computation of the integral of f(z) along Σ according to the theorem of the residues, while the residues in the poles P3 and P4 (in blue), lying outside Σ, are not to be included in the summation in Eq. (A.2).

The two poles in the positive imaginary semiplane are ν_1,2. We have $\varphi ({\nu }_2)={[\varphi ({\nu }_1)]}^{*}$, while ${\psi }^{\prime }({\nu }_2)=-{[{\psi }^{\prime }({\nu }_1)]}^{*}$. Considering Eq. (A.5), this implies that the residues in ν₁ and ν₂ are related by $r_{{\nu }_2}=-r_{{\nu }_1}^{*}$. Using this result and after some algebra, we find that the sum of the residues in Eq. (A.4) is: $${\sum }_{\Im {\nu }_k>0}{r}_{{\nu }_k}=r_{{\nu }_1}+r_{{\nu }_2}=2i\Im r_{{\nu }_1}=-\frac{i}{4b},$$(A.8)

where ℑz denotes the imaginary part of the complex number z. By applying Eq. (A.4), finally we find $${\displaystyle {\int }_{-\infty }^{\infty }\frac{{\nu }^2\text{d}\nu }{{(n^2-{\nu }^2)}^2+4b^2{\nu }^2}}=\frac{\pi }{2b}.$$(A.9)

A similar integral, that is useful to compute ⟨x²⟩ for our applications (cf. Eq. (42)), can be evaluated with the same method: $${\displaystyle {\int }_{-\infty }^{\infty }\frac{\text{d}\nu }{{(n^2-{\nu }^2)}^2+4b^2{\nu }^2}}=\frac{\pi }{2b{n}^2}.$$(A.10)

A.2 Computation of the integral in Eq. (51)

The integral in Eq. (51) can be cast in the form $${\displaystyle {\int }_{-\infty }^{\infty }\frac{a+{\gamma }^2}{{\left(b+{\gamma }^2\right)}^2+c}}\text{d}\gamma ,$$(A.11)

where a = n² + D²∕4 > 0, b = −n² + D²∕4, and c = a² − b² = D²n² > 0, are real numbers.

In this case, the functions ϕ(z) = a + z² and $\psi (z)={(b+z^2)}^2+c$. The poles of the complex integrand are the four roots of ψ(z) = 0. To find their expressions, we consider the complex number $w=-b+i\sqrt{c}$. Its trigonometric expression is w = ρexp(iθ), where $\rho =\sqrt{b^2+c}=a>0$, cos θ = −b∕a, and $\sin \theta =\sqrt{c}/a$.

The four roots of ψ(z) = 0 are ${\gamma }_{1,2}=\sqrt{w}$ and ${\gamma }_{3,4}=\sqrt{w^{*}}$. The square root of w can be obtained from its trigonometric expression as $\sqrt{w}=\sqrt{a}\exp [i(\theta +2k\pi )/2]$, where k is an integer. Taking into accountthat $\sin (\theta /2)=\sqrt{(1-\cos \theta )/2}$ and $\cos (\theta /2)=\sqrt{(1+\cos \theta )/2}$, we find the four poles off(z) $$\begin{array}{l}{\gamma }_{1,2}=\pm \frac{1}{\sqrt{2}}\left(\sqrt{a-b}+i\sqrt{a+b}\right),\\ {\gamma }_{3,4}=\pm \frac{1}{\sqrt{2}}\left(\sqrt{a-b}-i\sqrt{a+b}\right).\end{array}$$(A.12)

The two poles in the positive imaginary semiplane are γ₁ and γ₄. They are related by ${\gamma }_1=-{\gamma }_4^{*}$ and, by applying Eq. (A.5), where ψ′(z) = 4z(b + z²), this implies $r_{{\gamma }_1}=-r_{{\gamma }_4}^{*}$. Therefore, the sum of the residues in the positive imaginary semiplane is: $${\sum }_{\Im r_{{\gamma }_k}>0}{r}_{{\gamma }_k}=r_{{\gamma }_1}+r_{{\gamma }_4}=2i\Im r_{{\gamma }_1}=-\frac{\sqrt{2}}{2}i\frac{1}{\sqrt{a+b}}.$$(A.13)

Finally, by applying the theorem of the residues (Eq. (A.4)), we find $${\displaystyle {\int }_{-\infty }^{\infty }\frac{a+{\gamma }^2}{{\left(b+{\gamma }^2\right)}^2+c}}\text{d}\gamma =\frac{\pi \sqrt{2}}{\sqrt{a+b}}.$$(A.14)