© The Authors 2023

1 Introduction

Red giant stars with an initial mass similar to that of the Sun reach a luminosity of thousands of solar luminosities and a radius of a few hundred solar radii close to the tip of the red giant branch (RGB; e.g., Kippenhahn et al. 2012). Convective motions inside the extended envelopes of such stars transport virtually all those large luminosities from the energy-generation layers up to the photosphere. Density fluctuations are associated with convection because rising hotter columns of plasma are less dense than the average, while cool descending columns are denser. Such density fluctuations with a relative amplitude of up to ≈10⁻⁵ cause the red giant to have a slightly fluctuating outer gravitational potential that can affect the orbit of a close-by companion. This effect was proposed by Phinney (1992) and Phinney & Kulkarni (1994) to account for the nonzero eccentricity of binary millisecond pulsars where very small eccentricities (e ≲ 10⁻⁵) can be measured thanks to the exquisite precision of measurements of the orbital Doppler shift, which is made possible by the periodic pulses coming from the rapidly spinning neutron star.

According to these latter authors, such low-eccentricity millisecond pulsars are the end products of low-mass X-ray binaries, where the donor star of the binary evolves off the main sequence, becoming a red giant, and thus starts to transfer mass and angular momentum to a very old neutron star companion. This produces the X-ray emission and spins up the neutron star to a rotation period in the millisecond range. Once the mass transfer ceases, we observe the system as a millisecond binary pulsar at radio wavelengths (see, e.g., Tauris & Savonije 1999; Tauris & van den Heuvel 2006; Chen et al. 2021, for details and discussion of the different formation channels for those systems).

The evolution phase of such systems that is relevant to our line of enquiry occurs when the red giant approaches the end of its nuclear fuel reservoir and starts to contract, detaching from its Roche lobe, and thus ending the mass transfer to the neutron star. In this phase, the orbit of the pulsar evolves under the action of tides in the contracting red giant and the density fluctuations in its convective envelope, which produce a slightly fluctuating outer gravitational potential. This induces a small eccentricity in the orbit of the binary that tides are not capable of completely erasing because the radius of the secondary is drastically reduced on a timescale shorter than the damping timescale of the residual eccentricity. This residual eccentricity therefore remains as a relic of the final phase of the evolution of the detached secondary before it becomes a WD; it is a stochastic variable because it is produced by the random density fluctuations inside the convective envelope of the contracting red giant (cf. Phinney 1992).

The computation of the residual eccentricity in the case of millisecond pulsar binaries was revisited by Lanza & Rodonò (2001), who applied the hydromagnetic virial theorem of Chandrasekhar (1961) and methods of the stochastic differential equation theory to compute its statistical distribution. Moreover, Lanza & Rodonò (2001) suggested the relevance of the residual eccentricity for planets orbiting red giant stars, but did not investigate the topic (see a parallel treatment in Sect. 5 of Kissin & Thompson 2015).

Because the Sun will become a red giant, this kind of investigation can be used to predict the late phases of the evolution of our planetary system. The evolution of the orbits of the inner Solar System planets – including the Earth – during the remaining lifetime of the Sun on the main sequence cannot be predicted in a deterministic way because of the chaotic character of their long-term dynamics (e.g., Laskar 1996; Mogavero & Laskar 2021). Numerical integrations of thousands of Solar System realizations starting from its present dynamical state within the range of its uncertainty and simplified analytical models show that the planet Mercury has a probability on the order of 0.5–1% of colliding with Venus or with the Sun within the next 5 Gyr (Laskar & Gastineau 2009; Batygin et al. 2015). While most of the highly eccentric orbits of these models of Mercury do not destabilize the inner Solar System, a small fraction of them could result in the excitation of large changes in the orbit eccentricities of the Earth, Mars, or Venus, leading to close encounters or collisions between these planets. The probability of such a catastrophic outcome for our planet before the Sun leaves the main sequence is difficult to estimate, but it is likely to be below 0.1% (e.g., Laskar & Gastineau 2009; Zeebe 2015). Therefore, in our subsequent considerations, we assume that the eccentricity of the orbit of the Earth is not significantly modified by the perturbations of the other planets during the entire lifetime of the Solar System. In other words, the evolution of the Sun appears us to be the most relevant cause of changes in the orbit of the Earth and will most likely determine its final fate.

Sackmann et al. (1993) explored the evolution of our star in detail and considered the orbital expansion produced by the solar mass loss, but did not include the effects of tides in their model. We note that the orbital expansion when the star loses its mass is a consequence of the conservation of the angular momentum. On the other hand, tides counteract this expansion and damp the eccentricity of the orbit (cf. Sect. 2.1).

Rybicki & Denis (2001) re-examined the fate of the Solar System by discussing several other models for solar evolution and including tidal effects. These authors showed that the drag exerted by the solar wind, the mass accretion from the wind, and the evaporation of our planet due to solar irradiation have negligible effects on the evolution of the Earth’s orbit, even during the RGB phase and up to the tip of the asymptotic giant branch (AGB) phase. Therefore, the most relevant effect in determining whether the Earth will be finally engulfed by the Sun or not appears to be the solar mass-loss rate because it rules both the largest radii reached on the RGB and AGB phases and the orbital expansion (see also Guo et al. 2016).

Schröder & Smith (2008) adopted a specific model to predict solar mass loss and found that the Sun will reach its maximum radius at the end of the RGB phase. Including the effects of tides, these authors found that the Earth will not be able to escape engulfment at the tip of the RGB phase despite the expansion of its orbit.

The consideration of the residual eccentricity of the Earth’s orbit produced by the density fluctuations in the extended envelope of the red giant Sun will add an important conceptual difference to the final phases of the orbital evolution of our planet. While in principle the time of the engulfment of the Earth can be predicted by models such as those of Rybicki & Denis (2001) or Schröder & Smith (2008), provided both an accurate stellar evolution model and tidal theory, this is no longer possible with the introduction of the residual eccentricity because it is a stochastic variable for which we can only compute a statistical probability distribution. The difference from a practical point of view is negligible because the final fate of the Earth will not change, but the residual eccentricity may play a role in the final phases of the dynamical evolution of the Solar System because it can induce secular changes in the eccentricities of the orbits of small bodies such as asteroids or comets. This is not possible when the orbit of the Earth is circularized during the late phase of its evolution as predicted on the basis of tidal theory alone. Putting an asteroid on a highly eccentric orbit may produce its capture by the WD Sun, leading to its fall on the degenerate star, an event that will pollute its atmosphere with metals, as observed in at least ~25% of WDs (cf. Frewen & Hansen 2014; Veras 2021, and references therein, especially the list in their Fig. 6).

In the present work, we determine the statistical distribution of the residual eccentricity of the orbit of an Earth-like planet along the evolution of a Sun-like star up to its final phases. Moreover, if the planet is transiting across the disc of a red giant, the fluctuation of its outer gravitational potential affects the timing of the transits in a potentially measurable way. We revisited the methodology used to compute the residual eccentricity and the transit-time variations and developed a simpler approach than that of Lanza & Rodonò (2001) that allows a straightforward derivation of all the relevant equations and includes the effects of stellar mass loss and tides during the evolution of the star in its red giant phase. The model we use to compute the statistical distribution of the residual eccentricity and the transit-time variations is presented in Sect. 2, while an internal structure and evolution model of a Sun-like star is introduced in Sect. 3. We present an illustrative application of our approach to an Earth-Sun-like system in Sect. 4. We discuss our results and conclude in Sect. 5.

2 Model

In this section, we introduce the model to compute the residual eccentricity and the transit-time fluctuations, while its application to an Earth-like planet in orbit around a Sun-like star is presented in Sect. 4. The effects of the stellar mass loss and of the tides on the planet orbit are discussed in Sect. 2.1, while the computation of the outer gravitational potential of the star is presented in Sect. 2.2. The variations in the orbital elements of the planet under the perturbations produced by such a fluctuating potential are introduced in Sect. 2.3, that is, separately for the eccentricity (Sect. 2.3.1) and the mean longitude at the epoch from which the transit-time variations can be derived (Sect. 2.3.2). The power spectrum of the gravitational potential fluctuations, which is required to compute the fluctuations in the orbital elements, is computed in Sect. 2.4.

The evolution of the orbit of a planet around a RGB star has been investigated in order to account for the observed lack of giant planets with a semimajor axis of shorter than ~0.55 au around giant stars (cf. Sect. 4 of Villaver et al. 2014). The minimum initial orbital semimajor axis that ensures the survival of a planet of given mass has been determined under the effects of tides and stellar mass loss considering stars of different masses and adopting different prescriptions for their mass-loss rates and several evolution codes. Initially, most of the investigations were limited to solar-like or intermediate-mass stars evolving from the main sequence up to the RGB (e.g., Villaver & Livio 2009; Kunitomo et al. 2011; Villaver et al. 2014). Subsequent works investigated the fate of planets around AGB stars, including the thermal pulse phase (Mustill & Villaver 2012) and also considering companions up to the stellar mass range (e.g., Nordhaus & Spiegel 2013) in an attempt to predict the formation of planetary or stellar binary systems with a WD central star. Similar models have been computed to predict the fraction of planetary nebulae with central binary systems (cf. Madappatt et al. 2016). Most of those studies focused on systems consisting of only a single planet (however, see the discussion in Sect. 5.2 of Mustill & Villaver 2012), while an investigation of the role of the mutual gravitational perturbations in a system consisting of an inner Neptune-mass planet and an outer Jupiter-mass planet around a red giant star was carried out by (Ronco et al. 2020).

The most important effects to influence the orbits of late planets have been found to be the stellar mass loss and the tides raised by the planets inside their giant host stars (see Adams & Bloch 2013, for an analytic model of their joint effects). On the other hand, the impact of giant planets on the evolution of their host stars – in particular on their angular momentum, mass loss, and hydromagnetic dynamo action – has also been investigated (e.g., Soker & Harpaz 2000; Carlberg et al. 2009; Privitera et al. 2016a,b,c).

In the present investigation, we introduce a model to compute the residual eccentricity and the transit-time fluctuations and show a simple application to the case of a system consisting of a single Earth-mass planet orbiting a Sun-like star because a general investigation of the fates of planets and planetary systems around RGB and AGB stars is well beyond the scope of the present work. Therefore, we refer the interested reader to the above papers for information on such a wider scenario.

2.1 Stellar mass loss and tides

Let us consider a system consisting of a planet of mass m_p orbiting a Sun-like star of mass m_s; the semimajor axis of the orbit is a, while its eccentricity is e. The orbital angular momentum is given by $$J=m\sqrt{G{m}_{\text{T}}a\left(1-e^2\right)}=mn{a}^2\sqrt{1-e^2},$$(1)

where m ≡ m_pm_s/m_T is the reduced mass of the system, m_T ≡ m_s + m_p its total mass, and we make use of Kepler’s third law: $$n^2{a}^3=G{m}_{\text{T}},$$(2)

