Transitory tidal heating events and their impact on cluster isochrones

We perform a calculation from first principles (Moreno et al., 2011; Koenigsberger et al., 2021; Estrella-Trujillo et al., 2023). It consists of simultaneously solving the orbital motion and the equations of motion of a 3D grid of volume elements that cover the central region of the tidally perturbed star, which we refer to as the primary. This central region, referred to as the core, rotates as a solid body with constant angular velocity $\omega_{0}$ in an inertial reference frame, and its rotation axis is perpendicular to the orbital plane. It does not coincide with the central nuclear burning region, and in all our models it includes a significant portion of the radiative envelope. The equations of motion are solved in the reference frame with origin in the center of the primary and rotating at the rate of the companion’s orbital motion, and include gravitational, centrifugal, Coriolis, gas pressure gradient and viscous forces. A seventh order Runge-Kutta integrator is used to solve the set of equations. The companion is considered to be a point-mass source and its orbital plane is coplanar with the primary star’s equator. The code, in its latest version, is named TIDES-nvv and is publicly available.¹¹1Moreno, E., & Koenigsberger, G. 2024, Tidal Interactions with Dissipation of Energy due to Shear version nvv, Zenodo. https://doi.org/10.5281/zenodo.10799022

The rate of energy dissipation per unit volume, $\dot{E}_{\mathrm{V}}$ depends on the mass density $\rho$, viscosity $\nu$ as well as the velocity gradients in the radial, latitudinal and longitudinal directions across neighboring volume elements (McQuarrie, 1976). We use the formulation given in Moreno et al. (2011), in which only the gradients in angular velocity are implemented. Also, because the azimuthal perturbations are an order of magnitude larger than those in the radial and polar direction (Scharlemann, 1981; Harrington et al., 2009), the calculations performed in this paper only consider these perturbations.

The viscosity is non-isotropic and time-dependent. It is assumed to be dominated by a turbulent term whose value is computed for each volume element interface following the formulation of Landau & Lifshitz (1987):

The TIDES calculation reproduces the oscillatory properties of the tidal flows and the evolving differential rotation structure that is established as angular momentum is transported between stellar layers in response to the tidal torque. The time-marching TIDES algorithm is valid for binary stars with arbitrary rotation velocity and eccentricity, as long as neighboring grid elements retain contact over at least $\sim$80% of their surface and the centers of mass of two adjoining grid elements do not overlap. However, it neglects buoyancy effects, heat and radiation transfer, fluid dynamics microphysics, diffusion and advection. Also, because the energy dissipation rates are not fed back into the system, the model does not compute the modified internal structure that would result from this feedback. Hence, it is valid only for studying the short-term (i.e., orbital timescales) behavior under the influence of tides.

The input parameters are described in Table 1. The orbital period, eccentricity, primary and companion masses are labelled $P$ and $e$, $m_{1}$ and $m_{2}$, respectively. The primary’s structure is assumed to be polytropic with an index $n$. Its initial unperturbed radius is $R_{1}$, and its initial rotation structure is that of uniform rotation with an angular velocity $\omega_{0}$ in the inertial reference frame. We define the parameter $\beta_{0}=\omega_{0}/\Omega_{0}$ which establishes the initial rotation state. Here, $\Omega_{0}$ is the orbital angular velocity at a reference point in the orbit. In circular orbits, $\beta_{0}$ is constant but for $e\neq 0$ it varies over the orbital cycle, in which case, $\Omega_{0}$ is chosen to be the orbital angular velocity at the time of periastron. In general, the layers above the core respond to the forces in the system and over time develop a differential rotation structure. The oscillatory tidal field is superposed on the differential rotation structure (see, for example, Fig. 5 in Koenigsberger et al., 2021). The only case in which such a differential rotation structure does not develop is when $e=0$ and $\beta_{0}=1$, that is, in circular orbits in which the stellar rotation rate equals the orbital rate.

