Issue 
A&A
Volume 674, June 2023



Article Number  A176  
Number of page(s)  17  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202345860  
Published online  20 June 2023 
Residual eccentricity of an Earthlike planet orbiting a red giant Sun
^{1}
INAFOsservatorio Astrofísico di Catania,
Via S. Sofia 78,
95123
Catania, Italy
email: antonino.lanza@inaf.it
^{2}
LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris,
5 Place Jules Janssen,
92195
Meudon, France
^{3}
Université de Rennes, CNRS, IPR (Institut de Physique de Rennes) – UMR 6251,
263 av. Général Leclerc,
35042
Rennes cedex, France
^{4}
École Universitaire de Physique et d’Ingénierie, Université ClermontAuvergne,
Site des Cézeaux, 4 Av. Blaise Pascal,
63178
Aubière, France
Received:
7
January
2023
Accepted:
15
April
2023
Context. The late phases of the orbital evolution of an Earthlike planet around a Sunlike star are revisited in order to consider the effect of density fluctuations associated with convective motions inside the star.
Aims. Such fluctuations produce a random perturbation of the stellar outer gravitational field that excites a small residual eccentricity in the orbit of the planet. This counteracts the effects of tides, which tend to circularize the orbit.
Methods. We computed the power spectrum of the outer gravitational field fluctuations of the star in the quadrupole approximation and studied their effects on the orbit of the planet using a perturbative approach. The residual eccentricity is found to be a stochastic variable showing a Gaussian distribution.
Results. Adopting a model of the stellar evolution of our Sun computed with Modules for Experiments in Stellar Astrophysics (MESA), we find that the Earth will be engulfed by the Sun when it is close to the tip of the red giant branch phase of evolution. We find a maximum mean value of the residual eccentricity of ~0.026 immediately before engulfment. Considering an Earthmass planet with an initial orbital semimajor axis sufficiently large to escape engulfment, we find that the mean value of the residual eccentricity is greater than 0.01 for an initial separation of up to ~l.4 au.
Conclusions. The engulfment of the Earth by the red giant Sun is found to be a stochastic process instead of being deterministic as assumed in previous studies. If an Earthlike planet escapes engulfment, its orbit around its remnant white dwarf (WD) star will be moderately eccentric. Such a residual eccentricity of on the order of a few hundredths can play a relevant role in sustaining the pollution of the WD atmosphere by asteroids and comets, as observed in several objects.
Key words: planetstar interactions / stars: latetype / planetary systems
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
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 energygeneration 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^{−5} cause the red giant to have a slightly fluctuating outer gravitational potential that can affect the orbit of a closeby 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^{−5}) 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 loweccentricity millisecond pulsars are the end products of lowmass Xray 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 Xray 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 longterm 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) reexamined 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 massloss 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 Earthlike planet along the evolution of a Sunlike 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 transittime 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 transittime variations is presented in Sect. 2, while an internal structure and evolution model of a Sunlike star is introduced in Sect. 3. We present an illustrative application of our approach to an EarthSunlike 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 transittime fluctuations, while its application to an Earthlike planet in orbit around a Sunlike 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 transittime 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 massloss rates and several evolution codes. Initially, most of the investigations were limited to solarlike or intermediatemass 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 Neptunemass planet and an outer Jupitermass 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 transittime fluctuations and show a simple application to the case of a system consisting of a single Earthmass planet orbiting a Sunlike 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 Sunlike 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 (1)
where m ≡ m_{p}m_{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: (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 (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 starplanet system by producing variations in its orbital elements that can no longer be regarded as constant as in the case of the standard twobody problem. Osculating orbital elements can be introduced to describe the instantaneous orbit under the effects of the massloss perturbation (e.g., Dosopoulou & Kalogera 2016a). From a dynamical point of view, the massloss regime can be characterized through the parameter (cf. Veras et al. 2011) (4)
When Ψ_{ms} ≪ 0.1, the massloss timescale is much longer than the orbital period and the system is in the socalled adiabatic massloss 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^{2}), 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 solarmass star and a planet with an orbital period of 2 yr, the parameter Ψ_{ml} ~ 1.3 × 10^{−6} 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^{−2} 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)^{6}, 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 Earthmass planet on an Earthlike 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) (5)
where λ_{2} ≃ 0.019 α^{4/3} and are a nondimensional 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 Sunlike 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 λ_{2} 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 solarmass 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 (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 (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 solarlike 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 massloss rate between the stellar hemispheres of up to 1%, with an average massloss rate of ~3 × 10^{−7} M_{⊙} yr^{−1}. 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 (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 (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 (10)
where ρ is the density, r ≡ (x_{1}, x_{2}, x_{3}) 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 (11)
where ρ_{0}(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 (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 inertia tensor are different from zero. In other words, (13)
Therefore, only the fluctuating parts 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 (15)
Therefore, the gravitational potential becomes (16)
where we have separated the isotropic component of the quadrupole potential from the component depending on the coordinates x = x_{1}, y = x_{2}, z = x_{3} 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 transittime 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 righthand 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^{1}. 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: (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: (18) (19) (20)
The only two relevant Gauss equations for the variation in the orbital elements are (21) (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) (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 (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 (25)
where we make the dependence of the tidal damping timescale τ_{e} on the time explicit (cf. Eq. (7)) and define (26)
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 (28)
where C_{0} is a constant depending on the initial value of e and (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 (30)
The term containing the constant factor C_{0} 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_{0} have been dropped from Eq. (30). The ensemble mean commutes with the integration, allowing us to evaluate 〈e^{2}(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 convective motions. In other words, we can regard f(t) as being a deltacorrelated process, that is, (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 (32) (33)
where we have explicitly introduced the slow time dependence of D that is relevant on the timescales over which 〈e^{2}〉 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^{2}〉 can be obtained by a simple integration from Eq. (33) as (34)
for t ≫ τ_{e}/2. On the other hand, 〈e^{2}〉 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 (35)
where ω is the frequency and the imaginary unit. The inverse transform is (36)
Taking the Fourier transform of both sides of Eq. (25), we have (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 (38)
where is the power spectrum of the eccentricity and the asterisk indicates complex conjugation. Making use of Eq. (37), Eq. (38) becomes (39)
where 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 (40)
To evaluate P_{f}(0), we use Eq. (26), the Euler formula for the sine, and Eq. (35) to compute (41) (42)
where is the Fourier transform of Q(t), which, being a real function, verifies . Given that Q(t) is a stochastic variable, in a statistical sense, where P_{Q}(ω) is the power spectrum of the fluctuating quadrupole moment Q(t). By applying this result, we have (43)
which can be compared with Eq. (34) to find (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 (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^{2}〉 – according to the method proposed by Phinney (1992) – is presented in Appendix A.
The residual eccentricity 〈e^{2}〉 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 MaxwellBoltzmann 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 (46)
Therefore, the probability P(e ≤ e_{0}) that the residual eccentricity e is smaller than or equal to a given value e_{0} can be expressed in terms of the error function as (47)
where we define the error function as (48)
so that lim_{z→∞} erf(z) = 1.
2.3.2 Mean longitude at the epoch
The formal solution of Eq. (22) is (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 〈ϵ^{2}〉 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 ϵ^{2} is (50)
Considering much longer time intervals than the correlation time of the quadrupole moment fluctuations, we can assume a deltacorrelated process with 〈Q(t′)Q(t″)〉 = Eδ(t″ − t′), finding 〈ϵ^{2}〉 = 4K^{2} 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 (51)
where Q_{j}(t) is the contribution to the quadrupole moment fluctuation at the time t coming from the jth 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, (52)
The probability 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 (53) (54)
where we develop the integrand function 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 righthand side of Eq. (54) is 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 in a time series and retain only the firstorder term as (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 (56)
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 (58)
Therefore, the variance of ϵ_{j} is and increases linearly in time as expected in the case of a system performing a Brownian 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 (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^{2}(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 〈ϵ^{2}〉 can be measured.
The fluctuations in the time of midtransit can be derived from Eq. (24) and are (60) (61)
where Eq. (60) follows from 〈ϵ〉 = 0, while 〈ϵ^{2}〉 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 timedependent 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^{2}(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 (62)
The quadrupole Q can be expressed as (cf. Eq. (15)) (63)
The fluctuations , and are uncorrelated with each other and have the same autocorrelation function because of the spherical symmetry of the turbulent velocity field, that is, . This implies that (64)
Therefore, we are left with the computation of the autocorrelation function of the fluctuations of the principal moment of inertia . Considering the isotropy of the fluctuations in our model (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 becomes (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 mixinglength theory, that is, a local theory of convection, the density fluctuations at a given radius r are spatially correlated over a volume of order ℓ^{3}(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 mixinglength (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 ℓ^{3} are uncorrelated with each other, we can simply integrate all the local autocorrelations of the density fluctuations and obtain^{2} (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 mixinglength 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, (68)
where ɡ(r) = Gm(r)r^{−2} is the acceleration of gravity with m(r) being the mass of the star inside the radius r. Therefore, we find (69)
Substituting Eq. (69) into Eq. (67), we find (70)
where the quantities ρ_{0}(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 , we compute the Fourier transform of the exponential decorrelation function of the time that appears in Eq. (70): (71)
The power spectrum of 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), (73)
The power spectrum of Q follows from the Fourier transform of both sides of Eq. (64), giving (74)
Specifically, P_{Q}(n), which appears in Eq. (45), is given by (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 (76)
which can be substituted into Eq. (59) after computing the integral over ω to give (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 Earthlike 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 solarmass star from the premain 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 mixinglength 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 solarmass 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 semiconvective regions, which is a delicate procedure. To cope with these difficulties, we used the new convective premixing (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 massloss parameter, inferred by Miglio et al. (2012, 2021) from asteroseismic 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 photospheric 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^{3} at the solar age, provides the solar initial helium abundance in mass fraction Y_{0} = 0.2602 and a mixinglength 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 massloss 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^{−8}, even at the tip of the RGB or in the AGB phase, when the massloss rates are the largest. Therefore, we conclude that the adiabatic massloss 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 RGBtip radius. We also point out that the recent BaSTIIAC onesolarmass model of Hidalgo et al. (2018), which includes mass loss, predicts RGBtip 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) massloss 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 HertzsprungRussell diagram, Schröder & Smith (2008) progressively and linearly decreased the value of the mixinglength 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 Sunlike 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 massloss 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 convectiveinduced residual eccentricity and transittime 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 Earthmass planet at an initial separation of 1.0 au according to Eqs. (5) and (6), where we adopt λ_{2} = 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 massloss 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^{3} 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 constantperiod 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 Earthmass 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)^{6} 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^{2}〉^{1/2} = 2.65 × 10^{−2} 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} = 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^{2}〉^{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 constantperiod 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^{2}〉^{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^{2}〉^{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 Sunlike stellar model and orbital parameters of an Earthmass planet with an initial semimajor axis of 1.0 au. Top panel: stellar radius (red solid line) and orbital semimajor 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 〈e^{2} 〉^{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 Earthmass planet with an initial semimajor 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 〈e^{2}〉^{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 Earthlike planet orbiting a Sunlike 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^{2}〉^{1/2} ~ 0.026 using our model in the case of an Earthlike 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^{10} at an angular resolution of 1 microarcsec to detect an Earthsize 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 Earthlike planet with an initial orbital semimajor axis of 1.0 au by a Sunlike 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 Earthlike planet critically depends on the stellarmassloss 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 massloss 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 massloss rate can be increased by a binary companion or a massive closeby 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 massloss 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 massloss 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 massloss 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 Earthlike 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 lowmass 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 multiplanet 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 closeby orbiting planet. Considering a time baseline of 10 yr, we find deviations in the time of midtransit with respect to a constantperiod 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 Earthlike 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 INAFOsservatorio 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 x_{0} points toward the periastron of the orbit. 
In this Appendix, we derive the value of 〈e^{2}〉 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 starplanet system
To study the dynamics of our starplanet system, we apply the Lagrangian formalism. The Lagrangian ℒ of our system is defined as (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 starplanet system and whose x_{0}y_{0} 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_{s}m_{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 (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) (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 (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: (A.5) (A.6)
Equation (A.6) can be immediately integrated to give the conservation of the orbital angular momentum J of the system: (A.7)
Considering the definition of the reduced mass, the equation of the radial motion becomes: (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 (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 (A.10)
Substituting Eq. (A.9) into Eq. (A.10) and considering that and (1 + x)^{−q} ≃ (1 − qx), we find (A.11)
Making use of fact that J ≃ mna^{2} for e ≪ 1 (cf. Eq. 1), Kepler’s third law, and collecting all the terms in x, we obtain, to first order, (A.12)
The second term in brackets on the lefthand side is very small in comparison to the first because the quadrupole moment fluctuations Q ≪ m_{s}a^{2}. In conclusion, we can write the equation of the epicyclic motion as (A.13)
which is the equation of a forced harmonic oscillator of frequency n equal to that of the orbital motion. The timedependent 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 (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, (cf. Landau & Lifshitz 1969, § 25).
Taking the Fourier transform of both sides of Eq. (A.14) with , we obtain (A.15)
where the tilde indicates the Fourier transform and ω is the frequency. This can be recast as (A.16)
The power spectrum of the solution x(t) can be immediately computed as (A.18)
where the asterisk indicates complex conjugation and is the power spectrum of the quadrupole fluctuations. The expectation value of x^{2}, that is, 〈x^{2}〉, given that 〈x〉 = 0, is given by (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 (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) (A.21)
Making use of Eq. (A.21), expression (A.20) becomes (A.22)
The expectation value of the squared velocity can be similarly obtained by taking into account that , so that (A.23)
The equation analogous to Eq. (A.20) is then (A.24)
Applying again the method of the residues, the integral is found to be (cf. Appendix of Lanza 2021) (A.25)
Therefore, the average mechanical energy of the epicyclic oscillator is (A.27)
where m ≃ m_{p} is the reduced mass of the system.
As we see above, the stationary solution of the equation for the epicyclic radial motion can be approximated as x(t) = e(t) sin(nt), where e ≪ 1 is the eccentricity of the orbit. Because b ≪ n, the timescale upon which e varies is significantly longer than the orbital period. Therefore, 〈x^{2}〉 = (1/2)〈e^{2}〉, , and the average value of the total mechanical energy is (A.28)
Comparing Eq. (A.28) with Eq. (A.27), we find the average value of the squared residual eccentricity induced by the convective fluctuations of the quadrupole moment of the star (A.29)
It is interesting to compare this result with that given by Eq. (45) for a constant τ_{e}. Developing the exponential into a Taylor series and considering only the leadingorder term, the expression derived from Eq. (45) is the same as Eq. (A.29) with τ_{e} replaced by t, that is (A.30)
Therefore, Eq. (A.29) is the leadingorder term of the exact solution as given by Eq. (45) provided that we assume t = τ_{e}. Among the implicit hypotheses in the simplified approach adopted to derive Eq. (A.29) is the constancy of τ_{e} and the consideration of a weak eccentricity damping that is valid for t ≪ τ_{e} which is the same for the validity of the truncation of the general solution (45) to the leading order.
Equation (A.30) shows that 〈e^{2}〉 increases in proportion to time, as expected in the case of a Brownian motion. However, the maximum timescale for which such a model is valid is t ≲ τ_{e} because tides erase the eccentricity on timescales comparable to or longer than τ_{e}, and so the system loses memory of the previous accumulation steps. Therefore, regarding τe as the memory timescale for the growth of 〈e^{2}〉, Eq. (A.30) reproduces the result of Eq. (A.29).
The power dissipated by the damping of the radial oscillations due to the tides can be derived by multiplying Eq. (A.14) by ma^{2} and writing the resulting equation as (A.31)
The dissipated power is given by the first term on the righthand side of this equation because the other term represents the power of the forcing that maintains the oscillations. This implies that the average dissipated power per oscillation cycle P_{d} is (A.32)
thanks to Eq. (A.26); the value of P_{Q}(n) is given by Eq. (75). In our model, tidal dissipation occurs inside the star, and therefore this power makes a very small addition to the stellar luminosity. In the stationary regime, this power is independent of the eccentricity damping timescale and scales as a^{−8}, as can be seen by substituting from the Kepler’s third law into Eq. (A.33).
Appendix B Statistical distribution of the residual eccentricity
The statistical distribution of the residual eccentricity can be computed following the same method used to derive the statistical distribution of the velocities of the molecules in an ideal gas, the socalled MaxwellBoltzmann distribution.
As the direction of the periastron in the orbital plane is randomly oriented, the distribution of the argument of periastron of the planetary orbit is uniform in [0, 2π]. In other words, the distributions of the stochastic variables (B.1)
are the same because there is no preferred orientation of the line of the apsides. We indicate this distribution as φ(e_{x}) = φ(e_{y}), and it must be symmetric, that is, φ(z) = φ(−z), where z = e_{x} or e_{y}, and will therefore depend on or .
The probability that the eccentricity vector e ≡ (e_{x}, e_{y}) has components between e_{x} and e_{x} + de_{x} and e_{y} and e_{y} + de_{y} is . This probability must depend only on the square of the eccentricity , and it can therefore be expressed as f_{e}(e^{2}) de_{x} de_{y}, where f_{e} is the distribution function of e that depends on e^{2}. In conclusion, we find (B.2)
By taking the logarithms and deriving both sides of Eq. (B.2) with respect to and , respectively, we find (B.3) (B.4)
Equation (B.3) tells us that is independent of , while Eq. (B.4) tells us that it is independent of . Therefore, must be equal to a constant, which we denote −ζ, where ζ > 0 for normalization reasons, and (B.5)
This equation can be immediately integrated to give (B.6)
The distribution function must be normalized as (B.7)
On the other hand, the mean value of the square of the eccentricity is (B.8)
Equations (B.7) and (B.8) can be used to find the normalization constant and express ζ in terms of 〈e^{2}〉. This can be done considering that 〈e^{2}〉 ~ 10^{−5} − 10^{−4} is always very small, and so the exponential becomes almost zero well before reaching the upper limit of integration e = 1 . In other words, we can extend the upper limit of integration to +∞ without any appreciable error and make use of the notable integrals (B.9)
which give the distribution function of the residual eccentricity as (B.13)
In other words, it is a Gaussian distribution with zero mean and standard deviation 〈e^{2}〉^{1/2}.
Appendix C The final fate of the Earth inside the RGB solar envelope
When the Earth enters the envelope of the red giant Sun, it will experience hydrodynamical drag and gravitational drag. The ratio of the former to the latter is on the order of (υ_{orb}/υ_{esc})^{2} ~ 10, where υ_{orb} is the orbital velocity and υ_{esc} the escape velocity at the surface of our planet (cf. Sect. 2.1 of Yarza et al. 2022). Therefore, we neglect the gravitational drag and consider only the hydrodynamical drag given by , where C_{D} is the drag coefficient of order unity, ρ the density in the stellar envelope at the orbital radius of the planet, and R_{e} the radius of the Earth.
According to Jia & Spruit (2018), the disruption of the Earth takes place at the disruption radius r_{d}, where the drag pressure becomes equal to the density of the gravitational binding energy of the Earth, that is, where the ratio (C.1)
with being the mean density of the Earth and m_{e} the Earth’s mass. Once this condition is reached, the planet starts to fragment into smaller and smaller pieces over a characteristic fragmentation timescale , where ρ_{d} is the mean stellar density inside the disruption radius and P_{oгb} the orbital period around the center of the Sun.
Considering our internal structure model of the Sun at the tip of the RGB phase, we estimate that the disruption radius is r_{d}/R ~ 0.003, where R is the radius of the star, because of the very small density of the stellar plasma in the extended convective envelope. The disruption timescale is very short (t_{d} ~ 800 s), and therefore when the condition f ~ 1 is reached, the Earth will be rapidly fragmented and dissolved in the Sun’s envelope (cf. Sect. 5.1 of Jia & Spruit 2018). We note that thermal evaporation and ablation play a secondary role in the destruction of the Earth because the density in the solar red giant envelope is much lower than the mean density of the Earth (see Sect. 2.4 of Jia & Spruit 2018, for details).
The decay of the Earth’s orbit is ruled by the variation in its orbital angular momentum J, which can be written as (C.2)
where m ~ m_{e} is the reduced mass of the Earth orbiting the Sun, a the orbit semimajor axis, and n = 2π/P_{orb} the mean orbital motion. Making use of Kepler’s third law, we derive a differential equation for the evolution of the semimajor axis under the effect of the drag, which is assumed to be tangential to the circular orbit, that is, (C.3)
In the layers where the Earth orbits before its disruption, that is, where r ≥ r_{d}, we approximate the density with a power law as ρ(r) = ρ(R)(r/R)^{q}, where ρ(R) is the density at the surface of the star where r = R and q is a fixed exponent. Such an approximation is valid within one order of magnitude, and therefore it is suitable for our approximate treatment of the problem. By fitting the internal density stratification in the model of the Sun at the RGB tip for r ≥ r_{d}, we find q = −2.13.
Making use of the above approximation for the density, and considering that r_{d} ≪ R, we integrate analytically Eq. (C.3) and find the timescale t_{p} of the Earth spiralin from the engulfment at the surface r = R down to the disruption radius r = r_{d} as (C.4)
Adopting C_{D} = 1, we find t_{p} ~ 3 × 10^{3} years because the Earth’s spiralin is initially very slow owing to the very low density in the extended solar envelope.
The dissipation rate of the orbital energy E_{oгb} due to the drag force can be estimated by assuming that the mass inside the orbit of the planet is constant and equal to the total mass of the star as (C.5)
where da/dt is given by Eq. (C.3).
The minimum dissipation rate is that at the engulfment when a = R, and it is only ~ 10^{−8} of the stellar luminosity at the tip of the RGB phase. Immediately before the disruption, a = r_{d} and the maximum of the dissipation rate is reached, which amounts to about ten times the stellar luminosity. However, such an additional energy input is limited to a short burst occurring on a timescale on the order of a few times the shortest orbital period, that is, a few days at most. The excess heat is redistributed by convection over the entire convective envelope on a timescale R^{2}/η_{turb} ~ 200 years, where η_{turb} is the turbulent diffusivity of the envelope that is estimated as the average of the product of the mixing length by the convective velocity over the envelope itself. This excess heat is finally radiated at the stellar surface on a timescale comparable to the KelvinHelmoltz timescale of the envelope itself, that is, on the order of 25 years, and therefore its impact on the stellar luminosity is negligible.
References
 Adams, F. C., & Bloch, A. M. 2013, ApJ, 777, L30 [NASA ADS] [CrossRef] [Google Scholar]
 Antoniadou, K. I., & Veras, D. 2016, MNRAS, 463, 4108 [NASA ADS] [CrossRef] [Google Scholar]
 Antoniadou, K. I., & Veras, D. 2019, A&A, 629, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., et al. 2009, ARA&A, 47, 481 [Google Scholar]
 Assef, R. J., Gaudi, B. S., & Stanek, K. Z. 2009, ApJ, 701, 1616 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., Morbidelli, A., & Holman, M. J. 2015, ApJ, 799, 120 [CrossRef] [Google Scholar]
 Blocker’s, T. 1995, A&A, 297, 727 [Google Scholar]
 Bossini, D., Miglio, A., Salaris, M., et al. 2015, MNRAS, 453, 2290 [Google Scholar]
 Carlberg, J. K., Majewski, S. R., & Arras, P. 2009, ApJ, 700, 832 [Google Scholar]
 Chamandy, L., Blackman, E. G., Nordhaus, J., et al. 2021, MNRAS, 502, L110 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1961, Hydrodynamic and Hydromagnetic Stability, International Series of Monographs on Physics (Oxford: Clarendon) [Google Scholar]
 Chaussidon, M. 2007, in Lectures in Astrobiology, eds. M. Gargaud, H. Martin, & P. Claeys (Berlin, Heidelberg: Springer), 45 [CrossRef] [Google Scholar]
 Chen, H.L., Tauris, T. M., Han, Z., et al. 2021, MNRAS, 503, 3540 [NASA ADS] [CrossRef] [Google Scholar]
 Dosopoulou, F., & Kalogera, V. 2016a, ApJ, 825, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Dosopoulou, F., & Kalogera, V. 2016b, ApJ, 825, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Dréau, G., Lebreton, Y., Mosser, B., et al. 2022, A&A, 668, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Duncan, M. J., & Lissauer, J. J. 1998, Icarus, 134, 303 [NASA ADS] [CrossRef] [Google Scholar]
 Frewen, S. F. N., & Hansen, B. M. S. 2014, MNRAS, 439, 2442 [NASA ADS] [CrossRef] [Google Scholar]
 Goldstein, H. 1950, Classical Mechanics (Reading, MA: AddisonWesley) [Google Scholar]
 Grevesse, N., & Sauval, A. J. 1998, Space Sci. Rev., 85, 161 [Google Scholar]
 Guo, J., Lin, L., Bai, C., et al. 2016, Ap&SS, 361, 122 [NASA ADS] [CrossRef] [Google Scholar]
 Henyey, L., Vardya, M. S., & Bodenheimer, P. 1965, ApJ, 142, 841 [Google Scholar]
 Hidalgo, S. L., Pietrinferni, A., Cassisi, S., et al. 2018, ApJ, 856, 125 [Google Scholar]
 Jia, S., & Spruit, H. C. 2018, ApJ, 864, 169 [NASA ADS] [CrossRef] [Google Scholar]
 Khan, S., Hall, O. J., Miglio, A., et al. 2018, ApJ, 859, 156 [NASA ADS] [CrossRef] [Google Scholar]
 Kippenhahn, R., Weigert, A., & Weiss, A. 2012, Stellar Structure and Evolution (Berlin, Heidelberg: Springer) [Google Scholar]
 Kissin, Y. & Thompson, C. 2015, ApJ, 808, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Kunitomo, M., Ikoma, M., Sato, B., et al. 2011, ApJ, 737, 66 [NASA ADS] [CrossRef] [Google Scholar]
 Lagos, F., Schreiber, M. R., Zorotovic, M., et al. 2021, MNRAS, 501, 676 [Google Scholar]
 Landau, L. D., & Lifshitz, E. M. 1969, Course of Theoretical Physics, 2nd edn. (Oxford: Pergamon Press) [Google Scholar]
 Lanza, A. F. 2021, A&A, 653, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lanza, A. F., & Rodonò, M. 2001, A&A, 376, 165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Laskar, J. 1996, Celest. Mech. Dyn. Astron., 64, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 1997, A&A, 317, L75 [NASA ADS] [Google Scholar]
 Laskar, J., & Gastineau, M. 2009, Nature, 459, 817 [NASA ADS] [CrossRef] [Google Scholar]
 Ledoux’s, P. 1947, ApJ, 105, 305 [NASA ADS] [CrossRef] [Google Scholar]
 MacLeod, M., Cantiello, M., & SoaresFurtado, M. 2018, ApJ, 853, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Madappatt, N., De Marco, O., & Villaver, E. 2016, MNRAS, 463, 1040 [NASA ADS] [CrossRef] [Google Scholar]
 Maeder, A. 1975, A&A, 40, 303 [NASA ADS] [Google Scholar]
 Merlov, A., Bear, E., & Soker, N. 2021, ApJ, 915, L34 [NASA ADS] [CrossRef] [Google Scholar]
 Miglio, A., Brogaard, K., Stello, D., et al. 2012, MNRAS, 419, 2077 [Google Scholar]
 Miglio, A., Chiappini, C., Mackereth, J. T., et al. 2021, A&A, 645, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mogavero, F., & Laskar, J. 2021, A&A, 655, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mustill, A. J., & Villaver, E. 2012, ApJ, 761, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Nordhaus, J., & Spiegel, D. S. 2013, MNRAS, 432, 500 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
 Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [Google Scholar]
 Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
 Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [Google Scholar]
 Phinney, E. S. 1992, Philos. Trans. Roy. Soc. Lond. Ser. A, 341, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Phinney, E. S., & Kulkarni, S. R. 1994, ARA&A, 32, 591 [NASA ADS] [CrossRef] [Google Scholar]
 Privitera, G., Meynet, G., Eggenberger, P., et al. 2016a, A&A, 591, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Privitera, G., Meynet, G., Eggenberger, P., et al. 2016b, A&A, 593, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Privitera, G., Meynet, G., Eggenberger, P., et al. 2016c, A&A, 593, L15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rao, S., Meynet, G., Eggenberger, P., et al. 2018, A&A, 618, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reimers, D. 1975, Mem. Soc. Roy. Sci. Liege, 8, 369 [Google Scholar]
 Ronco, M. P., Schreiber, M. R., Giuppone, C. A., et al. 2020, ApJ, 898, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Roy, A. E. 1978, Orbital Motion (Bristol: Adam Hilger Ltd) [Google Scholar]
 Rybicki, K. R., & Denis, C. 2001, Icarus, 151, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Sabach, E., & Soker, N. 2018, MNRAS, 473, 286 [NASA ADS] [CrossRef] [Google Scholar]
 Sackmann, I.J., Boothroyd, A. I., & Kraemer, K. E. 1993, ApJ, 418, 457 [Google Scholar]
 Schröder, K.P., & Cuntz, M. 2005, ApJ, 630, L73 [Google Scholar]
 Schröder, K.P., & Smith, R. C. 2008, MNRAS, 386, 155 [CrossRef] [Google Scholar]
 Soker, N. 1998, MNRAS, 299, 1242 [NASA ADS] [CrossRef] [Google Scholar]
 Soker, N., & Harpaz, A. 2000, MNRAS, 317, 861 [NASA ADS] [CrossRef] [Google Scholar]
 Sun, M., Arras, P., Weinberg, N. N., et al. 2018, MNRAS, 481, 4077 [NASA ADS] [CrossRef] [Google Scholar]
 Tauris, T. M., & Savonije, G. J. 1999, A&A, 350, 928 [NASA ADS] [Google Scholar]
 Tauris, T. M., & van den Heuvel, E. P. J. 2006, Formation and Evolution of Compact Stellar Xray Sources (Cambridge University Press), 623 [CrossRef] [Google Scholar]
 Turyshev, S. G., & Toth, V. T. 2022, MNRAS, 515, 6122 [NASA ADS] [CrossRef] [Google Scholar]
 Ventura, P., Dell’Agli, F., Lugaro, M., et al. 2020, A&A, 641, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Veras, D. 2021, Planetary Systems Around White Dwarfs, Oxford Research Encyclopedia of Planetary Science, 1 [Google Scholar]
 Veras, D., Wyatt, M. C., Mustill, A. J., et al. 2011, MNRAS, 417, 2104 [NASA ADS] [CrossRef] [Google Scholar]
 Veras, D., Hadjidemetriou, J. D., & Tout, C. A. 2013, MNRAS, 435, 2416 [NASA ADS] [CrossRef] [Google Scholar]
 Veras, D., Georgakarakos, N., Mustill, A. J., et al. 2021, MNRAS, 506, 1148 [NASA ADS] [CrossRef] [Google Scholar]
 Verbunt, F., & Phinney, E. S. 1995, A&A, 296, 709 [NASA ADS] [Google Scholar]
 Villaver, E., & Livio, M. 2009, ApJ, 705, L81 [Google Scholar]
 Villaver, E., Livio, M., Mustill, A. J., et al. 2014, ApJ, 794, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Yarza, R., Razo Lopez, N., MurguiaBerthier, A., et al. 2022, ApJ, submitted [arXiv:2203.11227] [Google Scholar]
 Zahn, J.P. 1977, A&A, 57, 383 [NASA ADS] [Google Scholar]
 Zahn, J.P. 1989, A&A, 220, 112 [NASA ADS] [Google Scholar]
 Zeebe, R. E. 2015, ApJ, 811, 9 [NASA ADS] [CrossRef] [Google Scholar]
An alternative way to derive this result is to compute the autocorrelation of the quadrupole potential considering only the realizations in which the coordinates of the point P in the reference frame of the principal axes are (x, x, 0) so that the term in square brackets on the righthand side of Eq. (16) vanishes. Such realizations are a subset of all possible realizations of our dynamical system, but such a subset is statistically equivalent to the whole ensemble of possible realizations of our system because of the independence of its statistical properties on the orientation of the principal axes.
In the case of a very small mixing length (ℓ ≪ r), one could write the autocorrelation function by introducing a Dirac delta function of the spatial coordinates as 〈ρ′(r, t + τ)ρ′(r, t)〉 = ℓ^{3}(r)δ(r′ − r) exp[−τ/τ_{c}(r)]. In this way, the integration over dV′ is immediately performed and the double integral becomes a simple integral over dr.
We took the IAU 2015 resolution B3 values for the solar mass M_{⊙} = 1.9884 × 10^{30} kg, radius R_{⊙} = 6.957 × 10^{8} m, and luminosity L_{⊙} = 3.828 × 10^{26} W. The solar age is fixed to A_{⊙} = 4.57 × 10^{9} yr (Chaussidon 2007).
All Figures
Fig. 1 Late evolution of the main parameters of our Sunlike 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. 

In the text 
Fig. 2 Stellar massloss parameter Ψ_{ml} vs. the stellar age in the evolutionary phases characterized by a significant stellar mass loss rate. 

In the text 
Fig. 3 Late evolution of the radius of our Sunlike stellar model and orbital parameters of an Earthmass planet with an initial semimajor axis of 1.0 au. Top panel: stellar radius (red solid line) and orbital semimajor 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 〈e^{2} 〉^{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. 

In the text 
Fig. 4 Same as Fig. 3, but for an Earthmass planet with an initial semimajor 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. 

In the text 
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 〈e^{2}〉^{1/2}; bottom panel: tidal decay timescale of the eccentricity τ_{e}. 

In the text 
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 x_{0} points toward the periastron of the orbit. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.