where n = 2π/P_orb is the orbital mean motion, with P_orb being the orbital period. Given that m_p ≪ m_s and that the orbit is nearly circular (e ≪ 1), the conservation of the orbital angular momentum allows us to find the variation in the orbit semimajor axis produced by a variation in the mass of the star as $$\frac{1}{a}\frac{\text{d}a}{\text{d}t}\simeq -\frac{1}{m_{\text{s}}}\frac{\text{d}{m}_{\text{s}}}{\text{d}t},$$(3)

which can be used to compute the dynamical effect of the stellar mass loss (cf. Schröder & Smith 2008, Sect. 2).

Equation (3) is obtained by considering an isotropic stellar mass loss (Veras et al. 2011). Rigorously speaking, any mass loss perturbs the orbit of the star-planet system by producing variations in its orbital elements that can no longer be regarded as constant as in the case of the standard two-body problem. Osculating orbital elements can be introduced to describe the instantaneous orbit under the effects of the mass-loss perturbation (e.g., Dosopoulou & Kalogera 2016a). From a dynamical point of view, the mass-loss regime can be characterized through the parameter (cf. Veras et al. 2011) $${\text{Ψ}}_{\text{ml}}\equiv \frac{1}{n{m}_{\text{T}}}\frac{\text{d}{m}_{\text{T}}}{\text{d}t}\simeq \frac{1}{n{m}_{\text{s}}}\frac{\text{d}{m}_{\text{s}}}{\text{d}t}.$$(4)

When Ψ_ms ≪ 0.1, the mass-loss timescale is much longer than the orbital period and the system is in the so-called adiabatic mass-loss regime, which allows us to average the variations in the osculating elements along the orbit and obtain variation equations that are independent of the true anomaly. It is in this regime that Eq. (3), rigorously speaking, is valid. It is interesting to note that in the case of an anisotropic stellar mass loss in the adiabatic regime, the semimajor axis remains secularly constant because the effects of the anisotropic mass flux average to zero (Veras et al. 2013; Dosopoulou & Kalogera 2016b). Therefore, only the isotropic mass loss is to be included into Eq. (3).

The variation in the eccentricity due to an isotropic stellar mass loss, averaged along the orbit, becomes zero when computed in the adiabatic regime. Nevertheless, the eccentricity as an osculating element oscillates along the orbit with an amplitude of Ψ_ml(1 – e²), where e is the eccentricity in the absence of mass loss (cf. Veras et al. 2011, Sect. 2.4.1). Even assuming a mass loss of 0.1 M_⊙ in 1 Myr by a solar-mass star and a planet with an orbital period of 2 yr, the parameter Ψ_ml ~ 1.3 × 10⁻⁶ and the oscillation of the eccentricity along the orbit is negligible in comparison to the residual eccentricity induced by convective fluctuations in the stellar envelope, which reaches values of up to ~10⁻² in the late stages of stellar evolution (cf. Sect. 4). Therefore, we can safely assume for our purposes that the eccentricity is not significantly excited by isotropic mass loss from the star.

The second contribution to be considered in modeling the evolution of the orbit is the tidal torque that affects both the semimajor axis a and the eccentricity e. The tidal torque scales as (R_s/a)⁶, where R_s is the radius of the star, and therefore tides will be relevant only after the star has started to ascend the RGB, because it is only in that phase that its radius will become sufficiently large to produce a significant tidal torque in the case of an Earth-mass planet on an Earth-like orbit (a ~ 1 au). Following Zahn (1989), we write the tidal torque Γ acting on the planetary orbit by the equilibrium tide as (cf. Schröder & Smith 2008) $$\text{Γ=6}\frac{{\lambda }_2}{t_{\text{f}}}{\left(\frac{m_{\text{p}}}{m_{\text{s}}}\right)}^2{m}_{\text{s}}{R}_{\text{s}}^2{\left(\frac{R_{\text{s}}}{a}\right)}^6\left({\text{Ω}}_{\text{s}}-n\right),$$(5)

where λ₂ ≃ 0.019 α^4/3 and $t_{\text{f}}={\left(m_{\text{s}}{R}_{\text{s}}^2/L_{\text{s}}\right)}^{1/3}$ are a nondimen-sional parameter and the convective friction timescale appearing in the equilibrium tide theory of Zahn (1977, 1989), respectively, while Ω_s is the angular velocity of rotation of the star, L_s its luminosity, and α the ratio of the mixing length to the pressure scale height.

A Sun-like star on the RGB can be assumed to be almost nonrotating because of the loss of angular momentum during the evolution on the main sequence and the subgiant phases and the large increase in its moment of inertia due to the radius expansion. In other words, we can assume Ω_s ~ 0 for our purposes, and therefore Γ < 0, which produces a decrease in the orbit semimajor axis. Considering a planet with an orbital period of longer than ~ 150–200 days, the frequency of the equilibrium tide is smaller than the turnover frequency of convective eddies inside the star, which justifies the use of an unreduced turbulent viscosity in the convective envelope when estimating the tidal dissipation (Zahn 1989). Moreover, our simplified expression of λ₂ in comparison with that in Eq. (15) of Zahn (1989) is justified in the case of an orbital period of around 1 yr because t_f ≈ 1 yr for the Sun close to the tip of the RGB. The effects of dynamical tides can be neglected for a solar-mass star when the orbital separation is larger than 0.2–0.3 au, as in our case (see Rao et al. 2018; Sun et al. 2018, for a discussion of the role of dynamical tides).

Adding the effect of the tidal torque to that of the stellar mass loss, Eq. (3) becomes $$\frac{1}{a}\frac{\text{d}a}{\text{d}t}\simeq -\frac{1}{m_{\text{s}}}\frac{\text{d}{m}_{\text{s}}}{\text{d}t}+\frac{2}{J}\text{Γ,}$$(6)

which we use to compute the evolution of the semimajor axis of the planetary orbit in our application in Sect. 4.

Another effect of the equilibrium tide inside the red giant star is the damping of the orbital eccentricity of the planet. We compute the circularization time τ_e according to Verbunt & Phinney (1995), who calibrated the tidal theory of Zahn by means of a sample of binaries containing a giant component. Specifically, we shall adopt $$\frac{1}{{\tau }_{\text{e}}}\equiv -\frac{1}{e}\left(\frac{\text{d}e}{\text{d}t}\right)=\beta {\left(\frac{L_{\text{s}}}{m_{\text{env}}{R}_{\text{s}}^2}\right)}^{1/3}\left(\frac{m_{\text{env}}}{m_{\text{s}}}\right)\left(\frac{m_{\text{p}}}{m_{\text{s}}}\right)\left(\frac{m_{\text{T}}}{m_{\text{s}}}\right){\left(\frac{R_{\text{s}}}{a}\right)}^8,$$(7)

where ß is a nondimensional empirical factor of order unity and m_env the mass of the convective envelope of the giant star where the kinetic energy of the tidal flow is dissipated. By fitting their binary sample, Verbunt & Phinney (1995) found 0.5 ≲ ß ≲ 2, that is, Eq. (7) provides an estimate of τ_e within approximately a factor of two when we adopt ß = 1.

Equation (7) does not include the variation in the eccentricity produced by an anisotropic stellar mass loss, an effect that cannot be excluded in the case of solar-like stars in the late stages of their evolution (e.g., Soker 1998). The impact of this mass loss on the eccentricity depends on its time dependence and on the specific geometry. Several idealized cases have been explored by Veras et al. (2013) and Dosopoulou & Kalogera (2016b), for example. The most relevant occurs when there is a longitudinal anisotropy in the mass loss because of a difference in the mass-loss rate between the stellar hemispheres of up to 1%, with an average mass-loss rate of ~3 × 10⁻⁷ M_⊙ yr⁻¹. In this case, the anisotropy can excite an eccentricity on the order of 0.01 (e.g., Veras et al. 2013, Sect. 3.1.3). Such eccentricity adds to the residual eccentricity that we are considering in the present work and should be included as an additional effect when treating stars with an anisotropic stellar mass loss.

In addition to the effects of stellar mass loss and tides, the orbit of the planet can be affected by the frictional and gravitational drag forces. Frictional drag occurs because the planet is an obstacle that moves inside the stellar wind, while the gravitational drag arises because of the gravitational field of the planet acting on the wind itself. Their effect is small in comparison to those of the tides and stellar mass loss (cf. Duncan & Lissauer 1998; Villaver & Livio 2009; Mustill & Villaver 2012; Villaver et al. 2014; Rao et al. 2018; Yarza et al. 2022), and therefore we neglect them in our simplified model and assume that the mass of the planet is constant throughout its evolution. This amounts to neglecting planetary evaporation and mass accretion by the stellar irradiation and the stellar wind, respectively.

2.2 Gravitational potential of a nonspherically symmetric star

The gravitational potential Φ of a star at an external point P up to the quadrupole order can be expressed as $$\text{Φ\hspace{0.17em}=\hspace{0.17em}}-\frac{G{m}_{\text{s}}}{r}-\frac{3G}{2r^3}{\displaystyle \sum _{i,k}\frac{Q_{ik}{x}_i{x}_k}{r^2}},$$(8)

where G is the gravitation constant, m_s the mass of the star, r > R_s the distance of P from the barycenter O of the star, Q_ik the quadrupole moment tensor of the star, and x_i the coordinates of P in an orthogonal 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\hspace{0.17em}}I,$$(9)

where δ_ik = 1 for i = k and δ_ik = 0 for i ≠ k is the Kronecker δ tensor, and Tr I = 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 \left(r\right)x_i}{x}_k\text{\hspace{0.17em}d}V,$$(10)

where ρ is the density, r ≡ (x₁, x₂, x₃) the position vector, and V the volume of the star over which the integration is extended.

The density inside a stellar convective envelope can be written as $$\rho \left(r,t\right)={\rho }_0\left(r\right)+\rho \prime \left(r,t\right),$$(11)

where ρ₀(r) is the mean density at position r and ρ′ (r, t) is the fluctuation that depends on both the position and the time t. The components of the moment of inertia become $$I_{ik}\left(t\right)={\displaystyle {\int }_V\rho \left(r,t\right)x_i}{x}_k\text{\hspace{0.17em}d}V={\displaystyle {\int }_V\left[{\rho }_0\left(r\right)+\rho \prime \left(r,t\right)\right]}{x}_i{x}_k\text{d}V\equiv I_{ik}^{\left(0\right)}+I_{ik}^{\prime },$$(12)

where we have defined the time average of the component and its fluctuation. Because the average density is spherically symmetric, only the diagonal components of the $I_{ik}^{\left(0\right)}$ inertia tensor are different from zero. In other words, $$\begin{array}{l}{I}_{xx}^{\left(0\right)}=I_{yy}^{\left(0\right)}=I_{zz}^{\left(0\right)}=\frac{1}{3}{\displaystyle {\int }_V{\rho }_0\left(r\right)r^2\text{\hspace{0.17em}d}V,}\hfill \\ I_{ik}^{\left(0\right)}=0\text{\hspace{0.17em}for\hspace{0.17em}}i\ne k,\hfill \end{array}$$(13)

