© The Authors 2022

1 Introduction

The late 20th century saw the birth of the field of exoplanetology with the detection of the first exoplanet around a sun-like star, 51 Pegasi-b (Mayor & Queloz 1995). The improvement of planetary detection and characterisation techniques in the 21st century has allowed a confirmed detection of 5θ30 planets (Akeson et al. (2013) discussed NASA’s Exoplanet Archive, which provides data on these systems). Among all the planets discovered to date, hot, massive, close-in planets called hot Jupiters are in the most extreme environment and sometimes found with high eccentricities (Grunblatt et al. (2018, 2022), in paricular when their periods are longer than five days. These discoveries led to an improvement of the theoretical understanding and description of such systems, which can present different architecture to that of our Solar System. The study of the dynamical evolution and stability of planetary systems is fundamental in order to retrieve the formation of these systems in their variety and diversity, and the role played by tidal forces in this context is crucial.

Three channels of formation are favoured for hot Jupiters: in-situ formation, migration caused by the protoplanetary disc interaction, or high-eccentricity tidal migration (Dawson & Johnson 2018). The last scenario is particularly interesting from a dynamical point of view as the tidal dissipation plays a key role in determining the evolution of the system. Hot Jupiters are characterised by significantly inflated radii when compared with the one of Jupiter (Bodenheimer et al. 2001). Guillot & Showman (2002) suggested that irradiation can slow the Kelvin-Helmholtz contraction of hot Jupiters. Moreover, evidence that re-inflation can occur rapidly in post-main-sequence stars are presented in Grunblatt et al. (2016); Thorngren et al. (2021) and Grunblatt et al. (2022), while Thorngren et al. (2021) also came to the conclusion that a re-inflation could occur during the host star’s main sequence (MS). The drastic change of the planetary structure can have an important impact on its dynamics; in particular, a change in the structure may lead to a change in the efficiency of tidal dissipation within the planet. Constraining the tidal dissipation efficiency is challenging, and only weak constraints have been derived from dynamical studies of extra-solar planets. In the Solar System, Irwin (2012) constrained the tidal dissipation of Jupiter with an accuracy of about 20%, while Saturn’s tidal dissipation was constrained with an accuracy of 30% (Lainey et al. 2012). A constraint on the tidal dissipation of Hot-Jupiters would be a significant improvement in the modelling of a multi-planetary system as such planets can strongly affect the dynamical evolution of smaller planets.

In this study, we investigated the dynamical evolution of close-in planets orbiting red giant branch (RGB) stars, aiming to constrain the tidal dissipation of the single hot Jupiter orbiting the RGB star Kepler-91 (KIC 8219268). The particular interest in the Kepler-91 system stems from the presence of a close-in orbiting hot-Jupiter-like planet, namely Kepler-91b, which orbits its host star at one of the smallest semi-major axes over stellar radius ratio ever seen in any catalogue. This planet is expected to undergo strong tidal interaction and to be engulfed by its host star in about 55 million years (Lillo-Box et al. 2014a,b).

Asteroseismology has recently become the most accurate method to precisely characterise and model solar oscillating stars. The advent of long space-based photometric surveys such as CoRoT (Auvergne et al. 2009), Kepler (Borucki et al. 2010), and TESS (Ricker et al. 2014, 2015) gave the material to conduct precise asteroseismic studies on various types of stars and in particular exoplanet host stars. As knowledge of the dynamics of a planetary system is proportional to knowledge concerning its host star, asteroseismology appears to be an ideal way to characterise stars allowing the precise dynamical study of the planetary systems orbiting it. Among all observed RGB stars, which are evolved solar stars, are particularly interesting physics laboratories. The extended envelope and contracted core of RGB stars allows mixed modes to propagate inside all the layers, meaning that the very deep rotation profile and structure of the star can be probed. The extended convective zone of RGB stars enhances the interaction between the planet and the star through the tidal force (Villaver & Livio 2009; Villaver et al. 2014; Rao et al. 2018). Moreover, during the stellar RGB phase, the whole planetary evolution is drastically speeded up due to the fast and dramatical changes of the stellar properties (Gamow 1939).

The strategy used in this work to study Kepler 91b is based on a detailed modelling of the entire system. First, the star is modelled independently with the Code Liégeois d’Évolution Stellaire (CLES, Scuflaire et al. 2008b) and the Liège OScillation Code (LOSC, Scuflaire et al. 2008a) taking particular care in reproducing the observed stellar classical constrains and asteroseismic quantities. Then the evolution of the star is coupled to the one of the planets with our orbital evolution code (Privitera et al. 2016a,b; Meynet et al. 2017; Rao et al. 2018) to study the newly implemented evolution of the planet orbital eccentricity. In particular, the stellar equilibrium tides are expected to play a key role in the evolution of the system during the stellar RGB phase (Villaver & Livio 2009; Villaver et al. 2014; Rao et al. 2018). The physics and the dissipation of the planetary tides are closely linked to dissipation factors: k_2,pΔt_p (Darwin 1879) or Q′_p (Goldreich 1963; Goldreich & Soter 1966). In our study, we tried to determine boundaries to these quantities by estimating the maximum eccentricity reachable by the planet for different values of k_2,p∆t_p or Q′_p.

As a good description of the structure of the star and its evolution is required to precisely study the dynamical evolution of its planet, we started by carrying out a detailed stellar modelling of Kepler 91, which is presented Sect. 2. Section 3 is dedicated to a brief presentation of the planetary system and properties. The orbital evolution code used to model the evolution of the planet and the equations describing the evolution of the orbital eccentricity are presented in Sect. 4. The detailed study of the past dynamical evolution of Kepler 91b and the consequence of our study on the tidal dissipation factors are presented in Sect. 5. Finally, a review of our results and the potential improvements to our method are presented in Sect. 6.

2 Stellar models and properties

To study the evolution of the planetary system around Kepler-91, a reliable model of the evolution of the star should be obtained. To that end, we combined global and local minimisation techniques using seismic and non-seismic constraints to derive a robust model of our target. The coupling between the star and the planetary system is comprehensively taken into account by the orbital evolution code (Privitera et al. 2016a,b; Rao et al. 2018); thus, the star was considered as isolated and modelled independently.

2.1 Observational constraints

The star was modelled based on asteroseismic and classical constraints as in Lillo-Box et al. (2014b). In Table 1, we summarise the global seismic indexes as well as the classical constraints for Kepler-91, while a table with the full seismic data adopted is given in A.1.

In Table 1, ν_max denotes the frequency of the oscillation spectrum where the higher amplitude is reached. The luminosity presented in Table 1 was computed from the following formula (Casagrande & VandenBerg 2014): $$\log \left(\frac{L}{L_{\odot }}\right)=-0.4\left(m_{\lambda }+B{C}_{\lambda }-5\log d+5-A_{\lambda }-M_{\text{bol,}\odot }\right),$$(1)

where m_λ, BC_λ, and A_λ are the apparent magnitude, bolometric correction, and extinction in a given band λ,respectively. We used the 2MASS K-band magnitude properties. The bolometric correction is estimated using the code written by Casagrande & VandenBerg (2014, 2018), and the extinction is inferred with the Green et al. (2018) dust map. A value of M_bol,⊙ = 4.75 is adopted for the solar bolometric magnitude (Mamajek et al. 2015). Using Gaia DR2 (Gaia Collaboration 2018), we obtain a luminosity value of L = 16.745 ± 2.578 L_⊙.

2.2 Stellar modelling

The stellar modelling procedure follows that of Fellay et al. (2021). Our modelling is divided in two distinct phases; a first modelling step is carried out with the AIMS software (Reese 2016; Rendle et al. 2019) using classical and asteroseismic constraints. The second modelling step is carried out with a Levenberg-Marquardt minimisation technique coupled with the Code Liégeois d’Evolution Stellaire (CLES, Scuflaire et al. 2008b) and the Liège OScillation Code (LOSC, Scuflaire et al. 2008a), aimed at reproducing the whole oscillation spectrum of Kepler-91 by computing models on the fly rather than using a pre-defined grid.

The first modelling step is required to determine a first estimate of the stellar properties. It is then supplemented by a deeper local exploration of the parameter space to reproduce the full oscillation spectrum. We used the AGSS09 (Asplund et al. 2009) abundances for both AIMS and the Levenberg–Marquardt method. For both modelling steps, the models use the FreeEOS equation of state (Irwin 2012), the OPAL opacities (Iglesias & Rogers 1996), the T(τ) relation from Model-C of Vernazza et al. (1981) for the atmosphere, the mixing length theory of convection implemented as in Cox & Giuli (1968), and the nuclear reaction rates of Adelberger et al. (2011). The constraints and free parameters for the AIMS modelling are summarised in the second column of Table 2, while the fourth column presents the constraints and free parameters of the Levenberg-Marquardt modelling. In Table 2, L denotes the luminosity of the star, [Fe/H] its observed metallicity, X₀ the initial hydrogen mass fraction, and Z₀ the initial metal mass fraction. Convection is controlled by two parameters, the classical mixing-length parameter α_MLT and α_over characterising the length of the core overshooting region with the same formalism as the MLT The temperature gradient in the overshooting region is assumed adiabatic in both modelling steps. The mixing-length parameter was kept at a solar-calibrated value for all tracks of the AIMS grid. While fitting the individual radial modes, we took into account the impact of surface effects by using the two-term surface correction of Ball & Gizon (2014) in AIMS.

During the last modelling step, we ensure that our models reproduce the seven asteroseismic quantities d₀₂ defined as $$d_{02,n}=v_{n,0}-v_{n-1,2},$$(2)

where v_n,0 is the frequency of the oscillation mode with a radial order n and a degree ℓ = 0. The d₀₂ are less affected by the so-called surface effects. The reproduction of the d₀₂ coupled with the first radial mode ensures the reproduction of the behaviour of the ℓ = 0 and ℓ = 2 modes; however, due the mixed character of the ℓ = 1 modes we took particular care to verify that the global seismic pattern was reproduced with an Echelle diagram.

Fig. 1 Echelle diagram of Kepler 91 comparing observed oscillation frequencies (listed in Table A.1) represented as squared points and theoretical oscillation frequencies represented as circle points. The whole modelled oscillation spectrum for the ℓ = 1 modes is included with shaded points. The difference between observed and modelled data is quantified, in the legend, as a χ2 for each spherical order ℓ. Finally, the error bars on the observed frequencies are included in the data points and not clearly distinguishable here.

2.3 Stellar models

The final solutions of our structural modelling are given in Table 3. The agreement obtained in terms of pulsation frequencies is illustrated by the Echelle diagram shown in Fig. 1 for the whole pulsation spectrum, while the Hertzsprung-Russell diagram of the star modelled with the Levenberg-Marquardt minimisation technique presented in Fig. 2.

The non-seismic stellar properties obtained with a Levenberg-Marquardt minimisation technique lay within the 1 − σ errors of all the models and spectroscopic data obtained in Lillo-Box et al. (2014b). By fitting only the lowest frequency modes, we ensured a limited impact of the surface effect correction on the final results, while avoiding additional parameters for the modelling. The best model found lays in a valley of minima that was also found by Lillo-Box et al. (2014b), with a direct fitting of the individual frequencies proving the robustness and consistency of the two asteroseismic modelling approaches. Overall, this model shows a good agreement with both the non-seismic stellar properties and the asteroseismic properties.

We note that the Echelle diagram shown in Fig. 1 includes the empirical surface corrections of Ball & Gizon (2014) calibrated from acoustic oscillations of main-sequence stars and added posteriori to the stellar modelling on the ℓ = 0 and the p-dominated ℓ = 1 modes. The surface corrections were only used for display purposes once our final model was obtained, and they were not included in our detailed modelling to avoid any mis-modelling attributed to surface effects.

Fig. 2 Hertzsprung-Russell diagram of the evolution of Kepler 91 from the pre-Main Sequence to the modelled actual state of the star and the its future evolution with dotted line. The grey area corresponds to the constraints on the effective temperature and luminosity listed in Table 1.

3 Kepler-91b

In this section, we briefly present the properties of our targeted planet, namely Kepler-91b, which is the only planet known to date orbiting Kepler-91 Lillo-Box et al. (2014b,a). Three authors successfully characterised the planet Kepler 91b, namely Lillo-Box et al. (2014b,a); Sato et al. (2015) and Barclay et al. (2015).We now detail and explain the choice of the planetary parameters of Kepler 91b adopted in this study. While Lillo-Box et al. (2014b,a) and Barclay et al. (2015) found very similar system parameters, Sato et al. (2015) found a different system with a lighter planet and a much smaller a/R_*. Each author used different sets of radial velocities obtained from their respective instruments and the Kepler light curves. Assuming that the orbital period found by the different authors of 6.24668005 ± 0.00002647 days is valid, one can derive, using Kepler’s third law, the radius of the host star needed to fulfil both constraints on the semi major axis: the period and the measured a/R_*. With the stellar mass derived from our modelling, we found that the stellar radius necessary to fulfill both constraints on the semi-major axis is R_* = 6.92 ± 0.14 R_⊙ for the parameters of Sato et al. (2015) and $R_{*}={6.37}_{-0.37}^{+0.88}{R}_{\odot }$ for the parameters derived by Lillo-Box et al. (2014b). Since our stellar model has a radius of R_* = 6.299 R_⊙, we choose to not use the planetary parameters of Sato et al. (2015) due to the impossibility to fit both the period and a/R_*. Our choice was also guided by the stellar modelling and characterisation, which is absent from the study of Sato et al. (2015).

We thus chose to adopt the planetary and orbital parameters obtained by Lillo-Box et al. (2014b) in our analysis since the uncertainties encompass the parameters obtained in Barclay et al. (2015) and the analysis in Lillo-Box et al. (2014a) is only based on radial velocity measurements. The adopted characteristics of the system for the rest of the study are presented in Table 4.

Despite its unusually small inclination and high-impact parameter (Sato et al. 2015; Barclay et al. 2015), which would normally have resulted in a non-transiting planet, Kepler-91b is still observable thanks to the inflated radius of its host star. Regarding stellar rotation, Lillo-Box et al. (2014b) derived a surface rotation using spectroscopy and a global spectrum fit based on a grid of parameters, the surface rotation obtained is v sin i = 6.8 ± 0.2 km s⁻¹ or ${\text{Ω}}_{*}={265}_{-16}^{+18}$ nHz (corresponding to a period of ${44}_{-3}^{+2}$ days) using the inclination, stellar radius, and uncertainties given in Lillo-Box et al. (2014b). Excluding this rotation from being the result of a previous engulfment, we could reproduce the observed stellar surface rotation using our orbital evolution model with an initial surface rotation at the beginning of the MS of Ω_* = 0.4 Ω_crit, where Ω_crit is the critical rotation rate after which the stellar envelope is not bounded to the star. Even if this rotation can be reproduced, we noticed that the errors on the V sin i are significantly small compared to similar stars such as Kepler-56 (Huber et al. 2013), and that Kepler-91 is indeed a fast rotator. The surface rotation could also be confirmed with aster-oseismology using rotational splittings; however, none of them were clearly identified by Lillo-Box et al. (2014b).

4 The orbital evolution code

In this section, we briefly describe the structure and the physics included in the orbital evolution code used in this work and present the modifications implemented. This orbital evolution code was first introduced in Privitera et al. (2016a,b); Meynet et al. (2017); Rao et al. (2018) and recently used in Pezzotti et al. (2021b). The principle of this code is to couple the evolution of the planet with the one of the host star by giving the stellar model as input. In particular, stellar models provide quantities such as the stellar momentum of inertia necessary to follow the evolution of the stellar rotation and the dynamical interactions between the star and the planets.

In our code, we account for the evolution of the planetary orbit due to the dissipation of tides. The tidal force is resultant of the gravitational force appearing in bodies not described as point-mass objects. The tidal force drives a mass redistribution into bodies by perturbing their gravitational potentials. Tides can be divided into two components: the equilibrium tide, which is a large-scale circulation resulting from hydrostatic adjustment; and the dynamical tides (Zahn 1977), which come from the excitation of eigenmodes of oscillations. A part of the energy used to displace the fluid elements inside the body’s gravitational potential is directly extracted from the orbital energy of the system and the tidally perturbed body’s rotation. This energy is then dissipated inside bodies through friction or heat diffusion in radiative zones. The tidal dissipation affects the orbital evolution of a system through the evolution of a few orbital parameters, the semi-major axis, the eccentricity, the obliquity, the precession, and the rotation periods of the involved celestial bodies.

In the orbital evolution code, we can study the orbital evolution of single planet around its host star. We assume that the planet is on a coplanar orbit with respect to stellar equatorial plane, and we do not account for the planetary obliquity. In this new version of the code, we extend the study of the evolution of the orbit to planets on eccentric orbits. The code computes the evolution of the planetary semi-major axis caused by the tidal force; gravitational, frictional and drag forces; and migration linked to planetary mass loss (Zahn 1966, 1977, 1978, 1989; Villaver & Livio 2009; Villaver et al. 2014). Initially, our orbital evolution code took into account the equilibrium and dynamical tides dissipated in the stellar convective envelopes, while the dissipation in the planet was not accounted for. Finally, the planetary mass loss is computed using the energy limited formalism (Watson et al. 1981; Lammer et al. 2003; Pezzotti et al. 2021b). We explored the effect’s mass loss and found it to be negligible for this target.

4.1 Equilibrium tides

A precise description of the dynamical evolution of Kepler 91-b is key to constraining the tidal dissipation coefficients of the planetary equilibrium tides. To that end, and moreover to study the case of eccentric orbits, we decided to implement the evolution of the eccentricity inside our orbital evolution code (Privitera et al. 2016a,b; Meynet et al. 2017; Rao et al. 2018). Since this planet is on an eccentric orbit, as a good fraction of the hot Jupiters observed around RGB stars are, we decided to implement the evolution of the orbital eccentricity in our orbital evolution code.

From the formalism of the stellar equilibrium tides included in our orbital evolution code (see Villaver & Livio 2009; Villaver et al. 2014), derived from that of Zahn (1966, 1977, 1978, 1989), we generalised the expression of the semi-major-axis variation to include the contribution of the eccentricity. In the same formalism, the variation of eccentricity was obtained directly from the expressions of Zahn (1966, 1977, 1978, 1989). The expressions of the semi-major axis and eccentricity evolution only accounting for the equilibrium tides raised in the star are as follows: $$\begin{array}{l}{\left(\frac{1}{\text{a}}\frac{\text{da}}{\text{dt}}\right)}_{\text{eq,tides,*}}=\hfill \\ \frac{f_{\text{orb}}}{\tau }\frac{M_{\text{env}}}{M_{*}}q\left(1+q\right){\left(\frac{R_{*}}{\text{a}}\right)}^8\left[N_{\text{a2}}\left(\text{e}\right)\frac{{\text{Ω}}_{*}}{\text{n}}-N_{\text{a1}}\left(\text{e}\right)\right],\hfill \end{array}$$(3) $$\begin{array}{l}{\left(\frac{1}{\text{e}}\frac{\text{de}}{\text{dt}}\right)}_{\text{eq,tides,*}}=\hfill \\ \frac{11}{4}\frac{f_{\text{orb}}}{\tau }\frac{M_{\text{env}}}{M_{*}}q\left(1+q\right){\left(\frac{R_{*}}{\text{a}}\right)}^8\left[{\text{Ω}}_{\text{e}}\left(\text{e}\right)\frac{{\text{Ω}}_{*}}{\text{n}}-\frac{18}{11}{\text{N}}_{\text{a1}}\left(\text{e}\right)\right],\hfill \end{array}$$(4)

with $$N_{\text{a2}}\left(e\right)=\frac{1+\frac{15}{2}{e}^2+\frac{45}{8}{e}^4+\frac{5}{16}{e}^6}{{\left(1-e^2\right)}^6},$$(5) $$N_{\text{a1}}\left(e\right)=\frac{1+\frac{31}{2}{e}^2+\frac{255}{8}{e}^4+\frac{185}{16}{e}^6+\frac{25}{64}{e}^8}{{\left(1-e^2\right)}^{15/2}},$$(6) $${\text{Ω}}_{\text{e}}\left(e\right)=\frac{1+\frac{3}{2}{e}^2+\frac{1}{8}{e}^4}{{\left(1-e^2\right)}^5},$$(7) $$N_{\text{e}}\left(e\right)=\frac{1+\frac{15}{2}{e}^2+\frac{15}{8}{e}^4+\frac{5}{64}{e}^6}{{\left(1-e^2\right)}^{13/2}},$$(8)

where a denotes the semi major axis of the planet, e its eccentricity, M_env its mass, M_* the mass of the star, M_env the mass of its envelope, R_* its total radius, q is the planet-to-star mass ratio $\left(q=\frac{M_{\text{pl}}}{M_{\text{*}}}\right)$, Ω_* the stellar surface rotation rate, and n the planet orbital frequency.τ corresponds to the eddy turnover timescale (Rasio et al. 1996; Villaver & Livio 2009) defined as $$\tau ={\left[\frac{M_{\text{env}}{\left(R_{*}-R_{\text{env}}\right)}^2}{3L}\right]}^{\frac{1}{3}},$$(9)

with R_env being the radius at the base of the convective envelope, and L the bolometric luminosity.

Finally, f_orb is a numerical factor equal to 1 when $\tau <\frac{P_{\text{orb}}}{2}$ (with P_orb being the orbital period of the planet); otherwise, f_orb is expressed as $$f_{\text{orb}}={\left(\frac{P_{\text{orb}}}{2\tau }\right)}^2.$$(10)

The equilibrium tides are only considered when the star has a convective envelope. The term $N_{\text{a}2}\left(e\right)\frac{{\text{Ω}}_{*}}{n}-N_{\text{a}1}\left(e\right)$ in Eq. (3) defines whether the planet is inside the co-rotation radius, the radius at which the orbital period of the planet is equal to the stellar rotation period. The impact of tides is a widening of the orbit when the planet is beyond the co-rotation radius and shrinkage when the planet is inside it.

The equilibrium tides raised inside the planets are also included in the code with the constant time lag formalism (Darwin 1879). Assuming a planet in pseudo-synchronisation (Hut 1981) with a rotation rate of ω_pl = ω_pl,ps, the expressions for the evolution of the semi major axis and eccentricity caused by planetary equilibrium tides are $$\begin{array}{l}{\left(\frac{1}{\text{a}}\frac{\text{da}}{\text{dt}}\right)}_{\text{eq,tides,pl}}=\hfill \\ 6k_{2,\text{p}}\text{Δ}{t}_{\text{p}}\frac{1}{q}{\left(\frac{R_{\text{pl}}}{\text{a}}\right)}^5{n}^2\left[N_{\text{a2}}\left(\text{e}\right)\frac{{\omega }_{\text{pl,ps}}}{\text{n}}-{\text{N}}_{\text{a1}}\left(\text{e}\right)\right],\hfill \end{array}$$(11) $$\begin{array}{l}{\left(\frac{1}{\text{e}}\frac{\text{de}}{\text{dt}}\right)}_{\text{eq,tides,pl}}=\hfill \\ \frac{33}{2}{k}_{2,\text{p}}\text{Δ}{t}_{\text{p}}\frac{1}{q}{\left(\frac{R_{\text{pl}}}{\text{a}}\right)}^5{n}^2\left[{\text{Ω}}_{\text{e}}\left(\text{e}\right)\frac{{\omega }_{\text{pl,ps}}}{\text{n}}-\frac{18}{11}{\text{N}}_{\text{a1}}\left(\text{e}\right)\right],\hfill \end{array}$$(12)

where R_pl denotes the radius of the planet, k_2,p is the potential Love number of degree 2 of the planet, and Δt_p its constant time lag. For a rotationally synchronised planet on a circular orbit, the tidal dissipation inside the planet vanishes even if a whole spectrum of frequencies, not considered in our code, could dissipate energy.

Following Hut (1981), in the case of a planet in a state of tidal locking, the ratio $\frac{{\omega }_{\text{pl,ps}}}{n}$ can be expressed as $$\frac{{\omega }_{\text{pl,ps}}}{n}=\frac{1+\frac{15}{2}{e}^2+\frac{45}{8}{e}^4+\frac{5}{16}{e}^6}{\left(1+3e^2+\frac{3}{8}{e}^4\right){\left(1-e^2\right)}^{3/2}}.$$(13)

The tides raised in the planet can dissipate energy in the considered body through friction. In the case of a pseudo-synchronised planet, the energy dissipated inside the planet Ė_pl (Hut 1981) is given by $${\dot{E}}_{\text{pl}}=6k_{2,\text{p}}\text{Δ}{t}_{\text{p}}\frac{G{M}_{*}^2}{R_{\text{pl}}}{\left(\frac{R_{\text{pl}}}{a}\right)}^6{n}^2\left(N_{\text{a1}}\left(e\right)-\frac{N_{\text{a2}}{\left(e\right)}^2}{\text{Ω}\left(e\right)}\right),$$(14)

where Ω(e) is defined as $$\text{Ω}\left(e\right)=\frac{1+3e^2+\frac{3}{8}{e}^4}{{\left(1-e^2\right)}^{9/2}}.$$(15)

This dissipated energy is then added to the energy balance of the planet, assuming that all the tidal flux is instantly and homogeneously deposited in the planet (see Eq. (22)).

Although in our orbital evolution code the dissipation of dynamical tides in stellar convective zones is taken into account when the planet is on a circular and coplanar orbit with respect to the stellar equatorial plane, we do not take them into account here. The formalism adopted in our code (Ogilvie 2013) does not apply to the case of dynamical tides raised in the star by planet on eccentric orbits. The dynamical tides were not expected to play a key role for this particular system during the stellar MS as the dissipation of stellar tides is negligible compared to the planetary tides in the regimes considered here.

The total change of semi-major axis caused by tides is thus given by $${\left(\frac{1}{\text{a}}\frac{\text{da}}{\text{dt}}\right)}_{\text{tides}}={\left(\frac{1}{\text{a}}\frac{\text{da}}{\text{dt}}\right)}_{\text{eq.tides,*}}+{\left(\frac{1}{\text{a}}\frac{\text{da}}{\text{dt}}\right)}_{\text{eq.tides,pl}},$$(16)

and the same expression can be adopted for the change of eccentricity: $${\left(\frac{1}{\text{e}}\frac{\text{de}}{\text{dt}}\right)}_{\text{tides}}={\left(\frac{1}{\text{e}}\frac{\text{de}}{\text{dt}}\right)}_{\text{eq.tides,*}}+{\left(\frac{1}{\text{a}}\frac{\text{de}}{\text{dt}}\right)}_{\text{eq.tides,pl}}.$$(17)

4.2 Computation of the stellar rotation

As specified in the introduction, the orbital evolution code has the responsibility of coupling the stellar evolution with the planet evolution, and in particular computing the stellar rotation. In this study, we assumed that the star rotates as solid body along its evolution; this assumption is supported by helioseismic measurements as well as by asteroseismic observations of the internal rotation of solar-type stars.

The initial value of the stellar surface rotation rate is given as input in the orbital evolution code and is chosen in order to reproduce the surface rotation rate observed for Kepler-91. The braking of the stellar surface by magnetized winds is neglected for this target, since it is expected to be inefficient for a star with a mass similar to the one of Kepler-91, which exhibits a very shallow convective envelope Matt et al. (2012). Using the conservation of the angular momentum of the tidal force, the evolution of the specific angular momentum of the star resulting from the migration and circularisation of the planet can be expressed as $$\begin{array}{l}\frac{{\text{dJ}}_{*}}{\text{dt}}=\hfill \\ -{\left(G\mu a\left(1-e^2\right)\right)}^{\frac{1}{2}}\left[\frac{1}{2}{\left(\frac{1}{a}\frac{\text{da}}{\text{dt}}\right)}_{\text{tides}}-\frac{e^2}{1-e^2}{\left(\frac{1}{e}\frac{\text{de}}{\text{dt}}\right)}_{\text{tides}}\right],\hfill \end{array}$$(18)

where µ denotes the reduced mass of the system, expressed as $$\frac{1}{\mu }=\frac{1}{M_{*}}+\frac{1}{M_{\text{pl}}}.$$(19)

Fig. 3 Diagram underlying the geometric increase of irradiation of the host star onto the planet surface caused by the increase of the stellar radius.

4.3 Evolution of the planetary radius

The particular case of close-in planets orbiting RGB stars requires several adjustments on the modelling of the planet and the treatment of the planetary mass loss. The planet equilibrium temperature should be adjusted to account for the increase of the irradiation received from the host star over the planetary surface, caused by the increase of the stellar radius. This situation is illustrated in Fig. 3. Physically, because of the increase of the stellar radius during the RGB phase, a larger portion of the planetary surface is irradiated by the star, and consequently its equilibrium temperature increases. For the case of Kepler-91, Lillo-Box et al. (2014b) found that the total surface of Kepler 91b irradiated by the star is above 70%, showing that this mechanism is non-negligible for close-in planets around RGB stars. In Fig. 3, the key parameter α denotes the angle between a situation with a normal planetary irradiation and the situation taking into account the geometrical increase of irradiation of the host star onto the planet surface. Using geometry, the expression of α is given by $$\begin{array}{lll}\alpha \hfill & =\hfill & {\sin }^{-1}\left(\frac{R_{*}}{\sqrt{a^2+R_{*}^2}}\right)-{\sin }^{-1}\left(\frac{R_{\text{pl}}}{\sqrt{a^2+R_{*}^2}}\right)\hfill \\ \hfill & \approx \hfill & {\sin }^{-1}\left(\frac{R_{*}}{\sqrt{a^2+R_{*}^2}}\right).\hfill \end{array}$$(20)

This angle is non-negligible for close-in planets around giant stars, or at least when the stellar radius has the same order of magnitude as the orbital separation. For planets with wide orbits, α ~ 0 and the planet recovers its 50% of irradiated surface. To compute the irradiation of the planet, we introduced a quantity β defined as the ratio of the actually surface irradiated to the surface irradiated without these geometrical considerations. This coefficient can be expressed as $$\beta =\frac{1}{2\pi R_{\text{pl}}^2}{\displaystyle {\int }_0^{2\pi }\text{d}\varphi }{\displaystyle {\int }_0^{\frac{\pi }{2}+\alpha }{R}_{\text{pl}}^2\sin \left(\theta \right)\text{d}\theta =1+\sin \left(\alpha \right).}$$(21)

Finally, the equilibrium temperature T_eq. of the planet including the geometrical correction and including the energy dissipated through planetary tides, presented in Eq. (14), can be expressed as $$T_{\text{eq}}={\left(\frac{L_{*}\left(1-A\right)\beta }{16\pi a^2\sqrt{1-e^2}{\sigma }_{\text{SB}}}+\frac{{\dot{E}}_{\text{pl}}}{4\pi R_{\text{pl}}^2{\sigma }_{\text{SB}}}\right)}^{\frac{1}{4}},$$(22)

where L_* denotes the stellar luminosity, σ_SB the Stefan–Boltzmann constant, and A the planetary bond albedo. In the previous formula, β is treated as a simple coefficient equal to 1 for planets with wide orbits; thus, this expression is a generalisation for the classical equilibrium temperature computation. For each of the evolutionary tracks presented below, the bond albedo was kept at 0.30, the maximum observed albedo comes from the sample of hot Jupiters of Adams et al. (2022).

Determining the mechanisms that concur in the inflation of the planetary radius for hot Jupiters represents a challenging task. The high quantity of energy deposited in the atmosphere leads to a significant increase in the planetary radius impacting the dynamical evolution and the planetary mass loss. In our cases, the classic mass-radius relations can be inaccurate and underestimate the planetary radius. Moreover, the radii of hot and ultra-hot Jupiters are dominated by their irradiations more than their masses; thus, we chose to compute the planetary radius with the temperature-radius law presented in Owen & Lai (2018) only valid for planets with masses greater than 0.2 M_Jup. The expression of the initial radius of the planet given in Owen & Lai (2018) is $$R_{\text{pl}}=0.8f_{\text{r}}{\left(\frac{T_{\text{eq}}}{1100K}\right)}^{0.7224}{R}_{\text{Jup}},$$(23)

with f_r being a free parameter with a value laying between 1 and 1.5; the parameter is adjusted in our code to reproduce the observed planetary radius. Following Owen & Lai (2018), the planetary radius evolves according to $$R_{\text{pl}}=0.57792f_{\text{r}}\frac{{\dot{T}}_{\text{eq.}}}{{\left(1100K\right)}^{0.7224}{T}_{\text{eq}.}^{0.2776}}{R}_{\text{Jup}},$$(24)

where ${\dot{T}}_{\text{eq}}$ denotes the time derivative of the planetary equilibrium temperature.

5 Orbital evolution of Kepler-91b

In this section, different observational constraints of the Kepler-91 system are used to constrain the physics of tides and their dissipation. The study of the Kepler-91 system is divided into three distinct parts. First, in Sect. 5.1, we start by carrying out a theoretical comparison between the dissipation of the planet and of the star. Then, we constrain the physics of tides and their dissipation using different observational constraints on the system. Section 5.2 is dedicated to the study of the past evolution of Kepler-91b to quantify and constrain the dissipations of tides in the planet during this period of time. In Sect. 5.3, we use the constraints obtained on the planetary equilibrium tides to quantify the dissipation of the stellar equilibrium tides through the different phases of stellar evolution, both for the past and future evolution of the system.

5.1 Comparison of the dissipation in the planet and that in the star

To illustrate the variation of the different components of the tidal dissipation, we carried out an initial comparison of the migration timescale τ_a and circularisation timescale τ_e. Figure 4 illustrates the evolution of these timescales over the lifetime of the whole system, only accounting for equilibrium tides.

Figure 4 shows that during the MS of the star (ending at an age of ~ 4.5 Gyr), the circularisation is totally dominated by the planetary tides in this case. The slow migration of the planet is led both by the stellar and planetary tides during the MS, while during the stellar sub-giant phase and RGB phase the stellar tides dominate. The particular peaks seen in τ_a,star, τ_tot,star and τ_e,star during the MS are caused by the co-rotation radius almost being reached. During the stellar RGB phase, the dissipation and the engulfment of the planet are driven by the stellar tides. We also found that the circularisation acts on shorter timescales than the migration during the stellar MS, while during the engulfment both timescales have the same order of magnitude even if the circularisation still happens faster. Concerning the stellar tides, all the results presented were expected (see Villaver et al. 2014; Van Eylen et al. 2016; Grunblatt et al. 2018), confirming the accuracy of our orbital evolution code.

5.2 Constraints on the maximum value of the planetary tidal dissipation factors

In this section, we study the past evolution of Kepler-91b to constrain the physics of tides and their dissipation. Before presenting our method, we should remind the reader of the different assumptions made during this study. First, on the stellar side, all the uncertainties presented here do not include possible uncertainties on the stellar modelling done in Sect. 3. Moreover, we ensure that the observed stellar rotation of ${\text{Ω}}_{*}={265}_{-16}^{+18}$ nHz is reproduced in our modelling at the age of the system. On the planetary side, we adopted the planetary characteristics presented in Table 4; each evolutionary track starts at the beginning of the stellar MS to avoid any uncertainties coming from the dynamical tides. The planet bond albedo was kept constant during the whole evolution with A = 0.30, and we also assumed that the planet was in pseudo-synchronisation throughout its evolution, meaning that the dissipation of the planetary equilibrium tides can be underestimated at the beginning of the evolution. Finally, we assumed that the planet’s orbital plane was co-planar with respect to the stellar equatorial plane and that the system is only composed of one planet during its whole lifetime.

To constrain the dissipation of the planet equilibrium tides, we built a grid of evolutionary tracks with different planetary equilibrium tides’ dissipation coefficients: k_2,p∆t_p. The properties and initial parameters used to compute the grid are resumed in Table B.1. The end of each evolutionary track corresponds to the age of the system given in Table 3. The purpose of this grid is to find the maximum eccentricity that Kepler-91b can reach at the age of the system given a tidal dissipation factor while fitting all the other planetary and stellar properties. The initial semi major axis is optimised in order to reproduce the observed period at the age of the system and for each k_2,p∆t_p. The maximum eccentricity reached by the system allows us to constrain the tidal dissipation of the planet by giving a limit to k_2,p∆t_p, after which point the observed eccentricity cannot be reproduced.

It is possible to notice that a degeneracy exists in between the initial parameters (orbital, stellar) and the value of the tidal dissipation factor in reproducing the current properties of the system. However, this degeneracy does not concern the tidal dissipation in the host star, which is computed by following the structural and rotational evolution of the host star, as in Eqs. (3) and (4). Therefore, while it is possible to tune the value of the dissipation factor for the planet, this is not the case for the tides dissipated in the host star. In that sense, the degeneracy on the stellar tidal dissipation factor is removed. There is a degeneracy in the initial orbital and planetary parameters and the planetary tidal dissipation factor, but not on the tidal dissipation factor of the host star. For this reason, we performed a study of the parameter space and found some constraints on k_2,p∆t_p. Instead of providing a set of degenerate initial parameters allowing us to reproduce the observations, we looked for the maximum value of the tidal dissipation factor that, by fine tuning of the rest of the parameters, can lead to the reproduction of the system properties.

We should mention that we are not trying to reproduce the observed eccentricity here. We are looking at the maximum eccentricity that the planet can have at the end of an evolutionary track while reproducing all the others orbital parameters (i.e. the period, planet radius, and stellar rotation). The initial eccentricities were set at the maximum values tolerated by the system to avoid any early engulfment, or to 0.90 to avoid any non-physical behaviour close to the singularity of e = 1.0. With this fixed initial eccentricity and each dissipation factor, we experimentally adjusted the initial semi-major axis to reproduce the planet period at the current age of the system. We define a limit track for the planet orbital parameters with a given dissipation factor.

At the age of the system, two possible scenarios can occur with respect to the final eccentricity value of the planet orbit:

If the eccentricity is higher than the observed one within the uncertainties, it is possible to adjust the initial semi-major axis and eccentricity to reproduce the observed planetary parameters.

If the eccentricity is lower than the observed one, it is not possible to find a combination of initial parameters to reproduce the observed orbital and planetary parameters; thus, the dissipation factor is not consistent with the observations.

Particular care was taken to reproduce the observed planetary radius to not bias the tidal dissipation. Finally, we made this grid for three different planetary masses, $M_{\text{pl}}={0.88}_{-0.33}^{+0.17}{M}_{\text{Jup}}$, to account for the uncertainties on the mass of the planet. The results of this model grid and the consequent constraints on the tidal dissipation coefficient are presented in Fig. 5 and Table B.1. Figure 5 shows that when taking into account both the uncertainties on the mass of the planet and on its observed eccentricity, the maximum value of k_2,p∆t_p allowed to reproduce the eccentricity is ~0.4 ± 0.25 s. For such values, the orbital migration undergone by the planet is quite significant, starting at an orbital distance of 0.33 AU with an eccentricity of 0.90; the planet finishes its evolution at 0.0721 AU with an eccentricity of 0.049. The damping of the eccentricity mostly results from the planet equilibrium tides during the early time of the system in this scenario.

Figure 5 also illustrates two regimes of past eccentricity damping, with the first one led by the planet equilibrium tides, where the past damping is proportional to 1/q (q = M_pl/M_*), as on the right side of Fig. 5. The second regime corresponds to a domination of the past damping by the stellar equilibrium tides. This regime occurs in the Kepler-91 system when k_2,p∆t_p ≲ 0.015 s. This transition in the regimes of dissipation can be seen with the transition in the dependency of q in Fig. 5, as stellar equilibrium tidal dissipation is proportional to q² (see Eqs. (3) and (4)), while the planet equilibrium tides dissipation is proportional to 1/q (see Eqs. (11) and (12)). Leconte et al. (2010) showed that for such systems, the planet tidal dissipation factor is of the order of k_2,p∆t_p ~ 10⁻²–10⁻³ s, meaning that in our case we expect the stellar equilibrium tides to have dominated the past evolution of Kepler 91b while planet equilibrium tides would have had a limited impact.

To express the dissipation of the planetary tides, we chose to use the constant time lag (CTL) model. However, one may choose to adopt the constant phase lag (CPL) model with the Q′ formalism (Goldreich 1963; Goldreich & Soter 1966). Staring from our Eqs. (11) and (12), we can estimate the dissipation in the CPL model by assuming that the planet is in pseudosynchronisation, as we did. In this case, the most dissipative tides are the annual eccentric tides, and we can link the dissipation of tides in the CPL model to their dissipation in the CTL model using (Leconte et al. 2010): $${Q^{\prime }}_{\text{p}}^{-1}\approx \frac{2}{2}{k}_{2,\text{p}}\text{Δ}{t}_{\text{p}}n,$$(25)

where Q′_p is the reduced tidal quality factor of the planet quantifying the efficiency of the tidal dissipation. Q′_p is calibrated depending on the type of planets studied and kept constant for the whole evolution. We should mention that the CPL model is expected to be accurate to e² (Leconte et al. 2010); therefore, one should be careful before using it for very eccentric orbits.

Equation 25 can be used to replace k_2,p∆t_p in Eqs. (11) and (12), with a new free parameter controlling the dissipation of tides: Q′_p. Repeating the method previously presented on a new grid of model with Q′_p, we obtain the equivalent of Fig. 5 in this new model (Fig. C.1).

All the remarks made concerning Fig. 5 can be applied to Fig. C.1. As accomplished previously, we can constrain the parameter controlling the planetary tidal dissipation. With our grid, we found that Q′_p should be higher than ${4.5}_{-1.5}^{+5.8}\times {10}^5$ to reproduce the observed eccentricity. Equivalently, a transition from a regime where the past tidal dissipation is dominated by planetary equilibrium tides towards a dissipation dominated by stellar equilibrium tides appears around Q′_p = 5.5 × 10⁶.

The minimum value of ${Q^{\prime }}_p={4.5}_{-1.5}^{+5.8}\times {10}^5$ found in this article is compatible with the common value of the dissipation factor for Jupiter found in literature (Q′_p = 10⁶; Goldreich 1963; Goldreich & Soter 1966; Leconte et al. 2010; Patra et al. 2017). We can mention that the value found is, however, not compatible with the tidal dissipation of the hot Jupiter HAT-P-13b (Batygin et al. 2009; Buhler et al. 2016), which is lower than Q′_p = 3 × 10⁵, found by Batygin et al. (2009) who neglected the impact of the stellar tidal dissipation. The results obtained here can also be compared to the empirical law to find the Q′_p provided in the work of Hansen (2010, 2012). When using this law by providing the properties of our system, we found a value for the dissipation factor Q′_p = 3.2 × 10⁷ that is consistent with our findings.

To illustrate the evolution of the Kepler-91 system in the different regimes of dissipation, in Fig. 6 we display the evolution of the semi-major axis and eccentricity of the planet with various values of the tidal dissipation factor k_2,p∆t_p used in Fig. 5 and a planet mass of M_pl = 0.88 M_Jup. Figure 6 highlights the significant difference between the regime where past dissipation was dominated by stellar tides with the black curve corresponding to k_2,p∆t_p = 0.01 s and the regime where past tides were dominated by planetary tides, where k_2,p∆t_p = 0.1 s and 0.4 s. In the regime dominated by stellar tides, the dissipation during the stellar MS is slow and the damping of the eccentricity almost null, even at high eccentricity. During the post-MS, the migration and damping eccentricity become particularly efficient at high eccentricity; this dissipation is attributed to stellar tides and is thus particularly dependent on the stellar parameter and the stellar model.

The opposite happens in the regime where past dissipation was dominated by planetary tides as illustrated in Fig. 6. The damping of eccentricity and migration is so efficient during the stellar MS due to planetary tides that most of the eccentricity is dissipated before the end of it. In such scenarios, stellar tides only have a limited impact on the evolution of the system during the stellar post-MS. Therefore, in this regime the stellar parameters only have a small impact on the final evolution of the system.

Fig. 4 Illustration of the impact of the equilibrium planetary and stellar tides on the orbital evolution of a system represented through the timescales of migration ${\tau }_{\text{a}}=\left(a/\dot{a}\right)$ and circularisation τe = (e/ė). The system evolves around the stellar model of Kepler 91 presented in Sect. 5. The initial mass of the planet is set at 0.88 Jupiter’s mass, its initial semi major axis is set at 0.10 AU, its initial eccentricity to 0.2 and k2,p∆tp is set at 10−2 s. The upper panel illustrates the different timescales, the dashed lines correspond to the circularisation timescales while the solid lines represent the migration timescales. Each modified component of these timescales (planetary, stellar and their total) are represented namely in green, red and black. The lower panel correspond to the associated migration of the planet, the engulfment of the planet occurs at tengulfment = 4.897 Gyr. Finally the shaded areas correspond to the different evolutionary phases of the star, in yellow the MS, in orange, the sub-Giant phase, and in red the RGB phase.

Fig. 5 Results from the grid of evolutionary tracks illustrating the maximum eccentricity that Kepler 91b can reach for given values of the tidal dissipation parameter k2,p∆tp. These results are compared to the observed eccentricity for different planetary masses considered. A summary of the initial and final stellar, orbital and planetary parameters of each evolutionary track used here is presented in Table B.l.

Fig. 6 Evolution of the semi major axis and eccentricity as a function of time on, respectively, the top and bottom panel. Each curve represents the path taken by the planet with the different k2,p∆tp used in Fig. 5. The blue shaded areas are to the corresponding observed quantities. The yellow area corresponds to the stellar MS, the orange one is the stellar sub-Giant phase, and the red one is the stellar RGB phase. The age of the system corresponds to the end of these evolutionary sequences.

5.3 Evolution of the dissipation of the stellar equilibrium tides

In this section, we investigate the tidal dissipation of Kepler-91 starting from the stellar model and the previously derived results. Usually, the dissipation of the stellar equilibrium tides is expressed with the same formalism as the planetary equilibrium tides, namely the CTL or CPL models. Each of these models are associated with a free parameter controlling the tidal dissipation, k_2,*∆t_* for the CTL model and Q′_* for the CPL model. In our modelling, the study of the dissipation of stellar equilibrium tides does not rely on a free parameter; on the contrary, the dissipation is deduced from dedicated stellar evolutionary models. We link the stellar tidal dissipation derived in our modelling to the classical dissipation models presented before. Bolmont & Mathis (2016) introduced a tidal dissipation factor, σ_*, similar to a k_2,*∆t_*, exclusively relying on the structure and rotation of the star. Following the approach of Rao et al. (2018), to derive σ_* we can start from Eq. (12) of Bolmont & Mathis (2016): $$\begin{array}{l}{\left(\frac{1}{a}\frac{\text{da}}{\text{dt}}\right)}_{\text{eq.tides,*}}=\hfill \\ 9{\sigma }_{*}{M}_{*}{R}_{*}q\left(1+q\right){\left(\frac{R_{*}}{a}\right)}^8\left[N_{\text{a}2}\left(e\right)\frac{{\text{Ω}}_{*}}{n}-N_{\text{a}1}\left(e\right)\right],\hfill \end{array}$$(26)

and we identify the tidal dissipation factor σ_* using the equivalent equation in our model, Eq. (3); σ_* can then be expressed as $${\sigma }_{*}=\frac{1}{9M_{*}{R}_{*}}\frac{f}{\tau }\frac{M_{\text{env}}}{M_{*}}.$$(27)

In our expression of the dissipation factor, σ_*, the tidal dissipation of the star depends on the orbital period of the planet trough the quantity f_orb defined in Eq. (10). This factor helps to account for the only convective cells that participate in the dissipation process in the regime of fast tides, and it depends on the position of the planetary companion at each evolutionary time step. Thanks to the constraint on the tides dissipated in the planet, its position is known at each time-step, and consequently it is also possible to compute the f_orb factor.

The tidal dissipation of the star, under the form of σ_*, can be linked to k_2,*∆t_* with the expression given in Bolmont & Mathis (2016): $${\sigma }_{*}=k_{2,*}\text{Δ}{t}_{*}\frac{2G}{3R_{*}^5},$$(28)

where G is the universal gravitational constant.

The constraints on the planetary tidal dissipation given in the previous section are used to constrain the factor f_orb to determine a range of values of σ_*, and then of k_2,*∆t_* using both Eqs. (27) and (28).

Finally, in the models used by Bolmont & Mathis (2016), σ_* is assumed constant through the evolution of the system for the stellar equilibrium tides. We also investigated the validity of this approximation with a similar model by setting σ_* to a constant value calibrated at the MS of Kepler-91.

The result from our modelling is illustrated in Fig. 7, where we represent the evolution of k_2,*∆t_* as a function of time for the different models considered. Using Eq. (25), we did the same study converting k_2,*∆t_* to Q′_*; the results are presented in Fig. D.1.

Figure 7 shows that a difference in tidal planetary dissipation and thus on the position of the planet for a given system only has a minor impact on the dissipation factor of the stellar equilibrium tides. We can mention, however, that in Fig. 7, the position of the planet dominates the variation of k_2,*∆t_* during the inward migration of the planet at the end of the stellar RGB phase.

As shown by Fig. 7, the dissipation model with a constant σ_* used by Bolmont & Mathis (2016) is a fair approximation during the stellar MS. However, such an approximation is inadequate during the PMS and post-MS phases, where changes of the tidal dissipation by orders of magnitude are obtained as a result of the evolution of the internal stellar structure (Rao et al. 2018). This illustrates the need to correctly follow the internal structure of the star during these evolutionary phases, which can either be done by computing full models of the star as in the present study or in our previous works (Pezzotti et al. 2021b,a, 2022; Bétrisey et al. 2022), or by using pre-calculated grids of stellar models that provide the needed quantities M_env and R_env as done for instance by Bolmont & Mathis (2016).

Fig. 7 Illustration of the evolution of k2,*∆t* as function of time for the different models considered. In orange our stellar model combined to the upper boundary of the planetary tidal dissipation found previously k2,p∆tp = 0.65. In black our stellar model is combined to the lower value of the planetary dissipation k2,p∆tp = 0. We can notice that the black and orange lines are close to be superimposed. In green is represented the model with σ* = const, as done in Bolmont & Mathis (2016), the constant was calibrated at value in the MS (σ* = 5.536 × 10−66 g−1 cm−2 s−1 at t = 0.83 Gyr). Finally the shaded areas correspond to the different evolutionary phases of the star, in yellow the MS, in orange, the sub-Giant phase, and in red the RGB phase.

6 Conclusion

In this study, we carried out a detailed analysis of the properties of Kepler-91, a well-known RGB exoplanet host star observed by Kepler (Lillo-Box et al. 2014a,b). Our goal was to couple the modelled evolution of the close-in planet to the evolution of the star to study and constrain the dissipation of the planetary equilibrium tides.

The strategy of our study is based on a detailed modelling of the entire system. First, the star was modelled independently with the Code Liégeois d’Evolution Stellaire (CLES, Scuflaire et al. 2008b) and the Liège OScillation Code (LOSC, Scuflaire et al. 2008a), taking particular care to reproduce the observed classical constrains and asteroseismic quantities. Then, the evolution of the star was coupled with the evolution of the planet by our orbital evolution code (Privitera et al. 2016a,b; Meynet et al. 2017; Rao et al. 2018). Our orbital evolution code was modified in order to account for the evolution of the eccentricity of the planet aiming at precisely describing the orbital evolution of Kepler 91b in the past and future. In particular, the tidal dissipation of both the star and the planet were taken into account to study the evolution of the planet eccentricity and orbital period.

In the modelling of the planetary equilibrium tides, the tidal dissipation is proportional to a free parameter, a dissipation factor, and either k_2,p∆t_p (Darwin 1879) or Q′_p (Goldreich 1963; Goldreich & Soter 1966) depending on the assumptions made.

The principle of our study was to determine the maximum eccentricity reachable by the planet, reproducing all the observational constraints and given a tidal dissipation factor.

With this methodology, we found that k_2,p∆t_p should be lower than 0.4 ± 0.25 s to be able to reproduce the observed eccentricity of the planet. We also found that Q′_p should be greater than ${4.5}_{-1.5}^{+5.8}\times {10}^5$. We then looked at the dissipation of the stellar equilibrium tide, which was computed based on the physical description of Zahn (1966, 1977, 1978, 1989); Villaver & Livio (2009) and Villaver et al. (2014) with the stellar structure given by our stellar modelling. The tidal dissipation from our modelling was compared to a more classical formalism such as the σ_* formalism (Bolmont & Mathis 2016). We confirm previous results of Rao et al. (2018) that a constant σ_* is a correct approximation during the stellar MS, but not during the PMS and post-MS phases, where the evolution of the internal stellar structure needs to be followed in detail.

More than putting a constraint on the planetary dissipation of the planet equilibrium tides, our study shows that hot Jupiters have tidal dissipation coefficients in the same range of values as colder, Jupiter-like planets. This result can be important for the theoretical study of multi-planetary systems as tidal interactions between planets play a key role in the understanding of their future evolution and stability. Moreover, our results show that planetary structural modifications caused by the extended atmosphere of hot Jupiters do not significantly increase the tidal dissipation factors.

The main limits to the present method come from the overall precision of the characterisation of the system, and in particular the planetary properties. In the future, the same type of study could be applied to similar targets characterised with more accuracy. The ideal target should be orbiting an evolved star with observed oscillation modes to conduct a similar precise stellar modelling; thus, we would favour a Kepler (Borucki et al. 2010) or TESS (Ricker et al. 2014, 2015) target. The studied planet needs to be in close orbit with respect to the star and have a high eccentricity; in this aspect, the target presented in Grunblatt et al. (2018, 2022) could be a good choice. Finally, the characterisation of the planet, and in particular a precise determination of the planetary mass, is necessary. The use of other constraints such as the planet period variation or precession variation can also provide important information about the tidal dissipation of the planet.

Acknowledgements

The authors are thanking the anonymous referee for his constructive comments. L.F. was supported by the Fonds de la Recherche Scientifique-FNRS as a Research Fellow. C.P. acknowledges support by the Swiss National Science Foundation (project number 200020-205154). P.E. has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 833925, project STAREX). G.B. acknowledges funding from the SNF AMBIZIONE grant No. 185805 (Seismic inversions and modelling of transport processes in stars). This article used an adapted version of InversionKit, a software developed within the HELAS and SPACEINN networks, funded by the European Commission’s Sixth and Seventh Framework Programmes. This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. The computations were performed at University of Geneva on the Baobab and Yggdrasil clusters. This research has made use of NASA’s Astrophysics Data System.