The computational input parameters are the grid size ($N_{r}\times N_{\varphi}\times N_{\theta}$), the thickness of each layer ($\Delta R/R_{1}$), the number of orbital cycles over which the computation is run ($N_{\mathrm{cycle}}$) and the tolerance for the Runge-Kutta integration. The triad $(r,\varphi,\theta)$ correspond, respectively, to the distance from the stellar center, the azimuth angle and the colatitude. The angle $\varphi=0$ is measured from the line joining the centers of the two stars in the direction of the orbit.²²2Note that because the solution of the equations of motion is performed in the reference frame rotating with the binary orbit ($S^{\prime}$ in the notation of Moreno et al. 2011), these variables correspond to the primed variables but for the purposes of this paper, this notation is not required and thus omitted.

For the nominal calculations presented in this paper, we chose the primary mass $m_{1}=1.2\,M_{\odot}$, which is slightly smaller than the M67 turnoff mass, $m_{\mathrm{MSTO}}=1.26\,M_{\odot}$. The orbital period is 1.44 d, the same as in Estrella-Trujillo et al. (2023). We probed stars having radii $R_{1}/R_{\odot}=0.97$, 1.2, 1.32 and 1.50, corresponding to different evolutionary times. The primary’s rotation velocity was assumed to always be subsynchronous with $\beta_{0}=0.7$, which corresponds to ${\mathrm{v}}_{\mathrm{rot}}/\mbox{(km\,s${}^{-1}$)}=24$, 29, 32, 36 for the above listed radii, respectively. These values are significantly larger than typical projected rotation speeds in M67 which lie in the range ${\mathrm{v}}_{\mathrm{rot}}\sin\mathrm{i}\sim$5–8 km s^-1 (Melo et al., 2001) but these cited values likely correspond to mostly single stars. Short-period binaries are expected to rotate much faster because synchronization times ($\sim 10^{7}$ yr) are orders of magnitude smaller than the main-sequence lifetime. Furthermore, Nine et al. (2020) find 113 stars with $\mathrm{v}\sin\mathrm{i}>120$ km s^-1 in NGC 7789. We note in passing that these very rapid rotators are hot and blue.

Energy dissipation rates were computed for $m_{2}/M_{\odot}=0.8$ and 0.4 companions while holding the other parameters constant. These models are labelled Case 1 and Case 2, respectively in Table 2. Each case is subdivided into cases b-d, corresponding to the different adopted values of $R_{1}$, which are listed in Column 3 of this table. We also computed a set of models with $m_{2}=0.4\,M_{\odot}$ and $\beta_{0}=1.1$, which is a super-synchronously rotating model. This is Case 3.

The energy dissipation rate per unit volume $\dot{E}_{V}$ is a time-dependent quantity due to the tidal oscillations and, in eccentric orbits, due to the varying conditions between periastron and apastron. Thus, $\dot{E}_{\mathrm{V}}=\dot{E}(r,\theta,\varphi,t)$. The algorithm computes this quantity at a set of user-defined orbital phases $t_{i}$ within several orbital cycles that are also chosen by the user. Because MESA performs only a one-dimensional structure calculation, the spatial distribution of the tidal perturbations is not relevant in this paper. Thus, we use the integral of $\dot{E}_{\mathrm{V}}$ over latitude and longitude for the energy dissipation rate profile in the radial direction:

In general, we chose to output data at $N_{\mathrm{tt}}=40$ orbital phases within each of five orbital cycles. An inspection of the five cycles allows assessment of the long-term variability patterns and whether the initial transitory state has ended. After the transitory state, the average angular velocities evolve very slowly as the star tries to synchronize its rotation rate with the orbital angular velocity. For the current MESA calculations we chose to use the last orbital cycle of the calculation and averaged $\dot{E}(r,t)$ over the 40 orbital phases of the cycle to obtain the radial energy dissipation rate profile, $\dot{E}(r)$, to be used as input for the MESA models:

Figure 1 shows that $\dot{E}(r)$ displays a pronounced increase between the lower boundary of the calculation (0.4 $R_{\odot}$) and the surface. The rate of increase depends primarily on the degree of departure from synchronicity, $|\beta_{0}-1|$, and less so on the mass of the companion. This can be seen by comparing Cases 1 and 2 ($m_{2}=0.8$ and $m_{2}=0.4$, respectively), which differ in $\dot{E}(r)$ by a factor $\sim 3$, while Cases 2 and 3 (both with $m_{2}=0.4$, but with $\beta_{0}=0.7$ and $\beta_{0}=1.1$, respectively) differ by $\sim$two orders of magnitude. Also noteworthy is the fact that $\dot{E}(r)$ is larger in the $\beta_{0}=0.7$ model than the corresponding $\beta_{0}=1.1$ model even though the former is rotating more slowly than the latter; for example, ${\mathrm{v}}_{\mathrm{rot}}=32$ km s^-1 for $\beta_{0}$=0.7 while ${\mathrm{v}}_{\mathrm{rot}}=51$ km s^-1 for $\beta_{0}=1.1$.