and $$\text{Tr\hspace{0.17em}}{I}^{\left(0\right)}=3I_{xx}^{\left(0\right)}=3I_{yy}^{\left(0\right)}=3I_{zz}^{\left(0\right)}.$$(14)

Therefore, only the fluctuating parts ${I^{\prime }}_{ik}$ contribute to the quadrupole moment tensor Q_ik. It is a symmetric tensor, and therefore it is possible to reduce it to a diagonal form by a suitable orientation of the coordinate axes (e.g., Goldstein 1950). Such axes are the principal axes of inertia of the star and their orientation varies randomly because of the density fluctuations.

Because the trace of the quadrupole moment tensor is zero according to Eq. (9), only two of its diagonal components are independent, and we can write its nonvanishing components in the reference frame of the principal axes by introducing two scalars, Q and T, as $$\{\begin{array}{l}{Q}_{xx}\equiv Q+T/2\hfill \\ Q_{yy}\equiv Q+T/2\hfill \\ Q_{zz}\equiv -2Q.\hfill \end{array}$$(15)

Therefore, the gravitational potential becomes $$\begin{array}{l}\text{Φ\hspace{0.17em}=\hspace{0.17em}}-\frac{G{m}_{\text{s}}}{r}-\frac{3G}{2r^3}\frac{Q_{xx}{x}^2+Q_{yy}{y}^2+Q_{zz}{z}^2}{r^2}\hfill \\ \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}}-\frac{G{m}_{\text{s}}}{r}-\frac{3G}{2r^3}Q-\frac{3G}{2r^5}\left[\frac{1}{2}T\left(x^2-y^2\right)-3Q{z}^2\right],\hfill \end{array}$$(16)

where we have separated the isotropic component of the quadrupole potential from the component depending on the coordinates x = x₁, y = x₂, z = x₃ of the point P in the reference frame of the principal axes. We note that such coordinates continuously change in time because of the random reorientation of the principal axes. Because the star is nonrotating, the statistical distribution of the orientation of the principal axes will be isotropic thanks to the spherical symmetry of the density fluctuations, the distribution of which will depend only on the distance r from the center of the star.

In Sect. 2.4, we compute the autocorrelation function of the quadrupole potential and use it to compute the residual eccentricity of the orbit and the transit-time variations. Such an autocorrelation cannot depend on the orientation of the principal axes because of the isotropy of the density fluctuations. Therefore, it cannot depend on the terms in square brackets on the right-hand side of Eq. (16) because they depend on the coordinates of the point P where the potential is evaluated in the reference frame of the principal axes¹. In other words, we can limit ourselves to considering only the isotropic part of the quadrupole potential in our model, which greatly simplifies our treatment: $$\text{Φ\hspace{0.17em}=\hspace{0.17em}}-\frac{G{m}_{\text{s}}}{r}-\frac{3G}{2r^3}Q.$$(17)

2.3 Perturbations of the planetary orbit

The fluctuating quadrupole gravitational potential induces perturbations of the orbital elements of the planet. As such perturbations are small, we assume an unperturbed circular orbit for reference and use the Gauss equations to compute the time derivatives of the orbital elements (e.g., Roy 1978, Sect. 6.7.4). From Eq. (17), it follows that the perturbative acceleration acting on the planet is purely radial. We indicate such a radial acceleration with S, while the components of the perturbative acceleration in the orbital plane and normal to the radius vector, T, and normal to the orbital plane, W, are both zero: $$S=-\frac{\partial \text{Φ}}{\partial r}=-\frac{9G}{2r^4}Q\left(t\right),$$(18) $$T=0,$$(19) $$W=0.$$(20)

The only two relevant Gauss equations for the variation in the orbital elements are $$\frac{\text{d}e}{\text{d}t}=\frac{S}{na}\sin nt-\frac{e}{{\tau }_{\text{e}}},$$(21) $$\frac{\text{d}ϵ}{\text{d}t}=-\frac{2S}{na},$$(22)

where e is the eccentricity and ϵ the mean longitude at the epoch of the planetary orbit, a its semimajor axis, n the mean motion, τ_e the tidal decay timescale of the eccentricity (cf. Eq. (7)), and t the time. Equations (21) and (22) are written for an unperturbed circular orbit with radius equal to a, while the action of the tides inside the star – which tends to damp the eccentricity – is represented by the term −e/τ_e. Such a term is not included in the Gauss equation for the eccentricity, but is added here to take into account the effect of the stellar tides.

The mean longitude of the planet 1_p is related to ϵ by (cf. Roy 1978) $$l_{\text{p}}=ϵ+nt,$$(23)

in the case of a circular reference orbit. As n can be assumed fixed on the timescales over which ϵ fluctuates (see below), the fluctuations in the mean longitude are ∆l_p = ∆ϵ, which can be used to compute the variations in the time of midtransit in the case of a planet that transits across the disk of our red giant star as $$O-C=\frac{\text{Δ}ϵ}{2\pi }{P}_{\text{orb}},$$(24)

where O – C is the difference in the time of midtransit with respect to an unperturbed orbit of constant orbital period P_orb = 2π/n.

2.3.1 Residual eccentricity

We recast Eq. (21) in the form $$\frac{\text{d}e}{\text{d}t}+\frac{e}{{\tau }_{\text{e}}\left(t\right)}=f\left(t\right),$$(25)

where we make the dependence of the tidal damping timescale τ_e on the time explicit (cf. Eq. (7)) and define $$f\left(t\right)\equiv -KQ\left(t\right)\sin nt,$$(26)

where $$K\equiv \frac{9G}{2n{a}^5}=\frac{9n}{2m_{\text{T}}{a}^2},$$(27)

with the second expression for K coming from the application of Kepler’s third law. By applying the method of the variation of the constants, the general solution of Eq. (25) is $$e\left(t\right)C_0\zeta \left(t\right)+\zeta \left(t\right){\displaystyle {\int }_0^{t}f\left(t\prime \right){\left[\zeta \left(t\prime \right)\right]}^{-1}\text{\hspace{0.17em}d}t\prime ,}$$(28)

where C₀ is a constant depending on the initial value of e and $$\zeta \left(t\right)\equiv \exp \left[-{\displaystyle {\int }_0^t\frac{\text{d}t\prime }{{\tau }_{\text{e}}\left(t\prime \right)}}\right].$$(29)

The ensemble mean value of the eccentricity 〈e〉 = 0 because f(t) is a stochastic function, while the mean value of its square can be obtained as $$\langle e^2\left(t\right)\rangle =\langle {\left[\zeta \left(t\right)\right]}^2{\displaystyle {\int }_0^t{\displaystyle {\int }_0^{t}f\left(t\prime \right)f\left(t″\right)}}{\left[\zeta \left(t\prime \right)\right]}^{-1}{\left[\zeta \left(t″\right)\right]}^{-1}\text{\hspace{0.17em}d}t\prime \text{d}t″\rangle .$$(30)

The term containing the constant factor C₀ in Eq. (28) becomes negligible after a sufficiently long interval of time because it decreases exponentially according to Eq. (29). For this reason, the terms containing C₀ have been dropped from Eq. (30). The ensemble mean commutes with the integration, allowing us to evaluate 〈e²(t)〉 if we know 〈f(t′)f(t″)〉.

The characteristic timescale for the evolution of the residual eccentricity is on the order of τ_e, that is, much longer than the correlation timescale of the fluctuations of the quadrupole moment Q(t) that determine the function f(t), because the quadrupole fluctuations occur on the timescale of the convec-tive motions. In other words, we can regard f(t) as being a delta-correlated process, that is, $$\langle f\left(t\prime \right)f\left(t″\right)\rangle =D\delta \left(t″-t\prime \right),$$(31)

where D is a slowly varying function of the time and δ is the Dirac delta function. By substituting Eq. (31) into Eq. (30), we find $$\langle e^2\left(t\right)\rangle ={\left[\zeta \left(t\right)\right]}^2{\displaystyle {\int }_0^t{\displaystyle {\int }_0^{t}D\delta \left(t″-t\prime \right){\left[\zeta \left(t\prime \right)\right]}^{-1}{\left[\zeta \left(t″\right)\right]}^{-1}\text{\hspace{0.17em}d}t\prime \text{d}t″}}$$(32) $$={\left[\zeta \left(t\right)\right]}^2{\displaystyle {\int }_0^{t}D\left(t\prime \right){\left[\zeta \left(t\prime \right)\right]}^{-2}\text{\hspace{0.17em}d}t\prime },$$(33)

where we have explicitly introduced the slow time dependence of D that is relevant on the timescales over which 〈e²〉 evolves.

We can evaluate D(t) considering a much longer time interval than the correlation timescale of the stochastic function f(t), yet short enough that the parameters a and n of the unperturbed orbit and the mass of the star, m_s, are almost constant, so that K, τ_e, and the slowly varying D(t) can be regarded as constant. In this case, 〈e²〉 can be obtained by a simple integration from Eq. (33) as $$\langle e^2\rangle =\frac{1}{2}D{\tau }_{\text{e}},$$(34)

for t ≫ τ_e/2. On the other hand, 〈e²〉 can be obtained from the power spectrum of e that can be computed by taking the Fourier transform of Eq. (25).

We introduce the Fourier transform of a function g(t) of the time t as $$\tilde{g}\left(\omega \right)\equiv {\displaystyle {\int }_{-\infty }^{\infty }g\left(t\right)\exp \left(-i\omega t\right)\text{\hspace{0.17em}d}t,}$$(35)

where ω is the frequency and $i=\sqrt{-1}$ the imaginary unit. The inverse transform is $$g\left(t\right)=\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }\tilde{g}\left(\omega \right)\exp \left(i\omega t\right)\text{\hspace{0.17em}d}\omega .}$$(36)

Taking the Fourier transform of both sides of Eq. (25), we have $$\tilde{e}\left(\omega \right)=\frac{\tilde{f}\left(\omega \right)\left({\tau }_{\text{e}}^{-1}-i\omega \right)}{{\tau }_{\text{e}}^{-2}+{\omega }^2}.$$(37)

The power spectrum of e is the Fourier transform of its autocorrelation function R_e(τ) ≡ 〈e(t + τ)e(t)〉, where τ is a time lag, so that $$\langle e^2\rangle =R_{\text{e}}\left(0\right)=\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }{P}_{\text{e}}}\left(\omega \right)\text{\hspace{0.17em}d}\omega ,$$(38)

where $P_{\text{e}}\left(\omega \right)=\tilde{e}\left(\omega \right){\tilde{e}}^{*}\left(\omega \right)$ is the power spectrum of the eccentricity and the asterisk indicates complex conjugation. Making use of Eq. (37), Eq. (38) becomes $$\langle e^2\rangle =\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }\frac{P_{\text{f}}\omega }{{\tau }_{\text{e}}^{-2}+{\omega }^2}\text{d}\omega ,}$$(39)