Appendix A Seismic data

Appendix B Additional tables

Appendix C Constraints on the tidal dissipation with the CPL model

Fig. C.1 Results from the grid of evolutionnary tracks where the maximum eccentricity of Kepler 91b is given for some values of the tidal dissipation parameter k2,pΔtp. These results are compared to the observed eccentricity for the three different planetary masses considered.

Appendix D Additional figures

Fig. D.1 Illustration of the evolution of Q′* as function of time for the different models considered. In orange our stellar model combined to the upper boundary of the planetary tidal dissipation found previously Q′* = 3 × 105. In black our stellar model is combined to the lower value of the planetary dissipation Q′* → ∞. We can notice that the black and orange line are close to be superimposed. In green is represented the model with σ* = const, as done in Bolmont & Mathis (2016), the constant was calibrated at value in the MS (σ* = 5.536 × 10−66 at t = 0.83 Gyr). Finally the shaded areas correspond to the different evolutionary phases of the star, in yellow the MS, in orange, the sub-Giant phase, and in red the RGB phase.

References

All Tables

All Figures

Fig. 1 Echelle diagram of Kepler 91 comparing observed oscillation frequencies (listed in Table A.1) represented as squared points and theoretical oscillation frequencies represented as circle points. The whole modelled oscillation spectrum for the ℓ = 1 modes is included with shaded points. The difference between observed and modelled data is quantified, in the legend, as a χ2 for each spherical order ℓ. Finally, the error bars on the observed frequencies are included in the data points and not clearly distinguishable here.