There are several factors that can affect the energy dissipation rates obtained for any given model. The model assumes that the central part of the star, below the smallest radius adopted for the TIDES calculation, is in solid body rotation. Even though the tidal perturbations in this region induce negligible velocity perturbations, the possibility exists that it has an intrinsic differential rotation. This would add shear energy dissipation to these inner zones and potentially affect the stellar structure as described in the next section. However, asteroseismological studies of pulsating stars in the mass-range 1.4-5 $M_{\odot}$ have led to the conclusion that, for rotation velocities up to 50% of breakup, the stars deviate from solid body rotation only mildly (Aerts et al., 2017). Another factor concerns the viscosity. Our calculations were performed with $\lambda=1$, which most likely is an overestimate. However, the fact that $\dot{E}_{\mathrm{V}}$ scales approximately linearly with $\nu_{\mathrm{turb}}$ which, in turn, scales linearly with $\lambda$ in the above prescription, allows scaling the values for smaller $\lambda$ values.

We use the open-source stellar structure and evolution code MESA, version 15140³³3https://docs.mesastar.org/en/r15140/, (Paxton et al., 2011, 2013, 2015, 2018, 2019) to study the effect of tidal heating on the stellar envelope and resulting observable properties. The basic model physics parameters (i.e., abundances, equation of state, opacities, reaction rates, diffusion, convection, overshooting, semiconvection, thermohaline, boundary conditions) are the same as those used by the MIST project (see Table 1 of Choi et al., 2016) for low-mass stars, which were calibrated using the old open cluster M67. We made the adjustments necessary to update the input file commands from MESA v7503 used by MIST (Dotter, 2016; Choi et al., 2016) to v15140 used in this paper. We do not include rotation in our models.

Stellar models were run for low-mass stars with masses close to the turn-off mass (around 1.26 $M_{\odot}$) for M67, which we use as a reference and source of observational comparisons. The energy dissipation rates as a function of radius, $\dot{E}(r)$, calculated by TIDES are used as inputs to the other_energy hook in the MESA code and allow us to calculate the effect of energy dissipation on the stellar envelope of the model stars. In particular, as our fiducial energy input, we chose the Case 1 energy dissipation profile (i.e., $m_{2}=0.8\,M_{\odot}$, $\beta_{0}=0.7$), which is shown in Fig. 1 and listed in Table 5 (see Appendix).

The TIDES tidal perturbation model is not an evolution code and the outer radius of a given model is fixed. TIDES models with different radii can therefore be taken to represent possible points along the evolutionary track of a star. At any given time, the radius of a MESA stellar evolution model can be compared with the list of TIDES model radii. Interpolation between TIDES models of different radii enables us to find the appropriate energy dissipation rates for each MESA model in an evolutionary sequence. The TIDES shell radial positions are at fixed fractions of the outer radius of the star, and so this interpolation is performed using the normalized TIDES shell radii. The mass of a given shell is the same between TIDES models and can be easily found from the appropriate polytropic models. These masses are used to calculate the final energy dissipation profile per unit mass, which is then mapped back onto the MESA grid. MESA uses a much finer grid of interior points than TIDES, and so linear interpolation is used to provide values for the energy dissipation rate per unit mass at intermediate positions. The heating profiles for MESA models that have grown beyond the radius of the largest TIDES model are taken to be that of $R_{1}=1.5\,R_{\odot}$.

As a reference, we compute the evolution of a $1.2\,M_{\odot}$ star with no heating using the MIST set-up detailed in Choi et al. (2016). This model, which we refer to as the “standard model”, was first evolved from the birthline to the zero-age main sequence (ZAMS), which it reaches after $2.6\times 10^{7}$ yr. The ZAMS model was then used as the starting point for all other models with the same stellar mass. The standard model was evolved for 5 Gyr, by which time the star has left the main sequence and has entered the subgiant stage. This is far longer than the age of M67, which is estimated to be around 4 Gyr, but allows us to study the effect of envelope heating on more than one stage of evolution of the same star.