where $P_{\text{f}}\left(\omega \right)=\tilde{f}\left(\omega \right){\tilde{f}}^{*}\left(\omega \right)$ is the power spectrum of the stochastic function f(t). As the tidal damping timescale τ_e is very long, the integrand in Eq. (39) is very large for ω = 0 and decreases sharply when |ω| > 0. In other words, it can be well approximated as $$\langle e^2\rangle \simeq \frac{P_{\text{f}}\left(0\right)}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }\frac{\text{d}\omega }{{\tau }_{\text{e}}^{-2}+{\omega }^2}=\frac{1}{2}{P}_{\text{f}}\left(0\right){\tau }_{\text{e}}.}$$(40)

To evaluate P_f(0), we use Eq. (26), the Euler formula for the sine, and Eq. (35) to compute $$\tilde{f}\left(0\right)={\displaystyle {\int }_{-\infty }^{\infty }-KQ\left(t\right)\frac{1}{2i}\left[\exp \left(int\right)-\exp \left(-int\right)\text{\hspace{0.17em}d}t\right]}$$(41) $$=\frac{K}{2i}\left[\tilde{Q}\left(n\right)-\tilde{Q}\left(-n\right)\right]=K\Im \left\{\tilde{Q}\left(n\right)\right\},$$(42)

where $\tilde{Q}$ is the Fourier transform of Q(t), which, being a real function, verifies ${\tilde{Q}}^{*}\left(n\right)=\tilde{Q}\left(-n\right)$. Given that Q(t) is a stochastic variable, ${\left[\Re \tilde{Q}\left(\omega \right)\right]}^2={\left[𝔍\tilde{Q}\left(\omega \right)\right]}^2=P_{\text{Q}}\left(\omega \right)/2$ in a statistical sense, where P_Q(ω) is the power spectrum of the fluctuating quadrupole moment Q(t). By applying this result, we have $$\langle e^2\rangle =\frac{1}{4}{K}^2{\tau }_{\text{e}}{P}_{\text{Q}}\left(n\right)=\frac{81}{16}\frac{n^2{\tau }_{\text{e}}{P}_{\text{Q}}\left(n\right)}{m_{\text{T}}^2{a}^4},$$(43)

which can be compared with Eq. (34) to find $$D\left(t\right)=\frac{81}{8}\frac{n^2{P}_{\text{Q}}\left(n\right)}{m_{\text{T}}^2{a}^4}.$$(44)

Substituting Eq. (44) into Eq. (33), we find the equation for the evolution of the expectation value of the square of the residual eccentricity as $$\langle e^2\left(t\right)\rangle =\frac{81}{8}{\left[\zeta \left(t\right)\right]}^2{\displaystyle {\int }_0^t\frac{n^2{P}_{\text{Q}}\left(n\right)}{m_{\text{T}}^2{a}^4}{\left[\zeta \left(t\prime \right)\right]}^{-2}\text{d}}t\prime ,$$(45)

where we have explicitly indicated the time dependence only for the function ζ given by Eq. (29), even if all the quantities appearing in the integrand are functions of the time. An alternative way to evaluate 〈e²〉 – according to the method proposed by Phinney (1992) – is presented in Appendix A.

The residual eccentricity 〈e²〉 is almost independent of the mass of the planet, because m_p ≪ m_s ≃ m_T and it enters into Eq. (45) only through the function ζ(t) via the timescale τ_e (cf. Eq. (29)). Considering the expression of τ_e given by Eq. (7), the factors with the planet mass cancel out in Eq. (45) because m_p is a constant that can be taken out of the integration. Moreover, the independence of the residual eccentricity on the planet mass is also true for the more general Eq. (30), provided that f(t) depends only on stellar quantities.

The statistical distribution of the residual eccentricity can be computed from the theory of stochastic processes (e.g., Phinney 1992; Lanza & Rodonò 2001). Nevertheless, a straightforward derivation exploits a method used to compute the Maxwell-Boltzmann distribution of molecular velocities in an ideal gas, as illustrated in Appendix B. Using p_e(e) de to indicate the probability that the residual eccentricity falls into the interval [e, e + de], the result is $$p_{\text{e}}\left(e\right)=\sqrt{\frac{2}{\pi \langle e^2\rangle }}\exp \left(-\frac{e^2}{2\langle e^2\rangle }\right).$$(46)

Therefore, the probability P(e ≤ e₀) that the residual eccentricity e is smaller than or equal to a given value e₀ can be expressed in terms of the error function as $$P\left(e\le e_0\right)=\text{erf}\left(\frac{e_0}{\sqrt{2\langle e^2\rangle }}\right),$$(47)

where we define the error function as $$\text{erf}\left(z\right)\equiv \frac{2}{\sqrt{\pi }}{\displaystyle {\int }_0^z\exp \left(-{\xi }^2\right)\text{d}\xi ,}$$(48)

so that lim_z→∞ erf(z) = 1.

2.3.2 Mean longitude at the epoch

The formal solution of Eq. (22) is $$ϵ\left(t\right)=-2K{\displaystyle {\int }_0^{t}Q\left(t\prime \right)\text{\hspace{0.17em}d}t\prime ,}$$(49)

assuming the initial condition ϵ(0) = 0. Given that the quadrupole moment Q is a stochastic variable with 〈Q〉 = 0, this must also mean that 〈ϵ〉 = 0, but 〈ϵ²〉 is different from zero and is a function of the time t because ϵ performs a Brownian motion around the point ϵ = 0. Specifically, the ensemble mean of ϵ² is $$⟨{ϵ}^2⟩=4K^2{\displaystyle {\int }_0^t{\displaystyle {\int }_0^t\langle Q\left(t\prime \right)Q\left(t″\right)\rangle \text{d}t\prime \text{d}t″.}}$$(50)

Considering much longer time intervals than the correlation time of the quadrupole moment fluctuations, we can assume a delta-correlated process with 〈Q(t′)Q(t″)〉 = Eδ(t″ − t′), finding 〈ϵ²〉 = 4K² Et as expected in the case of a Brownian motion.

To find the value of E, we use an approach based on the theory of the Brownian motion and consider that Q at any given time is the sum of the contributions coming from the different layers inside the convection zone of the star. Such contributions are uncorrelated with each other, while each of them has an autocorrelation timescale equal to the convective turnover time τ_c(r) at the radius r of the layer itself. We indicate the layers inside the stellar convection zone with an index j and assume that the layer centered at radius r_j contributes a variation ∆ϵ_j given by $$\text{Δ}{ϵ}_j\sim Q_j\left(t\right){\tau }_{\text{c}}\left(r_j\right),$$(51)

where Q_j(t) is the contribution to the quadrupole moment fluctuation at the time t coming from the j-th layer. We introduce the probability density function p_j(∆ϵ_j, t) giving the probability p_j(∆ϵ_j, t) d(∆ϵ_j) of having a fluctuation in the interval [∆ϵ_j, ∆ϵ_j + d(∆ϵ_j)] at the time t. By definition, $${\displaystyle {\int }_{-\infty }^{\infty }{p}_j\left(\text{Δ}{ϵ}_j,t\right)d\left(\text{Δ}{ϵ}_j\right)=1}.$$(52)

The probability ${𝒫}_j$ of having a fluctuation ϵ_j at the time t + dt as a consequence of the fluctuations of the quadrupole moment contribution of the jth layer is $${𝒫}_j\left({ϵ}_j,t+\text{d}t\right)={\displaystyle {\int }_{-\infty }^{\infty }{𝒫}_j\left({ϵ}_j-\text{Δ}{ϵ}_j,t\right)p_j\left(\text{Δ}{ϵ}_j,t\right)\text{d}\left(\text{Δ}{ϵ}_j\right)}$$(53) $$\begin{array}{l}={\displaystyle {\int }_{-\infty }^{\infty }[{𝒫}_j\left({ϵ}_j,t\right)p_j\left(\text{Δ}{ϵ}_j,t\right)-\frac{\partial {𝒫}_j}{\partial {ϵ}_j}\text{Δ}{ϵ}_j{p}_j\left(\text{Δ}{ϵ}_j,t\right)}\hfill \\ \hfill +\frac{1}{2}\frac{{\partial }^2{𝒫}_j}{\partial {ϵ}_j^2}{\left(\text{Δ}{ϵ}_j\right)}^2{p}_j\left(\text{Δ}{ϵ}_j,t\right)+\dots ]\text{d}\left(\text{Δ}{ϵ}_j\right),\hfill \end{array}$$(54)

where we develop the integrand function ${𝒫}_j$ in a Taylor series, retaining only the terms up to the second order, and consider that ∆ϵ_j is the variation in ϵ_j occurring along the small time step dt, whose probability distribution is given by p_j(∆ϵ_j, t) at the time t.

The integral of the first term in square brackets on the right-hand side of Eq. (54) is ${𝒫}_j\left({ϵ}_j,t\right)$ thanks to the normalization of the probability density p_j(∆ϵ_j, t) given by Eq. (52), while the integral of the second term vanishes because p_j(∆ϵ_j, t) = p_j(−∆ϵ_j, t). On the other hand, we can develop ${𝒫}_j$ in a time series and retain only the first-order term as $${𝒫}_j\left({ϵ}_j,t+\text{d}t\right)={𝒫}_j\left({ϵ}_j,t\right)+\frac{\partial {𝒫}_j}{\partial t}{\tau }_c\left(r_j\right)+\dots ,$$(55)

where we assume the small time increment dt to be equal to the autocorrelation time of the local quadrupole moment fluctuations, τ_c(r_j), in order to apply Eq. (51) in the later developments. By comparing Eqs. (54) and (55) and taking into account Eqs. (51) and (52), we find $$\frac{\partial {𝒫}_j}{\partial t}={𝒟}_j\frac{{\partial }^2{𝒫}_j}{\partial {ϵ}_j^2},$$(56)

with $${𝒟}_j\equiv \frac{1}{2}{\displaystyle {\int }_{-\infty }^{\infty }4K^2{Q}_j^2{\tau }_{\text{c}}\left(r_j\right)p_j\left(Q_j,t\right)\text{d}{Q}_j=2K^2\langle Q_j^2\rangle {\tau }_{\text{c}}\left(r_j\right),}$$(57)

where we make use of Eq. (51) to transform the integration variable ∆ϵ_j into the corresponding quadrupole moment fluctuation Q_j whose probability density function is p_j(Q_j, t). Equation (56) is a diffusion equation whose solution is $${𝒫}_j\left({ϵ}_j,t\right)=\frac{1}{\sqrt{2\pi {𝒟}_{j}t}}\exp \left(-\frac{{ϵ}_j^2}{2{𝒟}_{j}t}\right).$$(58)