In the text

	Fig. 2 Hertzsprung-Russell diagram of the evolution of Kepler 91 from the pre-Main Sequence to the modelled actual state of the star and the its future evolution with dotted line. The grey area corresponds to the constraints on the effective temperature and luminosity listed in Table 1.
In the text

	Fig. 3 Diagram underlying the geometric increase of irradiation of the host star onto the planet surface caused by the increase of the stellar radius.
In the text

Fig. 4 Illustration of the impact of the equilibrium planetary and stellar tides on the orbital evolution of a system represented through the timescales of migration ${\tau }_{\text{a}}=\left(a/\dot{a}\right)$ and circularisation τe = (e/ė). The system evolves around the stellar model of Kepler 91 presented in Sect. 5. The initial mass of the planet is set at 0.88 Jupiter’s mass, its initial semi major axis is set at 0.10 AU, its initial eccentricity to 0.2 and k2,p∆tp is set at 10−2 s. The upper panel illustrates the different timescales, the dashed lines correspond to the circularisation timescales while the solid lines represent the migration timescales. Each modified component of these timescales (planetary, stellar and their total) are represented namely in green, red and black. The lower panel correspond to the associated migration of the planet, the engulfment of the planet occurs at tengulfment = 4.897 Gyr. Finally the shaded areas correspond to the different evolutionary phases of the star, in yellow the MS, in orange, the sub-Giant phase, and in red the RGB phase.