We examined two general scenarios:

1. 1.

   Scenario 1: Tidal heating is present throughout the main sequence. This scenario is possible only if the main sequence lifetime is on the order of the synchronization timescale. Given the long lifetime of stars in the mass range we analyze here, this will not occur unless the binary is forced to remain in an eccentric orbit by external gravitational forces – for example a third component in the system.

2. 2.

   Scenario 2: Tidal heating occurs in a short burst ($10^{7}$ yr) at a time when the internal structure undergoes a relatively fast transition, such as when the primary star reaches the end of its core hydrogen-burning phase, and the binary may become asynchronous.

For Scenario 1, several tests were performed and the resulting models are labeled as follows:

• •

  h00: this is the standard model. The star evolves from the ZAMS to the stopping point with no added heat;

• •

  h10: after age $10^{8}$ yr heat is injected steadily using the radial profile of the energy dissipation rate for the 1.2$+$0.8 $M_{\odot}$ system shown in Fig. 1.

• •

  h05 and h01: same as h10 but only 50% and 10%, respectively, of the radial energy dissipation rates are injected into the MESA model. These tests provide insight into the range of heating rates that have a significant effect on the stellar structure.

• •

  hk10: same as h10 except that the heat is not injected into the outermost 10% of the star. The reason for this test is that very near the surface, where the optical depth $\tau\leq 2/3$, Zahn (1975) suggested that a significant amount of the excess energy could be rapidly radiated away. Also, this test serves to examine the impact of the large heating rates in the surface layers on the internal structure.

Scenario 2 is investigated for stars with 1.16, 1.2, 1.26, 1.29 and 1.32 $M_{\odot}$ with the aim of making direct comparisons with observed stars in a range of evolutionary states in the M67 open cluster. In this scenario, MESA stellar evolution models for each mass were run without any heating until an age of 3.97 Gyr at which time the h10 heating rate was injected for a timespan of $10^{7}$ yr. This timescale is so short compared to the stellar lifetime that it can be viewed as a heating pulse, as we refer to it in what follows.

For the heating pulse, we set up the calculation to have a maximum timestep of $10^{5}$ yr, which is much finer than our standard MESA model runs, in order to capture the details of this transient episode in the life of the star. The computations were halted at an age of 4.25 Gyr but the main focus is around 3.98 Gyr, which the MIST isochrone shows to be a good fit to the observed stars in M67.

Figure 2 shows the manner in which the tidal heating affects the basic parameters of the perturbed star compared to the standard model. The heated star grows to a larger radius and this bloating scales with the amount of injected energy. Similarly, the effective temperature and surface luminosity are also larger than in the standard model. This description applies to all ages $\leq$5 Gyr, which is when the MESA model was terminated. The model with the highest heating rate (model h10) develops oscillations in radius, effective temperature and luminosity. The onset occurs around 3 Gyr, which is when the stellar radius has exceeded 1.5 $R_{\odot}$ and the heating rate becomes saturated at the limit of the largest TIDES model. Such instabilities are not present in model hk10, in which tidal heating is injected only in layers below the surface, $r\leq 0.85R_{1}$, where $R_{1}$ is the surface radius. This case could apply to a star that has already synchronized the outer layers with the binary orbit but not those below.

From an observational standpoint, eclipsing binaries are routinely used to determine the mass and radius of the two stars which, in turn, are used to constrain the age of the system. However, if tidal heating is active, stars of a given radius can have different ages. For example, a star evolving according to the standard model attains a radius $R_{1}=1.5\,R_{\odot}$ at an age of 4.1 Gyr while our h10 heated model attains this radius at an age of 3.1 Gyr (see Table 3). The internal structure of stars with the same radius but different heating rates will also be different. For example, Figure 3 presents a snapshot of the internal structure of the models at the point where the radius has expanded to 1.5 $R_{\odot}$. As was seen in Estrella-Trujillo et al. (2023), the internal luminosity of heated models is lower than that of the standard model for $r\leq 1.3\,R_{\odot}$ because the heated stars are younger. This is because the luminosity leaving the core increases with decreasing central hydrogen mass fraction, i.e. due to chemical abundance changes with age as a result of fusion reactions. On the other hand, the luminosity of the heated models rises rapidly in the outer layers because of the injected heat, while that of the standard model remains constant between the edge of the core and the stellar surface. The internal temperature and opacity do not change significantly.