Therefore, the variance of ϵ_j is $\langle {ϵ}_j^2\rangle ={𝒟}_{j}t$ and increases linearly in time as expected in the case of a system performing a Brow-nian motion. The variance of ϵ can be found by summing all the uncorrelated contributions coming from the single layers inside the stellar convection zone. This yields $$\langle {ϵ}^2\rangle ={\displaystyle \sum _j\langle {ϵ}_j^2\rangle =2K^{2}t{\displaystyle {\int }_{r_{\text{b}}}^{R_{\text{s}}}\langle Q^2\left(r\right)\rangle {\tau }_{\text{c}}\left(r\right)\text{d}r,}}$$(59)

where r_b is the radius at the base of the stellar convection zone and R_s is the stellar radius, and we have dropped the index j that specifies the layer inside the star because the sum over the layers has been substituted by the integration We do not consider the time dependence of K and 〈Q²(r)〉 in Eq. (59) because it is relevant only on evolutionary timescales, that is, on much longer time intervals than any observational interval over which 〈ϵ²〉 can be measured.

The fluctuations in the time of midtransit can be derived from Eq. (24) and are $$\langle O-C\rangle =0,$$(60) $$\langle {\left(O-C\right)}^2\rangle =\frac{\langle {ϵ}^2\rangle }{2\pi }{P}_{\text{orb}},$$(61)

where Eq. (60) follows from 〈ϵ〉 = 0, while 〈ϵ²〉 in Eq. (61) is given by Eq. (59). Therefore, O – C makes a Brownian motion around zero with a standard deviation that increases linearly in time. This Brownian motion can be measured provided that we have observations extending over a sufficiently long interval of time, as we consider in Sect. 4.

2.4 Power spectrum of the quadrupole moment fluctuations

To compute the residual eccentricity according to Eq. (45), we need to determine the power spectrum of the time-dependent quadrupole moment Q(t). From the power spectrum of the local contribution to the quadrupole moment fluctuations at radius r, the mean squared value 〈Q²(r)〉 appearing in Eq. (59) follows, which allows us to compute the variance of the mean longitude at the epoch.

First we consider the power spectrum of the total quadrupole moment of the star, Q(t), and make use of the property that the power spectrum is the Fourier transform of the autocorrelation function R_Q(τ) of Q(t). The autocorrelation function is the expectation value R_Q(τ) ≡ 〈Q(t + τ)Q(t)〉, where τ is a time lag, and we make use of the definition (35) of the Fourier transform to obtain $$P_Q\left(\omega \right)={\displaystyle {\int }_{-\infty }^{\infty }{R}_Q\left(\tau \right)\exp \left(-i\omega \tau \right)\text{\hspace{0.17em}d}\tau \text{.}}$$(62)

The quadrupole Q can be expressed as (cf. Eq. (15)) $$Q=-\frac{1}{2}{Q}_{zz}=-\frac{1}{2}\left(I_{zz}-\frac{1}{3}\text{Tr}I\right)=-\frac{1}{3}{I^{\prime }}_{zz}+\frac{1}{6}\left({I^{\prime }}_{xx}+{I^{\prime }}_{yy}\right).$$(63)

The fluctuations ${I^{\prime }}_{xx}$, ${I^{\prime }}_{yy}$ and ${I^{\prime }}_{zz}$ are uncorrelated with each other and have the same autocorrelation function because of the spherical symmetry of the turbulent velocity field, that is, $R_{{I^{\prime }}_{xx}}\left(\tau \right)=R_{{I^{\prime }}_{yy}}\left(\tau \right)=R_{{I^{\prime }}_{zz}}\left(\tau \right)$. This implies that $$R_Q\left(\tau \right)=\frac{1}{6}{R}_{{I^{\prime }}_{xx}}\left(\tau \right).$$(64)

Therefore, we are left with the computation of the autocorrelation function of the fluctuations of the principal moment of inertia ${I^{\prime }}_{xx}$. Considering the isotropy of the fluctuations in our model $$I{\prime }_{xx}=\frac{1}{3}\left(I{\prime }_{xx}+I{\prime }_{yy}+I{\prime }_{zz}\right)=\frac{1}{3}{\displaystyle {\int }_V{\rho }^{\prime }}\left(r,t\right)r^2\text{d}V,$$(65)

where r is the position vector, the modulus r of which is the radial distance from the center of the star. The autocorrelation of ${I^{\prime }}_{xx}$ becomes $$\begin{array}{l}{R}_{I{\prime }_{xx}}\left(\tau \right)=\langle I{\prime }_{xx}\left(t+\tau \right)I{\prime }_{xx}\left(t\right)\rangle \hfill \\ =\frac{1}{9}\langle {\displaystyle {\int }_V{\rho }^{\prime }}\left(r,t+\tau \right)r^2\text{d}V\times {\displaystyle {\int }_V{\rho }^{\prime }}\left(r\prime ,t\right)r{\prime }^2\text{d}V\prime \rangle \hfill \\ =\frac{1}{9}{\displaystyle {\int }_V{\displaystyle {\int }_V\langle {\rho }^{\prime }\left(r,t+\tau \right)\rho \prime \left(r\prime ,t\right)\rangle }{r}^{2}r{\prime }^2\text{d}V\text{\hspace{0.17em}d}V\prime .}\hfill \end{array}$$(66)

The autocorrelation function of the density fluctuations can be written as the product of a spatial autocorrelation and a temporal autocorrelation. In the framework of the mixing-length theory, that is, a local theory of convection, the density fluctuations at a given radius r are spatially correlated over a volume of order ℓ³(r), where ℓ(r) is the mixing length at radius r, while their correlation time is τ_c(r) = ℓ(r)/υ_c(r), where υ_c(r) is the convective velocity at radius r as provided by the mixing-length (or any, more sophisticated) convection model.

If we assume that ℓ ≪ r in the convective envelope and that the fluctuations outside each convective cell of volume ℓ³ are uncorrelated with each other, we can simply integrate all the local autocorrelations of the density fluctuations and obtain² $$\begin{array}{l}{R}_{I{\prime }_{xx}}\left(\tau \right)\simeq \frac{1}{9}{{\displaystyle {\int }_V\left[\rho \prime \left(r,0\right)\right]}}^2{\ell }^3\left(r\right)r^4\exp \left[-\left|\tau \right|/{\tau }_c\left(r\right)\right]\text{d}V\hfill \\ =\frac{4\pi }{9}{{\displaystyle {\int }_{r_{\text{b}}}^{R_{\text{s}}}\left[\rho \prime \left(r\right)\right]}}^2{\ell }^3\left(r\right)r^6\exp \left[-\left|\tau \right|/{\tau }_c\left(r\right)\right]\text{d}r,\hfill \end{array}$$(67)

where r_b is the radius at the lower boundary of the stellar convection zone, R_s the radius of the star, and we assume that the local fluctuations do not explicitly depend on the time and are spherically symmetric, on the average.

The local density fluctuations can be estimated by means of the mixing-length theory (see Ch. 7 of Kippenhahn et al. 2012) by equating the average work per unit volume done by the buoyancy force along the path of a convective element (equal to the mixing length ℓ(r)) to the kinetic energy density in the convective motions, $$\frac{1}{2}{\rho }_0\left(r\right){\upsilon }_{\text{c}}^2\left(r\right)\sim \frac{1}{2}g\left(r\right)\left|\rho \prime \left(r\right)\right|\ell \left(r\right),$$(68)

where ɡ(r) = Gm(r)r⁻² is the acceleration of gravity with m(r) being the mass of the star inside the radius r. Therefore, we find $$\left|\rho \prime \left(r\right)\right|\sim \frac{{\rho }_0\left(r\right){\upsilon }_{\text{c}}^2\left(r\right)}{g\left(r\right)\ell \left(r\right)}.$$(69)

Substituting Eq. (69) into Eq. (67), we find $$R_{I{\prime }_{xx}}\left(\tau \right)\simeq \frac{4\pi }{9}{{\displaystyle {\int }_{r_{\text{b}}}^{R_{\text{s}}}\left[\frac{{\rho }_0\left(r\right){\upsilon }_{\text{c}}^2\left(r\right)}{g\left(r\right)\ell \left(r\right)}\right]}}^2{\ell }^3\left(r\right)r^6\exp \left[-\left|\tau \right|/{\tau }_c\left(r\right)\right]\text{d}r,$$(70)

where the quantities ρ₀(r), υ_c(r), ℓ(r), τ_c(r), and ɡ(r) follow from a model of the stellar structure at the given time t.

To facilitate the computation of the power spectrum of ${I^{\prime }}_{xx}$, we compute the Fourier transform of the exponential decorrela-tion function of the time that appears in Eq. (70): $$E_{\text{D}}\left(\tau \right)\equiv \exp \left(-\left|\tau \right|/{\tau }_0\right).$$(71)

We have $$\begin{array}{l}{\tilde{E}}_{\text{D}}\left(\omega \right)\equiv {\displaystyle {\int }_{-\infty }^0\exp \left(-\left|\tau \right|/{\tau }_0\right)}\exp \left(-i\omega \tau \right)\text{d}\tau \hfill \\ ={\displaystyle {\int }_{-\infty }^0\exp \left(\tau /{\tau }_0\right)}\exp \left(-i\omega \tau \right)\text{d}\tau +{\displaystyle {\int }_0^{\infty }\exp \left(-\tau /{\tau }_0\right)}\exp \left(-i\omega \tau \right)\text{d}\tau \hfill \\ ={\displaystyle {\int }_{-\infty }^0\exp \left[\left(\frac{1}{{\tau }_0}-i\omega \right)\tau \right]\text{d}\tau +{\displaystyle {\int }_0^{\infty }\exp \left[-\left(\frac{1}{{\tau }_0}+i\omega \right)\tau \right]\text{d}\tau }}\hfill \\ =\frac{1}{\left({\tau }_0^{-1}-i\omega \right)}+\frac{1}{\left({\tau }_0^{-1}+i\omega \right)}=\frac{2{\tau }_0^{-1}}{{\tau }_0^{-2}+{\omega }^2}=\frac{2{\tau }_0}{1+{\tau }_0^2{\omega }^2}.\hfill \end{array}$$(72)

The power spectrum of ${I^{\prime }}_{xx}\left(t\right)$ follows from the Fourier transform of Eq. (70). It can be immediately computed because the Fourier transform implies an integral over τ that can be taken out of the integral over the radius r giving, thanks to Eq. (72), $$P_{I{\prime }_{xx}}\left(\omega \right)\sim \frac{8\pi }{9}{{\displaystyle {\int }_{r_{\text{b}}}^{R_{\text{s}}}\left[\frac{{\rho }_0\left(r\right){\upsilon }_{\text{c}}^2\left(r\right)}{g\left(r\right)\ell \left(r\right)}\right]}}^2{\ell }^3\left(r\right)r^6\frac{{\tau }_{\text{c}}\left(r\right)}{{\tau }_{\text{c}}^2\left(r\right){\omega }^2+1}\text{d}r.$$(73)

The power spectrum of Q follows from the Fourier transform of both sides of Eq. (64), giving $$P_Q\left(\omega \right)=\frac{1}{6}{P}_{I{\prime }_{xx}}\left(\omega \right).$$(74)