In the text

Fig. 5 Results from the grid of evolutionary tracks illustrating the maximum eccentricity that Kepler 91b can reach for given values of the tidal dissipation parameter k2,p∆tp. These results are compared to the observed eccentricity for different planetary masses considered. A summary of the initial and final stellar, orbital and planetary parameters of each evolutionary track used here is presented in Table B.l.

In the text

Fig. 6 Evolution of the semi major axis and eccentricity as a function of time on, respectively, the top and bottom panel. Each curve represents the path taken by the planet with the different k2,p∆tp used in Fig. 5. The blue shaded areas are to the corresponding observed quantities. The yellow area corresponds to the stellar MS, the orange one is the stellar sub-Giant phase, and the red one is the stellar RGB phase. The age of the system corresponds to the end of these evolutionary sequences.

In the text

Fig. 7 Illustration of the evolution of k2,*∆t* as function of time for the different models considered. In orange our stellar model combined to the upper boundary of the planetary tidal dissipation found previously k2,p∆tp = 0.65. In black our stellar model is combined to the lower value of the planetary dissipation k2,p∆tp = 0. We can notice that the black and orange lines are close to be superimposed. In green is represented the model with σ* = const, as done in Bolmont & Mathis (2016), the constant was calibrated at value in the MS (σ* = 5.536 × 10−66 g−1 cm−2 s−1 at t = 0.83 Gyr). Finally the shaded areas correspond to the different evolutionary phases of the star, in yellow the MS, in orange, the sub-Giant phase, and in red the RGB phase.

In the text

	Fig. C.1 Results from the grid of evolutionnary tracks where the maximum eccentricity of Kepler 91b is given for some values of the tidal dissipation parameter k2,pΔtp. These results are compared to the observed eccentricity for the three different planetary masses considered.
In the text

Fig. D.1 Illustration of the evolution of Q′* as function of time for the different models considered. In orange our stellar model combined to the upper boundary of the planetary tidal dissipation found previously Q′* = 3 × 105. In black our stellar model is combined to the lower value of the planetary dissipation Q′* → ∞. We can notice that the black and orange line are close to be superimposed. In green is represented the model with σ* = const, as done in Bolmont & Mathis (2016), the constant was calibrated at value in the MS (σ* = 5.536 × 10−66 at t = 0.83 Gyr). Finally the shaded areas correspond to the different evolutionary phases of the star, in yellow the MS, in orange, the sub-Giant phase, and in red the RGB phase.

In the text