The evolutionary tracks of the tidally heated models are compared to that of the standard model in Fig. 4. The heated models are characterized by higher temperatures and luminosities than the reference case. The difference between heated and standard models increases with increasing heating rates. The largest heating rate (model h10) presents large oscillations shortly after the star reaches a radius 1.5 $R_{\odot}$. Since MESA dynamically adjusts the integration timestep, these oscillations do not appear to be due to numerical instabilities. Instead, we tentatively interpret the oscillations as the outer layers of the star adjusting to the heating rate becoming saturated in the limit of the largest TIDES model ($R_{1}=1.5\,R_{\odot}$). A dynamical treatment of this model, beyond the scope of the present paper, would be able to ascertain whether such oscillations could lead to mass loss.

We note that the excess luminosity of the MESA models is significantly smaller than the total energy dissipation rate in the shearing layers; that is, $L_{\mathrm{exc}}=\alpha\dot{E}$, where $L_{\mathrm{exc}}=L_{\mathrm{heated}}-L_{\mathrm{standard}}$ is the difference between the MESA heated and standard models for the same radius and $\alpha<1$ is the fraction of dissipated energy that escapes as radiation and contributes to the luminosity. Comparing the results of model h00 with the heated models we find that $\alpha\leq 10\%$.⁴⁴4For example, at a time when the stellar radius reaches 1.5 $R_{\odot}$, the luminosity of the MESA models are, respectively, 3.01 $L_{\odot}$ (standard model), 3.29 $L_{\odot}$ (h01), 3.65 $L_{\odot}$ (h05), and 3.83 $L_{\odot}$ (h10) (see Table 3). Thus, $L_{\mathrm{exc}}=0.28$, 0.64 and 0.82 $L_{\odot}$, respectively. The tidal shear energy dissipation rates that were injected into the MESA model are 2.6, 6, and 21 $L_{\odot}$, respectively. The remaining $\sim$90% of the shear energy goes into sustaining the star against self-gravity and inflating it.

The evolutionary paths shown in Fig. 4 correspond to continuous heating, a rather unrealistic state because the main sequence lifetime of solar-type stars is $>10^{9}$ yr, while synchronization and orbital circularization processes are believed to last $\sim 10^{7}$–$10^{8}$ yr. Hence, even if a star reaches the main sequence with non-synchronous rotation or an eccentric orbit, it is expected to attain an equilibrium state early in its main sequence lifetime. A more realistic situation is one in which the system departs from the equilibrium condition at some point during its evolutionary trajectory. A clear example is when a star nears the end of core hydrogen burning and the core contracts and spins up while the envelope expands and spins down. At this stage, the rotation becomes asynchronous and a differential rotation structure is established, thus triggering a tidal heating event. This is described by our Scenario 2, the pulsed-heating scenario.

The evolution over time of luminosity and radius of the pulse-heated models are summarized in Fig. 5. As soon as the heat is injected, the stellar response is a radius and luminosity increase. The response time of the 1.20 and 1.26 $M_{\odot}$ models is $\sim 10^{5}$ yr, with the luminosity and radius first increasing and then decreasing to an intermediate level where they settle at a plateau. At this time they undergo a series of oscillations about the mean level. As soon as the heating is turned off, the luminosity and radius decrease to that of the standard model on a timescale that is as short as when the heat was first injected. For the 1.20 and $1.26M_{\odot}$ models, the mean luminosity and radius increases are modest, reaching a factor 1.5 increase in luminosity and only a factor 1.15 in radius. Also shown in Fig. 5 is the response of the $1.32M_{\odot}$ model, which shares similar features with the lower-mass models but without the oscillations. In addition, in this higher mass model, the luminosity increases by a factor $\sim 7.7$, while the radius grows by a factor $\sim 3.8$ over the $10^{7}$ yr duration of the heat injection.

Also shown in Fig. 5 is the response of the 1.32 $M_{\odot}$ model, which shares similar features with the lower-mass models but, instead of remaining at a plateau with oscillations, it continues brightening slowly until the heating is turned off. The overall luminosity increase is nearly a factor of three, and the radius increase mimics that of the luminosity, nearly doubling in size over the $10^{7}$ yr duration of the heat injection.