Specifically, P_Q(n), which appears in Eq. (45), is given by $$P_Q\left(n\right)\sim \frac{4\pi }{27}{{\displaystyle {\int }_{r_{\text{b}}}^R\left[\frac{{\rho }_0\left(r\right){\upsilon }_{\text{c}}^2\left(r\right)}{g\left(r\right)\ell \left(r\right)}\right]}}^2{\ell }^3\left(r\right)r^6\frac{{\tau }_{\text{c}}\left(r\right)}{{\tau }_{\text{c}}^2\left(r\right)n^2+1}\text{d}r.$$(75)

To evaluate the integral that appears in Eq. (59), we compute the mean variance of the local contribution to the quadrupole moment fluctuations at radius r from Eq. (74) and the integrand in Eq. (73), that is, the power spectrum of the local fluctuations of the moment of inertia. We have $$\langle Q^2\left(r\right)\rangle {\tau }_c\left(r\right)\rangle =\frac{2}{27}{\left[\frac{{\rho }_0\left(r\right){\upsilon }_c^2\left(r\right)}{g\left(r\right)\ell \left(r\right)}\right]}^2{\ell }^3\left(r\right)r^6{\displaystyle {\int }_{-\infty }^{\infty }\frac{{\tau }_{\text{c}}^2\left(r\right)}{{\tau }_{\text{c}}^2\left(r\right){\omega }^2+1}}\text{\hspace{0.17em}d}\omega ,$$(76)

which can be substituted into Eq. (59) after computing the integral over ω to give $$\langle {ϵ}^2\rangle =\frac{3\pi n^2}{m_{\text{T}}^2{a}^4}\left\{{{\displaystyle {\int }_{r_{\text{b}}}^{R_{\text{s}}}\left[\frac{{\rho }_0\left(r\right){\upsilon }_c^2\left(r\right)}{g\left(r\right)\ell \left(r\right)}\right]}}^2{\ell }^3\left(r\right){\tau }_{\text{c}}\left(r\right)r^6\text{d}r,\right\}t,$$(77)

where we make use of Eq. (27) to express K appearing in Eq. (59).

3 Stellar structure and evolution

To apply the above model to an Earth-like planet orbiting a Sunlike star, we first introduce a model for the internal structure and evolution of the host star. We used the r21.12.1 release of the MESA stellar evolution code (Paxton et al. 2011, 2013, 2015, 2018, 2019) to calculate the structure and evolution of a one solar-mass star from the pre-main sequence phase to the tip of the AGB.

We considered the standard inputs and recipes recommended by the MESA team for the equation of state, nuclear reaction rates, and screening factors (Paxton et al. 2019), as well as for the radiative and conductive opacities (Paxton et al. 2011). Regarding the mixing-length theory of convection, we adopted the formulation by Henyey et al. (1965). The Ledoux's (1947) criterion was used to assess convective stability.

The occurrence of core overshooting during the core helium burning phase (CHeB) was suggested by Bossini et al. (2015) to satisfy both asteroseismic constraints on the CHeB phase and the observed luminosity at the AGB bump. Then, considering the observed constraints on the AGB bump of solar-mass stars of solar metallicity, Dréau et al. (2022) estimated that the ratio of the extent of overshooting to the pressure scale height is in the range α_ov,CHeB = 0.25–0.50, following a step scheme (Maeder 1975). In our model, we therefore take into account the moderate amount of convective core overshooting, α_ov,CHeB = 0.25. The temperature gradient in the overshooting region is taken to be radiative. During CHeB, the convective core grows with time; the convective boundaries must therefore be correctly located to allow the emergence of semi-convective regions, which is a delicate procedure. To cope with these difficulties, we used the new convective pre-mixing (CPM) scheme recommended by the MESA team (Paxton et al. 2019). Furthermore, we considered step overshooting from the bottom of the convective envelope, with a value α_ov,env = 0.3. This value was obtained by Khan et al. (2018) on the basis of constraints on the position of the RGB bump of ~3000 stars. Moreover, concerning mixing, we did not take into account microscopic diffusion, nor rotational mixing, nor thermohaline mixing.

We took mass loss into account both on the RGB and on the AGB. From the base of the RGB up to the end of the core helium burning phase, we used Reimers’ prescription (Reimers 1975) with the moderate, maximum value η_R = 0.3 of the mass-loss parameter, inferred by Miglio et al. (2012, 2021) from asteroseis-mic constraints on red giants in the Kepler field. On the AGB, we used Blocker's (1995) prescription with η_B = 0.02 as prescribed by Ventura et al. (2020). As for the atmosphere boundary condition, we used an Eddington’s gray T(τ) relation with τ being the optical depth here.

We adopted the GS98 solar mixture (Grevesse & Sauval 1998) as a reference. It corresponds to a present solar photo-spheric metal Z to hydrogen X mass fraction ratio (Z/X)_⊙ = 0.0229. The solar model calibration, that is, the requirement that the solar model reaches the solar radius and luminosity³ at the solar age, provides the solar initial helium abundance in mass fraction Y₀ = 0.2602 and a mixing-length parameter value of α_MLT,⊙ = 1.950.

The late evolution of the stellar radius, total mass, mass of the convection zone, and luminosity is shown in Fig. 1. At the tip of the RGB, the maximum radius is 181.8 R_⊙, while the mass and luminosity are 0.896 M_⊙ and 2776 L_⊙, respectively. On the other hand, the maximum radius at the tip of the AGB is 163.5 R_⊙ with a mass of 0.886 M_⊙ and a luminosity of ~2477 L_⊙, respectively.

In Fig. 2, we plot the mass-loss parameter Ψ_ml introduced in Eq. (4) for a planet initially at a distance of 1.5 au from the star, focusing on the late stages of the stellar evolution when Ψ_ml reaches its maximal values. We see that Ψ_ml is always smaller than 2 × 10⁻⁸, even at the tip of the RGB or in the AGB phase, when the mass-loss rates are the largest. Therefore, we conclude that the adiabatic mass-loss approximation is very well verified in our case and our model equations for the evolution of the semimajor axis and the eccentricity given in Sect. 2.1 hold. We note that, in our stellar evolution modeling, we always adopt an isotropic stellar mass loss because our stellar structure model is unidimensional.

We explored several reasonable options for input parameters in MESA, such as the use of the AGSS09 solar mixture of Asplund et al. (2009). We find minor differences of less than 5% in the RGB-tip radius. We also point out that the recent BaSTI-IAC one-solar-mass model of Hidalgo et al. (2018), which includes mass loss, predicts RGB-tip radii of only 3% smaller than ours.

Our stellar parameters differ from those of the evolution model adopted by Schröder & Smith (2008); in particular, these authors find a larger radius of the Sun at the tip of the RGB phase of 256 R_⊙ and a lower mass of 0.668 M_⊙. We implemented the Schröder & Cuntz (2005) mass-loss formalism used by Schröder & Smith (2008) in MESA, both on the RGB and AGB branches. The mass loss is higher and at the RGB tip, the mass, luminosity, and effective temperature are lower by 21%, 1.3%, and 95 K respectively, leading to an RGB tip radius of ~191 R_⊙, which is ~5% higher than our reference value. Furthermore, we point out that, in order to better fit the observed position of evolved giants in the Hertzsprung-Russell diagram, Schröder & Smith (2008) progressively and linearly decreased the value of the mixing-length parameter in their models starting from their solar calibrated value α_MLT,⊙ = 2.0 at log ɡ = 1.94, which is at about half way from the base of the RGB to the RGB tip, down to a value of 1.67 at log ɡ ≈ 0.0 at the RGB tip. We verified that if we mimic this Schröder & Smith (2008) recipe by taking a value of α_MLT,⊙ ≈ α_MLT,⊙ − 0.30 when ascending the RGB, we get higher values of the tip radius and smaller values of the tip T_eff, which is similar to what these authors found.

Fig. 1 Late evolution of the main parameters of our Sun-like stellar model. Top panel: stellar radius (red solid line) vs. the time. Middle panel: stellar mass (red solid line) vs. the time; the mass of the convective zone is overplotted as a green dashed line. Bottom panel: stellar luminosity vs. the time (red solid line). We note the discontinuities in the radius and luminosity when the He flash occurs at the tip of the RGB evolution.

Fig. 2 Stellar mass-loss parameter Ψml vs. the stellar age in the evolutionary phases characterized by a significant stellar mass loss rate.

4 Applications

We apply our model for the convective-induced residual eccentricity and transit-time variation to the case of a planet of one Earth mass orbiting a star of initial mass equal to that of the Sun, the structure and evolution of which have been computed as described in Sect. 3.

The top panel of Fig. 3 shows the evolution of the stellar radius according to our model and the corresponding change in the orbit semimajor axis for an Earth-mass planet at an initial separation of 1.0 au according to Eqs. (5) and (6), where we adopt λ₂ = 0.046. In our model, the maximum radius at the tip of the RGB phase is 29% smaller than in the Schröder & Smith (2008) model, but this is not enough to allow the Earth to avoid the engulfment, even without including the effect of the residual eccentricity. This happens because the lower mass-loss rate in our model is not capable of overcoming the effects of tides leading to the eventual decay of the Earth’s orbit.

The fate of the Earth after entering the envelope of the RGB Sun is to slowly spiral toward the center of the star where it is finally destroyed close to the core at a radius of r_d/R ~ 0.003 – where R is the stellar radius at the RGB tip – after a time interval on the order of ~3 × 10³ yr following its engulfment, as we show in Appendix C. The effects on the stellar structure, luminosity, and rotation are negligible. Even in the case of the engulfment of a giant planet, those effects are not expected to be large given the large luminosity of the giant star and the slow spiral in of the planet (cf. MacLeod et al. 2018; Yarza et al. 2022). However, in some systems, the added effects of successive planetary engulfments or of the additional energy inputs during the helium flash or the AGB thermal pulses may lead to the ejection of the common envelope and the formation of systems such as WD 1856+534, where a WD is accompanied by a planet with a mass of ≲ 14 Jupiter masses and an orbital period of 1.4 days (e.g., Lagos et al. 2021; Chamandy et al. 2021; Merlov et al. 2021).

Reverting to the evolution of our planet before its engulfment, any initial eccentricity of the orbit of the Earth is rapidly erased close to the tip of the RGB, because the tidal damping timescale becomes shorter than 0.1 Myr close to that evolution point, as shown by the second panel of Fig. 3, where τ_e has been computed with Eq. (7), adopting ß = 1. Including the effects of the residual eccentricity, the engulfment may occur slightly earlier, because the most probable value of the eccentricity becomes close to 0.01 during the final phase of the RGB ascent. Nevertheless, the difference is of very little practical relevance because the residual eccentricity becomes significant when the orbit of the planet is already on its final tidal decay path (cf. the third panel of Fig. 3). The variation in the time of midtransit with respect to a constant-period ephemeris is very small and reaches only a few seconds over a time baseline of 10 yr during a phase occurring only 1–2 Myr before the engulfment (cf. the bottom panel of Fig. 3). This phase is so short in comparison to the evolutionary timescale of the star, and the variation in the time of midtransit is so small, that it is extremely unlikely that it is observable.

The role of the residual eccentricity becomes relevant if we consider a planet on an initially wider orbit that allows it to escape engulfment. The top panel of Fig. 4 shows the evolution of the solar radius according to our model and the evolution of the orbital semimajor axis for an Earth-mass planet with an initial separation of 1.02 au, which is sufficient to avoid engulfment assuming only the tidal orbital decay; this is a consequence of the strong dependence of the tidal decay on the relative orbital separation as expressed by the factor (R_s/a)⁶ in Eq. (5). Now, when the star reaches the tip of the RGB, the orbit semimajor axis is 1044 times larger than the radius of the star. The sudden drop in the stellar radius by more than one order of magnitude produced by the internal structure changes associated with the helium flash puts the planet in safety, halting its tidal orbital decay. Any initial orbital eccentricity is effectively erased by tides because τ_e becomes shorter than ~1 Myr close to the RGB tip (cf. the second panel of Fig. 4). Nevertheless, the situation changes if we take into account the residual eccentricity, because 〈e²〉^1/2 = 2.65 × 10⁻² at the tip of the RGB evolution. The planet can avoid engulfment if the eccentricity of its orbit is smaller than 0.044; otherwise its periastron distance becomes smaller than the radius of the star. The probability of such a condition is given by Eq. (47) with e₀ = 0.044, and is 0.903. In other words, taking into account the residual eccentricity, there is a ~10% probability that the planet will be engulfed during its periastron passage close to the tip of the RGB.

If the planet is capable of escaping engulfment at the RGB tip, its residual eccentricity will be frozen until the end of the stellar evolution because the maximum of the stellar radius reached at the tip of the AGB phase is remarkably smaller than the maximum reached at the RGB tip in our model and does not affect the value of 〈e²〉^1/2 in any significant way. In this case, the eccentricity of the orbit of the planet is maintained over the WD phase of the system and may play a relevant role in the pollution of the remnant WD by residual planetesimals in the system (Frewen & Hansen 2014). The variation in the time of midtransit with respect to a constant-period ephemeris over a time baseline of 10 yr is plotted in the bottom panel of Fig. 4 and reaches a maximum of 1.2 s, making it extremely difficult to detect.

The maximum residual eccentricity is reached in our model close to the tip of the RGB and remains almost constant throughout the final evolution phases, provided that the planet can avoid engulfment. This is a consequence of the fact that the maximum stellar radius is reached at the tip of the RGB phase followed by a sudden drop at the helium flash, while the maximum radius reached during the tip of the AGB phase is significantly smaller in our model. Therefore, the damping of the residual eccentricity becomes negligible after the RGB tip − cf. the dependence of the damping timescale on (R_s/a) in Eq. (7) − while the excitation of the eccentricity by the quadrupole moment fluctuations is negligible because they also decrease rapidly with decreasing (Rs/a).

Figure 5 shows the residual eccentricity 〈e²〉^1/2 (top panel) and the minimum of the eccentricity damping timescale τ_e (lower panel) that are reached at the tip of the RGB phase for different initial values of the orbital semimajor axis. As already noted in Sect. 2.3, 〈e²〉^1/2 is almost independent of the mass of the planet. We see that the residual eccentricity at the tip of the RGB and in later evolutionary phases is greater than 0.01 for an initial semimajor axis up to about 1.4 au and decreases rapidly with increasing semimajor axis because of the rapid decrease in quadrupole potential. Any initial orbital eccentricity is rapidly erased on timescales shorter than ~10 Myr for an initial semimajor axis of < 1.2 au, leaving only the residual eccentricity in those final phases of the stellar evolution up to the WD phase. On the other hand, for planets with an initial semimajor axis of a ≳ 1.2–1.3 au, the minimum tidal damping timescale increases remarkably and any initial eccentricity is not erased during the red giant phase of the stellar evolution.

In planetary systems with a central star in the WD stage of its evolution, the presence of an eccentric planet is a necessary condition for the perturbation of the orbits of the asteroids leading to their accretion onto the WD. Extensive analytical studies (e.g., Antoniadou & Veras 2016, 2019) and numerical simulations (e.g., Frewen & Hansen 2014; Veras et al. 2021) showed that a planet on a circular orbit is not capable of inducing accretion of the asteroids inside the Roche sphere of the WD, while this is possible in the case of an eccentric planetary orbit, even for eccentricities as small as a few hundredths. A detailed exploration of the role played by circularization of the planetary orbit and residual eccentricity occurring close to the tip of the RGB in WD accretion is beyond the scope of the present work, but the results summarized in Fig. 5 can provide a useful starting point for such investigations.

Fig. 3 Late evolution of the radius of our Sun-like stellar model and orbital parameters of an Earth-mass planet with an initial semimajor axis of 1.0 au. Top panel: stellar radius (red solid line) and orbital semi-major axis (green solid line) vs. time. Second panel: timescale for the damping of the orbital eccentricity τe vs. time (red solid line) computed according to Eq. (7) with ß = 1. Third panel: residual eccentricity 〈e2 〉1/2 vs. time (red solid line) as computed from Eq. (45). Bottom panel: observed minus calculated time of midtransit 〈(O − C)2〉1/2 vs. time (red solid line) computed according to Eqs. (61) and (77) over a time interval of 10 yr. The green line in the top panel and the plots in the second, third, and bottom panels are truncated at the time of the engulfment of the planet by the red giant star.

Fig. 4 Same as Fig. 3, but for an Earth-mass planet with an initial semi-major axis of 1.02 au, which is sufficient to escape engulfment at the tip of the RGB phase when only the tidal orbital decay is considered. Therefore, the plots are extended until the time when the star reaches the AGB tip.

Fig. 5 Relevant parameters at the tip of the RGB phase vs. the initial semimajor axis of the planetary orbit. Top panel: mean residual eccentricity 〈e2〉1/2; bottom panel: tidal decay timescale of the eccentricity τe.

5 Discussion and conclusions

We introduce a model to compute the residual eccentricity of the orbit of a planet around a star with an outer convective envelope. The eccentricity is maintained against the tidal damping by the random fluctuations of the stellar outer gravitation field, which are caused by the small density perturbations that drive convective motions. In the case of an Earth-like planet orbiting a Sun-like star, this eccentricity is negligible during the main sequence phase of the evolution of the star, but becomes relevant during its red giant phase, especially close to the tip of the RGB when the star becomes hundreds of times larger and thousands of times more luminous than on the main sequence, making the amplitude of convective density fluctuations remarkably greater.

The residual eccentricity is a random variable with a Gaussian probability density distribution. We find a maximum mean value of the eccentricity of 〈e²〉^1/2 ~ 0.026 using our model in the case of an Earth-like planet with an initial semimajor axis of 1.02 au, which is the minimum value required to have a probability of ~0.9 of escaping engulfment at the tip of the RGB phase.

Direct observation of the orbit of the planet in order to measure its residual eccentricity would require imaging with a contrast on the order of 10¹⁰ at an angular resolution of 1 microarcsec to detect an Earth-size planet around a red giant star at a distance of ~1000 pc. Such observations can become possible, in principle, by using the Sun as a gravitational lens (see Turyshev & Toth 2022, and references therein).

The introduction of the residual eccentricity demonstrates that the engulfment of an Earth-like planet with an initial orbital semimajor axis of 1.0 au by a Sun-like star during its RGB evolution is a stochastic process, which is in contrast to the findings of previous studies where it was considered to be a deterministic process ruled by the stellar evolution and the tidal decay of the orbit (cf. Schröder & Smith 2008).

The late orbital evolution scenario of an Earth-like planet critically depends on the stellar-mass-loss rate during the red giant and the asymptotic branch phases because it significantly affects the evolution of the stellar radius and its maximum value. This mass-loss rate is still poorly known and may depend on physical processes that are not solely under the control of stellar evolution. Specifically, Sabach & Soker (2018) speculate that the mass-loss rate can be increased by a binary companion or a massive close-by planet that produce a remarkable tidal interaction or are engulfed during the ascent of the RGB branch, leading to a faster stellar rotation that may enhance the mass loss. As our mass-loss rate parameterizations are calibrated with stars whose binarity status is not known, they suggest that they may suffer significant contamination by binarity effects, while the mass-loss rate for a star without massive companions, such as our Sun, would be significantly smaller. As a consequence, considering for example a reduction in the mass-loss rate by a factor of ~7, the maximum radius reached at the tip of the AGB becomes more than twice larger than that at the tip of the RGB, leading to the tidal engulfment of a telluric planet even if its initial orbital semimajor axis is as large as ~ 1.4–1.5 au.

If the planet escapes engulfment during the RGB and AGB phases, the residual eccentricity can play a remarkable role during the WD phase of a planetary system because it can excite the eccentricity of the orbits of residual planetesimals leading them to collide with the WD and pollute its atmosphere; even residual eccentricities as small as a few hundredths are sufficient for this (cf. Frewen & Hansen 2014). However, our analysis does not take into account the role of other possible planets in the evolution of the eccentricity of an Earth-like planet surviving the RGB and the AGB phases of its star.

Angular momentum exchanges among several planets can lead to an increase in the eccentricity of the orbit of an inner low-mass planet as shown, for example, in the angular momentum deficit model by Laskar (1997). We shall not discuss the possible role of other planets in a system because this is outside the scope of the present work, but we note that the residual eccentricity produced by convective fluctuations in a red giant can play a role in the exchanges of angular momentum between planets in a multi-planet system and should be considered in the model of their orbital evolution.

Another consequence of the fluctuating gravitational potential of a red giant are the small changes in the mean longitude at the epoch of a close-by orbiting planet. Considering a time baseline of 10 yr, we find deviations in the time of midtransit with respect to a constant-period ephemeris that do not exceed a few seconds. These are completely undetectable because planetary transits are extremely shallow in the case of a red giant star owing to its large radius. Observations in the core of chromospheric spectral lines have been proposed to improve the observability and the timing of transits across a red giant star ascending the RGB because its chromosphere has a much thinner radial extension than the stellar radius, which enhances the photometric depth of the transits (Assef et al. 2009). Nevertheless, this method is not useful in our case because of the extended chromospheres of stars close to the RGB tip and the small radius of Earth-like planets.

Acknowledgements

The authors are grateful to an anonymous Referee for a careful reading of the manuscript and several suggestions that helped them to substantially improve the presentation of their work. This investigation was begun during an internship of CS at INAF-Osservatorio Astrofísico di Catania held remotely in 2021. Research in the field of exoplanetary astronomy at Osservatorio Astrofísico di Catania, that belongs to the Italian National Institute for Astrophysics (INAF), is supported by the Italian Ministry for University and Research.

Appendix A Residual eccentricity according to the method of Phinney (1992)

Fig. A.1 Reference frame in the plane of the orbit to describe the relative motion of the planet P around the star and define the true anomaly ƒ and the relative separation of the barycenters of the star and the planet r = OP. The origin of the reference frame is at the barycenter О of the star, while the axis x0 points toward the periastron of the orbit.