Estrella-Trujillo et al. (2023) found that the evolutionary track of a star in which the tidal heating is introduced late during its main sequence lifetime rapidly reaches the evolutionary track of a star with continuous heating throughout the main sequence. We confirm this result in our current calculations. Figure 6 shows that the evolutionary track of the pulse-heated model intersects that of of the h10 continuously heated model during the duration of the pulse. Thus, the continuously heated track serves to represent the locus of locations on the HRD that will be occupied during a tidal heating event. The precise location will depend on the age at which the heating event is triggered.

The theoretical HRD for masses in the range 1.16 $M_{\odot}$ and 1.29 $M_{\odot}$ with heat injection starting at an age of 3.97 Gyr and ending in 3.98 Gyr is displayed in Fig. 7. In all cases, the star crosses through the standard tracks of more massive stars during the heat pulse.

Another important feature of the pulses is that as soon as the heating is turned off, the star rapidly returns to the standard evolutionary track. Thus, the permanence of a tidally heated star on the blue side of the standard evolutionary track is as brief as the the time during which the heating is injected. Considering the long lives of these types of stars, the heating event can be considered to be extremely brief, thus explaining the relatively small number of stars that would be in this state at any given time.

Further insight into the timescales involved can be gained by a detailed look at the 1.32 $M_{\odot}$ model undergoing a tidal heating episode that lasts 10 Myr, as illustrated in Fig. 8. The heat is first injected at an age of 3.97 Gyr and causes an immediate, transitory spike in luminosity and surface temperature, but the star then settles onto a new trajectory higher up the HRD with increasing luminosity and decreasing temperature. The first $\sim$1 Myr of this trajectory exhibits oscillatory behaviour superimposed on the general trend and by the end of this stage (point 3.9719 in the figure) the surface temperature has returned to its original value but the luminosity is four times higher. For the remainder of the heat-pulse duration, $\sim 8.9$ Myr, the star’s luminosity continues to rise and the temperature continues to decrease, until the trajectory intercepts the red-giant branch shortly before the heating switches off (point 3.980). Once the heating is turned off, the star quickly descends the red-giant branch, decreasing in luminosity while its surface temperature remains roughly constant. It takes $\sim$8 Myr for the star to return to its original location in the HRD, with the luminosity decreasing from $45L_{\odot}$ to $8L_{\odot}$ in the first $\sim$1 Myr but taking a further 7 Myr to fall to $7L_{\odot}$. At this point the star retakes the original trajectory and the luminosity begins to increase again as the star climbs back up the red-giant branch. The increases (decreases) in the luminosity due to the heat pulse are accompanied by increases (decreases) in the stellar radius.

The appearance of oscillations in the heated 1.2 $M_{\odot}$ and 1.26 $M_{\odot}$ models merits further comment. The layers into which the heat is injected lie very far from the nuclear processing zone which, for the ages considered, is limited to the central region of the star. They do, however, include the surface convection zone. Kippenhahn diagrams of the models show that the oscillations are present at the base of the convection zone, within the radiative zone (see Fig. 11). These oscillations are forced and for this reason only last the duration of the heat pulse. When heat is added to the outer layers of the model stars, these layers expand and their densities fall. This reduces the opacity in these layers (which correspond to temperatures $T<10^{6}$ K, hence the opacities are Kramers-type) and so the radiative region of the star now extends to larger radii. The partially ionized zones of H and He are also lifted to larger radii. Although physical oscillations may occur under such circumstances, we do not model them since our minimum time step is of order 1000 years, far longer than any physical oscillation. The oscillatory behaviour exhibited in Figure 5 is a result of fluctuations in the forcing term due to the mapping between the stellar model and the TIDES model. This affects the base of the heating region, since its position depends on the total radius of the star, and as a consequence all the layers above receive different heat forcing at each timestep, and in particular the base of the convection zone. The behaviour appears oscillatory because the mapped heating at each position diminishes as the stellar radius grows, i.e. the forcing becomes less and this leads to damping. Determining whether true oscillations during a heating pulse can be excited requires a separate investigation which is beyond the scope of this paper.