In this Appendix, we derive the value of 〈e²〉 by revisiting the approach introduced by Phinney (1992), which is based on the equation for the epicyclic motion with respect to an unperturbed circular orbit. We derive the equations of the motion in Sect. A.1 following the Lagrangian formalism and solve the equation for the radial epicyclic motion in Sect. A.2 to find the mean squared residual eccentricity.

Appendix A.1 Lagrangian function of the star-planet system

To study the dynamics of our star-planet system, we apply the Lagrangian formalism. The Lagrangian ℒ of our system is defined as $$ℒ=𝒯-\text{Ψ},$$(A.1)

where 𝒯 is the kinetic energy of the system and Ψ its potential energy expressed as functions of the coordinates and their time derivatives in an inertial reference frame. We consider a reference frame with its origin at the barycenter Z of the star-planet system and whose x₀y₀ plane is the orbital plane of the system. Given that m_p ≪ m_s, the barycenter of the system Z almost perfectly coincides with the barycenter О of the star. We consider the star as nonrotating and the planet as a point mass m_p, and therefore we neglect the kinetic energy of their rotation.

The kinetic energy of the relative orbital motion of our system can be expressed as the energy of the motion of a body with 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 О of the star. Therefore, the expression of the kinetic energy of the system is $$𝒯=\frac{1}{2}m\left({\dot{r}}^2+r^2{\dot{f}}^2\right),$$(A.2)

where r is the relative separation between the two bodies and ƒ is the true anomaly of the relative orbital motion (see Fig. A.1); the dot over a variable indicates its time derivative.

The gravitational potential energy Ψ is given by (cf. Eq. 17) $$\text{Ψ}={\text{Φ}}_{m_{\text{p}}}=-\frac{G{m}_{\text{s}}{m}_{\text{p}}}{r}-\frac{3GQ\left(t\right)m_{\text{p}}}{2r^3},$$(A.3)

where Q is the axisymmetric component of the quadrupole moment of the star as defined in Sect. 2.2; because of the density fluctuations, it is a function of the time t.

The equations of motion of our system can be derived from its Lagrangian as $$\frac{d}{dt}\left(\frac{\partial ℒ}{\partial \dot{q}}\right)-\frac{\partial ℒ}{\partial q}=0,$$(A.4)

where t is the time and q = r, ƒ is any of the coordinates adopted to describe the motion of the system. In this way, we obtain the following equations of motion: $$m\ddot{r}-mr{\dot{f}}^2+\frac{G{m}_{\text{s}}{m}_{\text{p}}}{r^2}+\frac{9GQ\left(t\right)m_{\text{p}}}{2r^4}=0,$$(A.5) $$m\frac{d}{dt}\left(r^2\dot{f}\right)=0.$$(A.6)

Equation (A.6) can be immediately integrated to give the conservation of the orbital angular momentum J of the system: $$m{r}^2\dot{f}=J.$$(A.7)

Considering the definition of the reduced mass, the equation of the radial motion becomes: $$\ddot{r}-r{\dot{f}}^2+\frac{G{m}_{\text{T}}}{r^2}+\frac{9G{m}_{\text{T}}Q\left(t\right)}{2m_{\text{s}}{r}^4}=0,$$(A.8)

where m_T ≡ m_s + m_p is the total mass of the system.

Appendix A.2 Radial motion

We consider an orbit of small eccentricity, so that $$r\left(t\right)=a\left[1+x\left(t\right)\right],$$(A.9)

where x(t) is the relative deviation from a circular orbit with |x(t)| ≪ 1. Making use of the conservation of the total angular momentum J as given by Eq. (A.7), the equation of the radial motion (A.8) becomes $$\ddot{r}-\frac{J^2}{m^2{r}^3}+\frac{G{m}_{\text{T}}}{r^2}+\frac{9GQ\left(t\right)m_{\text{T}}}{2m_{\text{s}}{r}^4}=0.$$(A.10)

Substituting Eq. (A.9) into Eq. (A.10) and considering that $\ddot{r}=a\ddot{x}$ and (1 + x)^−q ≃ (1 − qx), we find $$a\ddot{x}-\frac{J^2}{m^2{a}^3}\left(1-3x\right)+\frac{G{m}_{\text{T}}}{a^2}\left(1-2x\right)+\frac{9GQ{m}_{\text{T}}}{2m_{\text{s}}{a}^4}\left(1-4x\right)=0.$$(A.11)

Making use of fact that J ≃ mna² for e ≪ 1 (cf. Eq. 1), Kepler’s third law, and collecting all the terms in x, we obtain, to first order, $$\ddot{x}+n^2\left(1-\frac{18Q}{m_{\text{s}}{a}^2}\right)x=-\frac{9}{2}{n}^2\frac{Q}{m_{\text{s}}{a}^2}.$$(A.12)

The second term in brackets on the left-hand side is very small in comparison to the first because the quadrupole moment fluctuations |Q| ≪ m_sa². In conclusion, we can write the equation of the epicyclic motion as $$\ddot{x}+n^{2}x\simeq -\frac{9}{2}{n}^2\frac{Q\left(t\right)}{m_{\text{s}}{a}^2},$$(A.13)

which is the equation of a forced harmonic oscillator of frequency n equal to that of the orbital motion. The time-dependent nature of the forcing has been made explicit by indicating that Q is a function of time t. The relationship between x and the eccentricity e comes from the solution of Eq. (A.13), that is, x = e sin(nt).

Tides inside the star and the planet damp any initial orbital eccentricity over a characteristic timescale τ_e = e/|de/dt|, which can be calculated from the tidal theory (cf. Eq. 7). Therefore, Eq. (A.13) must be completed by adding a term that takes into account the damping of the radial motion by the tides, because they tend to restore a circular orbit. This term can be expressed as a dissipative term in the equation for the radial oscillations, which becomes $$\ddot{x}+2b\dot{x}+n^{2}x\simeq -\frac{9}{2}{n}^2\frac{Q\left(t\right)}{m_{\text{s}}{a}^2}.$$(A.14)

This is the equation of a damped harmonic oscillator with a random forcing because Q(t) is a stochastic function of time. Its solution for b ≪ n is x ~ exp(−bt) sin(nt), which can be compared with the undamped solution x = e sin(nt) to give the relationship between b and τ_e, that is, $b={\tau }_{\text{e}}^{-1}$ (cf. Landau & Lifshitz 1969, § 25).

Taking the Fourier transform of both sides of Eq. (A.14) with $i=\sqrt{-1}$, we obtain $$-{\omega }^2\tilde{x}\left(\omega \right)+2bi\omega \tilde{x}\left(\omega \right)+n^2\tilde{x}\left(\omega \right)=-\frac{9}{2}{n}^2\frac{\tilde{Q}\left(\omega \right)}{m_{\text{s}}{a}^2},$$(A.15)

where the tilde indicates the Fourier transform and ω is the frequency. This can be recast as $$\tilde{x}\left(\omega \right)=-\frac{9}{2}{n}^2\frac{\tilde{Q}\left(\omega \right)}{m_{\text{s}}{a}^2}\frac{1}{Z\left(\omega \right)},$$(A.16)

where the function $$Z\left(\omega \right)=\frac{1}{n^2-{\omega }^2+2bi\omega }.$$(A.17)

The power spectrum $P_x\left(\omega \right)\equiv \tilde{x}\left(\omega \right){\tilde{x}}^{*}\left(\omega \right)$ of the solution x(t) can be immediately computed as $$\begin{array}{l}{P}_x\left(\omega \right)=\frac{81}{4}{n}^4\frac{P_Q\left(\omega \right)}{m_{\text{s}}^2{a}^4}\frac{1}{Z\left(\omega \right)Z^{*}\left(\omega \right)}=\hfill \\ =\frac{81}{4}{n}^4\frac{P_Q\left(\omega \right)}{m_{\text{s}}^2{a}^4}\frac{1}{{\left(n^2-{\omega }^2\right)}^2+4b^2{\omega }^2},\hfill \end{array}$$(A.18)

where the asterisk indicates complex conjugation and $P_Q\equiv \tilde{Q}\left(\omega \right){\tilde{Q}}^{*}\left(\omega \right)$ is the power spectrum of the quadrupole fluctuations. The expectation value of x², that is, 〈x²〉, given that 〈x〉 = 0, is given by $$\begin{array}{l}\langle x^2\rangle =\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }{P}_x\left(\omega \right)\text{d}\omega =}\hfill \\ =\frac{81}{8\pi }{n}^4\frac{1}{m_{\text{s}}^2{a}^4}{\displaystyle {\int }_{-\infty }^{\infty }\frac{P_Q\left(\omega \right)}{{\left(n^2-{\omega }^2\right)}^2+4b^2{\omega }^2}d\omega .}\hfill \end{array}$$(A.19)

Because b ≪ n, the integrand is a very sharply peaked function around ω ~ n; that is, only the power spectrum of the quadrupole fluctuations very close to the resonance contributes to the integral. Therefore, without any significant loss of accuracy, we can write $$\langle x^2\rangle \simeq \frac{81}{8\pi }{n}^4\frac{1}{m_{\text{s}}^2{a}^4}{P}_Q\left(n\right){\displaystyle {\int }_{-\infty }^{\infty }\frac{d\omega }{{\left(n^2-{\omega }^2\right)}^2+4b^2{\omega }^2}}.$$(A.20)

The integral can be computed, for example, by means of the method of the residues of the theory of complex variables and its value is (cf. Appendix of Lanza 2021) $${\displaystyle {\int }_{-\infty }^{\infty }\frac{d\omega }{{\left(n^2-{\omega }^2\right)}^2+4b^2{\omega }^2}}=\frac{\pi }{2b{n}^2}.$$(A.21)

Making use of Eq. (A.21), expression (A.20) becomes $$\langle x^2\rangle =\frac{81}{16}\frac{n^2}{b}\frac{P_Q\left(n\right)}{m_{\text{s}}^2{a}^4}.$$(A.22)

The expectation value of the squared velocity $\langle {\dot{x}}^2\rangle$ can be similarly obtained by taking into account that $\tilde{\dot{x}}=i\omega \tilde{x}\left(\omega \right)$, so that $$P_{\dot{x}}\left(\omega \right)={\omega }^2{P}_x\left(\omega \right).$$(A.23)

The equation analogous to Eq. (A.20) is then $$\begin{array}{l}\langle {\dot{x}}^2\rangle =\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }{\omega }^2}{P}_x\left(\omega \right)d\omega \simeq \hfill \\ \simeq \frac{81}{8\pi }{n}^4\frac{1}{m_{\text{s}}^2{a}^4}{P}_Q\left(n\right){\displaystyle {\int }_{-\infty }^{\infty }\frac{{\omega }^{2}d\omega }{{\left(n^2-{\omega }^2\right)}^2+4b^2{\omega }^2}}.\hfill \end{array}$$(A.24)