Issue 
A&A
Volume 670, February 2023



Article Number  A44  
Number of page(s)  16  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202244971  
Published online  02 February 2023 
Structure and evolution of a tidally heated star
^{1}
Instituto de Ciencias Físicas, Universidad Nacional Autónoma de México, Ave. Universidad s/n, Chamilpa, 62210 Cuernavaca, Mexico
email: dtrujillo@icf.unam.mx; gloria@astro.unam.mx
^{2}
Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, Antigua Carretera a Pátzcuaro #8701, ExHda. San José de la Huerta, Morelia, CP 58089 Michoacán, Mexico
email: j.arthur@irya.unam.mx
^{3}
Instituto de Astronomía, Universidad Nacional Autónoma de México, Apdo. Postal 70264, Ciudad de México 04510, Mexico
Received:
13
September
2022
Accepted:
29
November
2022
Context. The shearing motion of tidal flows that are excited in nonequilibrium binary stars transform kinetic energy into heat via a process referred to as tidal heating.
Aims. We aim to explore the way tidal heating affects the stellar structure.
Methods. We used the TIDES code, which solves the equations of motion of the threedimensional (3D) grid of volume elements that conform multiple layers of a rotating binary star to obtain an instantaneous value for the angular velocity, ω″, as a function of position in the presence of gravitational, centrifugal, Coriolis, gas pressure, and viscous forces. The released energy, Ė, was computed using a prescription for turbulent viscosity that depends on the instantaneous velocity gradients. The Ė values for each radius were injected into a MESA stellar structure calculation. The method is illustrated for a 1.0 + 0.8 M_{⊙} binary system, with an orbital period of P = 1.44 d and departures from synchronous rotation of 5% and 10%.
Results. Heated models have a larger radius and surface luminosity, a smaller surface convection zone, and lower nuclear reaction rates than the equivalent standard stellar models, and their evolutionary tracks extend to higher temperatures. The magnitude of these effects depends on the amount of injected energy, which, for a fixed set of stellar, rotation and orbital parameters, depends on the perturbed star’s density structure and turbulent viscosity.
Conclusions. Tidal heating offers a possible alternative for describing phenomena such as bloated or overluminous binary components, age discrepancies, and aspherical mass ejection, as well as the extended main sequence turnoff in clusters. However, establishing its actual role requires 3D stellar structure models commensurate with the nonspherically symmetric properties of tidal perturbations.
Key words: binaries: general / stars: rotation / stars: interiors
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1. Introduction
Binary stars are generally assumed to be in equilibrium, thus allowing them to be modeled as single stars during most of their lifetime. In an equilibrium state, the rotation and orbital rates are equal, the orbit is circular and the orbital and rotation axes are parallel. There are, however, many binary stars that are not in equilibrium, such as those in eccentric orbits as well as those in circular orbits that have either not yet attained equilibrium or have departed from this state due to evolutionary changes. In such cases, the tidal interaction excites shearing flows where kinetic energy is converted into heat.
Kopal (1968) addressed the question of whether the energy dissipation rates, Ė, due to such tidal interactions could affect the internal stellar structure, concluding that if only molecular and radiative viscosities are active in the stellar material, then the energy dissipation rates due to shearing motions are too small to have any effect on the star’s structure (see also Kundu 1990). Hence, Ė is generally neglected in the stellar structure calculations. However, Kopal (1968) also noted that if the socalled “turbulent viscosity”, ν_{turb}, is present, then the internal stellar structure could be affected.
Turbulent viscosity is associated with the idea that turbulence in a fluid dynamical system induces an effective viscosity associated with the size of the eddies (Landau & Lifshitz 1987; Shakura & Sunyaev 1988). Numerical and laboratory experiments have been used to estimate its value in stars (Richard & Zahn 1999; Mathis et al. 2004, 2018; Penev et al. 2007). In binary stars with external convective layers, ν_{turb} is assumed to be related to ν_{cv}, which is associated with the convective eddies predicted by the mixinglength theory (Spiegel 1971). The scaling between the tidally induced ν_{turb} and ν_{cv} depends on the largest eddy turnover time and the tidal frequency (Zahn 1966, 1989; Goldreich & Keeley 1977; Goldreich & Nicholson 1977; Duguid et al. 2020). Values of ν_{cv} typically lie in the range between 10^{10}–10^{13} cm^{2} s^{−1} (Penev et al. 2007). However, determining the actual values of ν_{turb}, in addition to determining the conditions under which the shear instability is triggered remain challenging issues (see, e.g., Mathis et al. 2004, 2016; Garaud & Kulenthirarajah 2016, and references therein).
The value of ν_{turb} enters into the expression for calculating Ė and determines the amount of energy released by the motions of the shearing layers. Substantial ongoing efforts are being devoted to determining the impact of ν_{turb} and Ė on the rotation and orbital element evolution (see, e.g., Ogilvie 2014; Terquem 2021; Patel & Penev 2022, and references therein). Much less consideration has been given to evaluating the potential effects of turbulenceinduced heating on the stellar structure, although the dissipative heating rate could exceed the luminosity carried by convection and significantly alter the internal dynamics of stars and planets (Currie & Browning 2017). Also, the role of tidal heating is now well accepted as a mechanism responsible for the Jupiter moon Io’s volcanism (Segatz et al. 1988) and the subsurface liquid oceans in the icy moons of the Solar System (Greenberg et al. 1998; Pappalardo et al. 1999; Hussmann et al. 2002).
Koenigsberger & Moreno (2016, hereinafter KM2016) explored the possibility that tidally induced energy dissipation could cause a star to become more extended than a singlestar counterpart. Focusing on the postmain sequence phases, they found that Ė could act to rapidly expand the outer layers, resulting in a runaway process with each radius increase leading to a larger value of Ė and, therefore, an even larger radius increase. However, these results were based on three approximations. The first is that the tidal perturbation amplitudes were computed only for a surface layer, thus neglecting the contribution from deeper stellar regions and relied on a coarse estimate of the stellar density. The second is the use of a constant ν_{turb} value, while in reality it depends strongly on position and is timevariable. The third simplifying approximation is the assumption that the tidal shear energy dissipation leads only to a radius increase, with no other effect on the stellar structure.
In this paper we perform more realistic calculations for which we have implemented an nlayer calculation of the tidal perturbations, including a prescription for the turbulent viscosity allowing it to be computed as a function of the timevarying and locationdependent velocity fields within the perturbed star. In addition, the structure and evolution of the perturbed star are determined by injecting the tidal energy into a standard structure and evolution model.
The question of the degree to which tidal shear energy dissipation may affect the internal structure of asynchronously rotating binary stars and their evolution has implications for several important phenomena, including the onset of masstransfer and common envelope processes, the subsequent evolutionary path followed by the remnant star and the ejection of circumstellar material. In addition, it could have a bearing on the study of stellar populations if the observable properties of tidally heated stars differ significantly from those of their nonperturbed counterparts.
In Sect. 2, we describe the tidal shear energy dissipation method and the grid of models. In Sect. 3, we describe the stellar structure models that were constructed and the effects of injecting tidal shear energy dissipation into such models. In Sect. 4, we discuss the results and in Sect. 5, we present our conclusions.
2. Tidal shear energy dissipation
2.1. Method
The response of the stellar layers to the external gravitational field of a companion is computed with the nlayer TIDES code^{1} (Koenigsberger et al. 2021), which is based on the numerical method introduced in Moreno & Koenigsberger (1999) and Toledano et al. (2007), and upgraded in Moreno et al. (2011).
The TIDES code can be described as a quasihydrodynamic Lagrangian scheme which simultaneously solves the orbital motion of the companion and the equations of motion of a 3D grid of volume elements covering the inner, rigidly rotating “core” of the tidally perturbed primary star. The core is defined as the interior region that is rotating as a solid body and does not necessarily coincide with the nuclear burning region. The equations of motion include the gravitational acceleration of both stars as well as centrifugal, Coriolis, and gas pressure accelerations. The motions of individual elements are coupled to those of neighboring elements and to the core through viscous stresses. The method is fully described in Moreno et al. (2011) and Koenigsberger et al. (2021).
The initial dimensions of the volume elements are (, , ) determined by the selected 3D grid size. The mass contained in each volume element remains constant over time and is determined by the polytropic stellar structure that is selected. Once the calculation is initiated, the solution of the equations of motion determines the new location of the center of mass of each volume element and its distance from neighboring elements determines the new dimensions (ℓ_{r}, ℓ_{φ}, ℓ_{θ}). These, in turn, are used to compute the new gas pressure within each volume element and the corresponding acceleration term, which is then included in the new acceleration to compute the centers of mass locations in the next integration time step.
The typical depth of the volume elements is ΔR/R_{1} < 0.1, where R_{1} is the primary star radius. The equations are solved in the frame of reference with origin in the center of the primary star and that rotates with the binary orbit, using a seventh order RungeKutta integrator. The secondary is considered to be a pointmass source and its orbital plane is coplanar with the primary star’s equator.
The TIDES calculation captures the effects due to the oscillatory nature of the tidal flows, in addition to those that are due to a differential rotation structure. However, it neglects buoyancy effects and heat and radiation transfer, as well as detailed microphysical processes involving the diffusion and advection of chemical elements which, in most cases, occur on different timescales than the tidal forcing. Hence, the TIDES calculation provides the 3D internal rotation and energy dissipation structures that result from the particular instantaneous dynamical conditions in the system, but the impact of the above processes on the subsequent dynamical evolution cannot be assessed.
The timemarching algorithm is applicable for binary stars with arbitrary rotation velocity and eccentricity, as long as neighboring grid elements retain contact over at least ∼80% of their surface and the centers of mass of two adjoining grid elements do not overlap. Perturbations that depart from these conditions halt the computation. For these very strong perturbations, a full SPH calculation with a different scheme is required, which goes beyond the scope of our current investigation.
The rate of energy dissipation per unit volume is as given in Moreno et al. (2011):
with ν as the kinematical viscosity, ρ as the mass density, ω′ as the angular velocity, and r′, θ′, φ′ are, respectively, the radius, latitude, and longitude coordinates. The primes indicate that the variables are measured in the noninertial frame of reference that rotates with the companion star’s orbital velocity and with its origin at the center of the perturbed star.
For the new version of TIDES used in this paper, we implemented a selfconsistent calculation of turbulent viscosity instead of providing it as a fixed input parameter as in previous versions. Then we chose the simplest formulation (Landau & Lifshitz 1987):
where ℓ_{t} is the characteristic length of the largest eddies that are associated with the turbulence, Δu_{t} is the typical average velocity variation of the flow over the length, ℓ_{t}, and λ is the proportionality parameter that is analogous to α introduced by Shakura & Sunyaev (1988) in the αdisk model.
The largest amplitude perturbation caused by the tidal interaction is in the azimuthal direction (Scharlemann 1981; Harrington et al. 2009). Our nlayer numerical simulation considers only the azimuthal motions and their radial gradients. Thus, we set ℓ_{t} = ℓ_{r} and took Δu_{t} as the typical average velocity variation between a given volume element and the medium surrounding it within this distance, ℓ_{r}. This corresponds to velocity variations between any individual volume element and the elements above and below it.
In the numerical scheme, the value of Δu_{t} is obtained as follows. Let r be the radius of a particular element which we call e_{i} and its linear velocity. Here, is the angular velocity in the frame of reference, S′, that rotates with the binary orbit, and θ is the polar angle. The distance from the center of e_{i} to its top and bottom surfaces is, respectively, r + ℓ_{r}/2 and r − ℓ_{r}/2. The algorithm averages the angular velocity of the volume elements that lie above and below e_{i} and then multiplies this average by (r + ℓ_{r}/2)sinθ (above) and (r − ℓ_{r}/2)sinθ (below). This gives the typical velocities above and below e_{i}, which we call v_{top} and v_{bottom}, respectively. Hence, the average velocity of the elements surrounding e_{i} is v_{ave} = (v_{top} + v_{bottom})/2 and the typical average velocity variation of the flow over the length scale ℓ_{r} is Δu_{t} = v − v_{ave}.
The values of ℓ_{r} and Δu_{t} obtained as described above for each volume element in the grid are inserted in Eq. (2). The total viscosity is computed as ν = ν_{turb} + ν_{molec}, where ν_{molec} is the molecular viscosity and it is given as input. This value of ν at each volume element is now used by the TIDES algorithm to solve the equations of motion in the next time step of the calculation.
2.2. TIDES input parameters and model runs
The modeled binary star is called the “primary”. Its mass is m_{1} and it has an initial, unperturbed radius, R_{1}. Its initial rotation condition is one of rigid rotation with its rotation velocity given in terms of the synchronicity parameter, β_{0} = ω_{0}/Ω_{0}, where ω_{0} is the angular velocity of the inner core (assumed throughout to rotate as a rigid body) and Ω_{0} is the orbital angular velocity at a reference point in the orbit, for example, periastron in the case of nonzero eccentricity. The TIDES input parameters are described in Table 1, where we also list the values that were kept constant. The stellar masses and orbital period were chosen to be the same as in KM2016.
Description of TIDES input parameters.
In this paper, we mainly consider rotation rates that are close to corotation. Specifically, we analyzed cases with β_{0} = 0.95, 1.05 and 1.10. We also probed the effects for stellar radii in the range R_{1} = 0.97 − 2.25 R_{⊙}, which correspond to those of a 1 M_{⊙} star from the time it reaches the main sequence (MS) until shortly after the terminal age main sequence (TAMS). These radii are smaller than the Roche Lobe radius of the primary star (R_{RL} = 2.5 R_{⊙}). The radii that were probed and the corresponding surface equatorial velocities V_{rot} in the unperturbed star are listed in Cols. 2 and 6, respectively, of Tables 2 and 3.
TIDES fivelayer models.
TIDES tenlayer models.
The tidal shear energy dissipation rates depend on the tidal velocity gradient, the viscosity, and the density, with the latter depending on the assumed polytropic structure. The outer convective layers of lowmass stars are generally represented with a n = 1.5 polytrope. However, the inner layers correspond to polytropes with larger n values, depending on the evolutionary stage. For the early main sequence models, we chose n = 3, which provides a closer match to the values of the density in the outer layers. For later evolutionary states, we considered polytropic indices that provide a best approximation to the density structure of the MESA model in the outer layers.
The numerical experiments to determine the tidal shear energy dissipation rates are listed in Table 2 and are divided into three blocks. The aim of each block is to examine the dependence on particular input parameters.
In the first block, we probe the dependence on the stellar radius and λ (the turbulent viscosity coefficient; see Eq. (2)) while holding β_{0} = 0.95. Calculations were performed for a range in stellar radii R_{1}/R_{⊙} = [0.99, 2.25]. For each radius, three viscosity representations were tested: a constant value as described in Koenigsberger et al. (2021), with the same values that were used in KM2016 and listed here in Col. 4 of Table 2; values computed with λ = 0.1; and values computed with λ = 1. The value of λ is listed in Col. 5 of Table 2. The second block of numerical experiments listed in Table 2 was performed with β_{0} = 1.05. This means that the primary star’s rotation is slightly supersynchronous. Here, we fixed λ = 1 and tested different polytropic indices as well as different stellar radii within the range tested in block 1.
The third block of numerical computations was performed with a slightly greater departure from synchronicity than those in block 2, namely β_{0} = 1.10. All other input parameters were kept the same. Table 3 shows results of numerical computations with ten layers, instead of the five that are used in the models listed in Table 2. These models probe the behavior of layers that lie closer to the nuclear burning core. Since we retained the same layer thickness for all computations, adding more layers allows us to probe regions that lie deeper in the star. For example, the models in Table 2 include layers down to 0.7 R_{1}, where R_{1} is the unperturbed equilibrium radius of the star, while the corresponding model in Table 3 reaches down to 0.4 R_{1}. We have also performed experiments with 20 layers, arriving at nearly the same results as with ten layers because the perturbed velocities decline very rapidly at smaller radii. Increasing the grid size beyond ten layers significantly increases the processing time and is deemed unnecessary for our current purposes. For problems where the main concern is modeling the layers near the surface, we find that a reduced number of layers yields results that are comparable to a larger radial grid computation. This is illustrated in Table B.1, where we list energy dissipation rates in each layer obtained from a fivelayer and a tenlayer computation.
The TIDES computation provides the angular velocity, viscosity, and energy dissipation rates as a function of azimuth angle φ for each radius, r, and colatitude, θ, of the computational grid and as a function of time, t. Thus, the energy dissipation rates are represented as Ė_{r,θ,φ,t}. The temporal dependence cannot be neglected because asynchronous binaries undergo orbitalphase dependent variability. Therefore, for each model run, we chose to output the data at 40 orbital phases within each of five orbital cycles. An inspection of the five cycles allows for an assessment of longterm variability patterns and also whether the transitory state of the numerical integration has transpired. The latter has usually occurred within ∼200 orbital cycles. As the nature of the phasedependent variability is not the subject of this paper, Ė_{r,θ,φ,t} was averaged over the 40 orbital phases within an orbital cycle in which the calculation has reached the stationary state. This yields Ė_{r,θ,φ}. The radial profile Ė_{r} is obtained by integrating over θ and φ.
For a general characterization of the models in Tables 2 and 3, we opted to list the energy dissipation rates in the deepest layer of the calculation and also in the layer that neighbors the surface. The deepest layer is represented by Ė_{k = 1}, to indicate that it refers to the first layer, which in Table 2 corresponds to r ∼ 0.7 R_{1} and in Table 3, it corresponds to r ∼ 0.4 R_{1}. The layer that is contiguous to the surface layer is represented by Ė_{k = 4} in Table 2 and Ė_{k = 9} in Table 3, both of which correspond to r ∼ 0.9 R_{1}. We note that the k number index denotes the layer, with k = 1 corresponding always to the layer that lies closest to the core.
2.3. Results
An example of the results is illustrated in Fig. 1, where we plot the data obtained in model 63. The top panel illustrates the angular velocity, ω″, in the equatorial latitude (θ = 90°), measured in the rotating frame of the primary star^{2} as a function of the azimuth angle. The sinusoidal shape corresponds to the equilibrium tide. The peaktopeak amplitude of each curve decreases with decreasing distance from the stellar center, as expected from the tidal force. These characteristics are shared by all the models that were run for this paper, as shown in Appendix A. Another general feature is that the overall shape and the peaktopeak amplitude do not depend significantly on the viscosity for the range of viscosity values that were explored, a result that was already found in the onelayer calculations that were performed by KM2016.
Fig. 1. Angular velocity, viscosity and dissipated energy as a function of azimuth angle in the ten layers of model 63. The azimuth angle is measured in the direction of the companion’s orbital motion and φ = 0 corresponds to the subbinary longitude. Top: angular velocity ω″ at the equator in the perturbed star rest frame. The surface layers (largest amplitude curves) are more strongly coupled to the tidal field than the inner layers and thus are forced to lag behind the supersynchronously rotating core. Middle: turbulent viscosity values computed by the model. Bottom: tidal shear energy dissipation rate. 
The corresponding behavior of the viscosity in model 63 is illustrated in the middle panel of Fig. 1. Its value in the layers closest to the surface is a few orders of magnitude larger than near the core. Thus, the inner layers are significantly less coupled to the outer layers compared to the case in which a constant viscosity is used (models 2, 7, 11, and 23). Thus, the outer layers approach synchronous rotation more rapidly than the inner layers. This results in a radial gradient in the average velocity of each layer, a result that is consistent with the conclusion of Goldreich & Nicholson (1989), asserting that a star synchronizes from the surface inward.
The differential rotation structure that is obtained with a variable viscosity is in contrast with the uniform average rotation structure in models 2, 7, 11, and 23, where the viscosity is kept constant, illustrating the important role that is played by this parameter.
We list in Col. 7 of Table 2 the maximum viscosity, ν_{max}. It appears in the calculation always near the surface and around the equatorial latitude. For stars with radii between 0.99 R_{⊙} and 2.25 R_{⊙}, we find ν_{max} = 4 × 10^{13} cm^{2} s^{−1}–2 × 10^{15} cm^{2} s^{−1} (λ = 1).
The bottom panel of Fig. 1 shows Ė_{r,φ}, the energy dissipation rate for each radius as a function of azimuth. Its behavior mimics the φdependence of the viscosity as might be expected given that the density remains relatively constant within each layer and only the velocity gradients change as a function of azimuth angle.
We give the energy dissipation rates in the deepest layer, Ė_{k = 1}, and the layer immediately below the surface layer, Ė_{k = 4}, in Cols. 8 and 9, respectively, of Table 2. The last column of this table lists Ė_{tot}, the total energy dissipation rate obtained by integrating Ė_{r} over all layers. The relatively small difference between Ė_{tot} and Ė_{k = 4} is due to the fact that the layers near the surface are responsible for the largest share of the dissipated energy. In fact, the greatest contribution to the total energy dissipation rate for models with a polytropic index > 1.5 arises in the layer that is adjacent to the surface (in this case, Ė_{k = 4}), not in the surface layer. A priori, it would appear curious that the energy dissipation rate in the surface layer is not the greatest since it is the one subject to the strongest tidal perturbations. There are two explanations for this apparent contradiction. The first resides in the fact that only the inner boundary of the surface interacts with a neighbor, whereas all other layers interact with both a layer above and one below. The second factor is based on the interplay between the three terms that enter into the energy dissipation calculation (see Eq. (1)). When the polytropic index > 1.5, the density decrease near the surface is more pronounced than the increase in viscosity and in the velocity gradients. This effects is illustrated by comparing models 65−67, where Ė_{tot} decreases by over two orders of magnitudes due to the decreasing surface density that results from changing the polytropic index from 1.5 to 3.8.
The depth dependence of Ė_{r} for different density structures is most clearly illustrated with models that are computed with layers that are deeper, such as those listed in Table 3. These are tenlayer runs for which Ė_{k = 1} correspond to a layer that lies at ∼0.4 R_{⊙}. A layerbylayer comparison down to this depth clearly illustrates the densitydependence of Ė_{r}, as shown in Table 4 by comparing models 62 and 63.
Energy dissipation rate radial profiles from the TIDES tenlayer models.
The energy dissipation rate is very sensitive to the value of the synchronicity parameter, β_{0}. For example, models 66 and 73 have identical input parameters except for β_{0}: the latter model is slightly more supersynchronous than the former (β_{0} = 1.10 versus β_{0} = 1.05). The maximum viscosity value of model 73 is approximately twice as large as that of model 66, and Ė_{tot} is one order of magnitude larger. Model 71 illustrates the case in which the inner core rotates significantly slower than the other cases studied (β_{0} = 0.20), but because of the large departure from synchronicity, it has one of the largest values of Ė in all layers. Thus, increasing departures from synchronicity while holding other parameters constant does lead to increasing energy dissipation rates.
Finally, the nearzerovalues displayed in the viscosity plots merit a comment. The value of ν_{turb} is calculated using the instantaneous velocity gradients. There are longitudes in our calculations at which these gradients vanish and thus the minimum ν_{turb} → 0 (see Fig. 1). However, in reality, the transfer of the kinetic energy of the flow into eddies that then act as the viscosity source is not instantaneous, nor are the eddies expected to disappear instantaneously when the velocity gradient vanishes. Thus, there may be a minimum viscosity that is larger than the molecular viscosity even when the velocity gradients momentarily vanish. Furthermore, in convective envelopes, a baselevel viscosity associated with the convective eddies are expected to always be present.
3. Stellar structure models
The stellar structure models were computed with version 15140 of the MESA opensource 1D stellar evolution code^{3} (Paxton et al. 2011, 2013, 2015, 2018, 2019). We ran models for a 1 M_{⊙} star with solar abundances (Z = 0.02) and an initial surface rotation velocity of 40 km s^{−1} on the zero age main sequence (ZAMS). All the models were run until the stellar radius reached 2.25 R_{⊙}, which occurs shortly after the central H fraction falls below 10^{−4}. By this point, the surface rotation velocity has fallen to about 10 km s^{−1}.
We include the effect of energy dissipation on the stellar envelope using the other_energy hook provided in the MESA code. The extra heat we inject into the envelope is taken from the energy dissipation rate as a function of radius, Ė_{r}, of a TIDES binary model. In particular, we chose the energy dissipation profile of model 63, which is listed in Col. 4 of Table 4. The parameters of this model are given in Table 3. The TIDES code is not an evolution code and the radius of a given model is fixed, whereas the MESA simulations show that the radius of the star changes during the main sequence. We account for this in our models by normalizing the radius of the TIDES energy dissipation profile and apply the profile to the normalized stellar radii of the MESA models. In addition, MESA has many more interior radial points than the TIDES model and we find the energy dissipation rate at each of these by linear interpolation between the nearest TIDES model points.
We run three sets of models summarized in Table 5: (i) a standard model (h0) which evolves the star from the ZAMS to the stopping point with no added heat; (ii) set 2, in which heat is injected steadily into the stellar envelope with the TIDES energy dissipation rate starting at the ZAMS; and (iii) set 3, in which we emulate the change from synchronous rotation during the MS to asynchronous rotation as the primary star approaches the terminal age main sequence (TAMS; defined in this paper as the time when the central hydrogen mass fraction falls below 10^{−4}). In set 3, we multiply all of our dissipation rates by the timedependent factor of
Heat injection in MESA models.
which increases sharply to unity as t → t_{TAMS}, where it becomes saturated. Table 6 lists properties of the three sets.
Ages and events in MESA models.
It is evident from Table 4 that the energy dissipation rate has a strong dependence on radius, internal structure, and rotation rate. The aim of this paper is to simply explore the effect on the stellar structure of this type of energy input. Thus, in order to take into account larger energy input rates than given by the TIDES model 63, we also computed models with the model 63 profile multiplied by a factor of 10. These models are labeled set2h10 and set3h10 for the set 2 and set 3 cases, respectively. Such higher heating rates are obtained, for example, by increasing the β_{0} value from 1.05 to 1.10, while holding the polytropic index constant or by increasing the radius from 0.99 to 1.4 R_{⊙}; this can be seen by comparing models 66 and 73 and models 4 and 12, respectively, in Table 2.
3.1. Results for constant injected heat
The effect of adding heat into the stellar envelope on basic stellar physical properties is illustrated in Fig. 2. In the top panels, we plot the stellar radius and luminosity as a function of time for the standard model h0, and for models set2h1 and set2h10, which correspond to timeconstant heat injection throughout the MS, beginning at the ZAMS. For the same age, the h1 and h10 models have larger radii and luminosities, compared to the standard model, h0.
Fig. 2. Properties of the MESA models of set 2, in which the heating profile is introduced continuously throughout the evolutionary trajectory. Top: stellar radius (left) and luminosity (right). The blue curve corresponds to the nonheated models. The orange and green curves correspond, respectively, to the heating profile of model 63 and a heating profile that is ten times larger. Bottom: protonproton reaction rate power (left) and CNO reaction rate power (right). The blue curve corresponds to the nonheated models and its ordinate is on the right side of the plot and given in units of log(L_{⊙}). The orange and green curves correspond to the ratio of rates in the heated and nonheated models, with colors the same as in the top panels. 
The bottom panels show the total power generated by the protonproton chain and through the CNO cycle as a function of time. Throughout the MS, the pp chain dominates energy production but it is clear that the extra luminosity in the outer layers changes the energy transport in the envelope and has repercussions on the energy production in the stellar core, since the nuclear energy generation rate in the heated models is slightly smaller than in the standard model.
We can also examine the internal structure of the model stars at any point in their evolution. The internal luminosity structure for ages of 2 Gyr and ages between 9.257 and 11.644 Gyr is shown in Fig. 3. There is a significantly different structure in the heated models compared to the standard model. The luminosity inside ∼80% of the maximum radius is smaller in the heated models than in the standard model, while the outer ∼20% has a higher luminosity. Thus, the surface luminosities of the heated stars are higher than those of the standard model.
Fig. 3. Luminosity structure of MESA set 2 models at ages of (top left) 2 Gyr and ages between (top right) 9.257, (bottom left) 9.653 and (bottom right) 11.644 Gyr. The blue curves correspond to the models without added heating; orange to those with heating as given by model 63; green to heating by ten times that given by model 63. The tidally heated models have surface luminosities that are larger than the nonheated models, but with lower luminosities in deeper layers. 
In Fig. 4 we illustrate the Kippenhahn diagrams for the h0 and the h10 models^{4}. These enable us to study how energy transport in different regions of the stellar interior changes as a function of time. The most notable difference is the reduced size of the convective region (green shaded region) in the h10 models compared to the unperturbed models, both on the MS and afterward. A similar, but much less prominent reduction is present in the h1 models. This is because adding heat to the stellar envelope flattens the temperature gradient and energy can be transported by radiative diffusion almost up to the surface.
Fig. 4. Kippenhahn diagram for the standard model (top) and the set2h10 model (bottom). The panels from left to right represent: Left: evolution from ZAMS to TAMS (defined as abundance of H first drops below 10^{−4} in the stellar core). Center: evolution from the TAMS to 90% of the remaining time. Right: final 10% of the postTAMS evolution. The ordinate is the mass coordinate that runs from the center of the star to the surface. The nuclear energy generation rates are indicated by the different shades of blue. The green hatches indicate zones in which the energy transport is dominated by convection. 
3.2. Results for gradual injected heat
An overview of the results for gradual timeincreasing heat injection is illustrated in Fig. 5. In contrast to the set 2 models shown in Fig. 2, there is no significant difference throughout the MS between the heated and the unperturbed models; this is as expected since the injection of heat is very small until the TAMS is approached. However, at the TAMS, the set 3 models very rapidly approach the properties of the set 2 models at similar ages.
Fig. 5. Properties of the MESA models of set 3, in which the heating profile is introduced gradually, starting at the ZAMS and reaching maximum at the TAMS. Top: stellar radius (left) and luminosity (right). The blue curve corresponds to the nonheated models. The orange and green curves correspond, respectively, to the heating profile of model 63 and a heating profile that is ten times larger. Bottom: protonproton reaction rate power (left) and CNO reaction rate power (right). The blue curve corresponds to the nonheated models and its ordinate is on the right side of the plot and given in units of log(L_{⊙}). The orange and green curves correspond to the ratio of rates in the heated and nonheated models, with colors the same as in the top panels. 
The internal luminosity structure of set 3 models at ages 9.257, 9.653, and 11.644 Gyr is illustrated in Fig. 6, which shows similar characteristics as those of set 2. Specifically, the heated models have a lower internal luminosity than the standard model and a significantly larger surface luminosity. The similarity between the set 2 and set 3 heated models is best exemplified with a HertzsprungRussell Diagram (HRD), as shown in Fig. 7, where we see that the set2h10 and the set3h10 tracks are identical after the MS turnoff, which is when both sets are being equally heated.
Fig. 6. Luminosity structure of MESA set 3 models at ages 9.257, 9.653, and 11.644 Gyr. Blue curves correspond to the models without added heating; red to those with heating as given by model 63; green by heating ten times that given by model 63. The tidally heated models have a surface luminosity that is larger than the nonheated models, but lower luminosity in deeper layers. 
Fig. 7. Evolutionary tracks on the HRD of the standard model (light blue), the set3h1 model (orange), and the set3h10 model (green) showing the result of introducing a gradual heating near the end of the main sequence. The dark blue curve corresponds to the set2h10 model in which there was constant heating throughout the main sequence. 
4. Discussion
In this paper, we explore the manner in which extra energy that is injected into the external stellar layers affects its structure and its observational properties. The source of this extra energy is here assumed to be that which is released by the action of shearing layers whose motions are driven by the tidal interaction with a close companion. We adopted a simple prescription for the turbulent viscosity, ν_{turb} = λℓ_{t}Δu_{t}, where Δu is the instantaneous velocity difference between two contiguous layers, ℓ_{t} is a characteristic distance taken to be the separation of the layers, and λ is a proportionality factor. The instantaneous values of Δu are computed using the TIDES code, which solves the equations of motion in the rotating reference frame of the binary, taking into account the gravitational, centrifugal, Coriolis, gas pressure, and viscous accelerations. The stellar structure and evolution is computed with MESA.
The nominal calculation performed to probe the effects of tidal shear energy dissipation on the stellar structure is for a 1 M_{⊙} star with an initial radius of 0.97 R_{⊙} in a 1.44 d orbit and a 0.8 M_{⊙} companion. The TIDES computational grid consists of ∼14 300 volume elements (one hemisphere, northsouth symmetry is assumed) covering the outer ∼54% in radius of the star. The synchronicity parameter β_{0} = 1.05 is chosen such that the core of the star rotates at a rate that is 5% faster than corotation. Additional models are computed for stellar radii up to 2.25 R_{⊙}, polytropic indices 1.5 ≤ n ≤ 3.8, and synchronicity parameters 0.20 ≤ β_{0} ≤ 1.1.
4.1. Viscosity and tidal shear energy dissipation
For the ∼1 R_{⊙} models, we find maximum turbulent viscosity values near the stellar surface at the equator in the range 5×10^{−4} to 10^{−3}/d in calculations with λ = 1. These values correspond to 3−6×10^{13} cm^{2} s^{−1}, respectively. The corresponding total energy dissipation rates integrated over the entire star Ė_{tot} are in the range 2−20×10^{30} erg s^{−1}. These values decrease in proportion to the value of λ.
As a star evolves and becomes larger, the turbulent viscosity grows and, assuming that the density structure does not significantly change, the energy dissipation rate also increases. The maximum value of ν_{turb} obtained with λ = 1 in this paper is not very different from the values used in Koenigsberger & Moreno (2016, henceforth KM2016), except for the largest radii listed in block 1 of Table 2, where we see a difference of a factor of ∼3. The difference arises because the value of ν in the KM2016 is an input parameter which, broadly speaking, is unconstrained. To avoid using a completely arbitrary ν value, it was estimated assuming it to scale with the characteristic spatial and velocity dimensions. This was not a very precise estimate. In our current formulation, ν is computed internally, allowing for a more consistent estimate.
Our current Ė_{tot} values are significantly smaller compared to those listed in KM2016 (their Table 3). For example, they list Ė_{n = 1.5} ∼ 10^{34} ergs s^{−1} for the R_{1} = 0.99 R_{⊙} model, while we obtain Ė = 2 × 10^{30} ergs s^{−1} for our model 4. The dominant source of this difference is that KM2016 assumes a n = 1.5 polytropic index, while all the models in block 1 of our Table 2 were computed with n = 3. The density at a radius of 0.965 R_{⊙} is ∼3 orders of magnitude higher for n = 1.5 than for n = 3. We performed an analogous onelayer computation to that of KM2016 using their same model, but this time with n = 3 and obtained Ė_{n = 3} = 6 × 10^{30} ergs s^{−1}, which is only a factor ∼3 larger than what we obtain in the current calculations for model 4. This remaining difference is due to a combination of factors. The first is that in the KM2016 model, the viscosity is constant over the entire surface layer, while in our current model its value significantly decreases at various locations along the azimuthal coordinate and, especially, in the polar direction. The second is that the KM2016 model includes radial motions which, although they are approximately ten times smaller than the azimuthal motions, also contribute toward the energy dissipation rates.
The lower energy dissipation rates that we find in this paper could impact the interpretation of the V1309 Sco merger phenomenon in terms of a tidallyinduced runaway process, as discussed in KM2016. However, the higher energy dissipation rates needed for a rapid orbital evolution timescale can still be obtained in our current model by increasing the departure from synchronicity. While it is beyond the scope of this paper, a reanalysis using our current model, but relaxing the conditions on β_{0} and the polytropic index, is warranted before abandoning the tidal runaway scenario.
4.2. Implications for stellar structure
The energy dissipation rates computed by TIDES for each layer were injected into MESA stellar structure and evolution calculations. With even a small (∼5%) departure from synchronicity, the radius and luminosity values of the tidally heated stars are larger than those of the equivalent unperturbed model at all evolutionary times. The star also has a smaller surface convective region and lower nuclear processing rates, the latter allowing the tidally perturbed star to live longer. The differences between an asynchronous binary and its unperturbed counterpart depend on the amount of injected energy which, for a fixed set of stellar and orbital parameters, depends on how much the stellar rotation departs from synchronicity and the value of turbulent viscosity.
From an observational perspective, determining whether a star is truly in synchronous rotation is a challenging problem as the only available information is the projected surface equatorial speed. Because the synchronization time scales as and the turbulent viscosity is largest near the surface, the star tends to synchronize from the surface inward (Goldreich & Nicholson 1989; Koenigsberger et al. 2021). Thus, stars in circular orbits that are thought to be synchronized may actually retain an internal angular velocity gradient upon which tidally excited oscillations are superposed. Furthermore, all eccentric binary systems are asynchronously rotating during most of their orbital trajectory; hence, they suffer from tidal perturbations regardless of their age. This would be particularly true of recentlyformed binaries in dense regions of stellar clusters where close encounters of single stars are believed to frequently occur.
In this context, it is interesting to note that radii of lowmass binary stars in shortperiod orbits have been determined to be as much as 10% larger than their counterparts in longperiod orbits (Hoxie 1973; Lacy 1977; Popper 1997; Clausen et al. 1999; Torres et al. 2006), consistent with the radius increase that we find in the tidally heated stars explored in this paper. Many lowmass stars are associated with significant surface activity, as evidenced by their light curves and the emission cores of Ca II H and K lines (Torres et al. 2006; LópezMorales 2007; Morales et al. 2008). This activity is explained in terms of the presence of strong surface magnetic fields and these fields are thought to inhibit efficient convection (Torres et al. 2006; Chabrier et al. 2007; Clausen et al. 2009; Feiden 2016). Tidal flows and magnetic fields are not mutually exclusive, but the manner in which these two physical processes might interact is an open question.
Our scenario for tidal shear energy dissipation may also have a bearing on the mass discrepancy problem in massive stars that has been known for several decades, but that has still eluded explanation. The discrepancy, first noted by Herrero et al. (1992), consists of the fact that the masses derived from spectroscopic analysis are systematically lower than those found from evolutionary models; alternatively, they are more luminous than predicted by the evolutionary models. In their analysis of a set of eclipsing binary stars, Massey et al. (2012) found them to be on average 11% less massive or conversely, 0.2 dex more luminous, as compared to stellar structure models. Because it is now generally accepted that a large majority of massive stars are in binary systems, it is tempting to suggest that a possible solution to the massdiscrepancy problem may reside in the phenomena we discuss in this paper.
4.3. Implications for stellar evolution and population synthesis
The effects of binary interactions on stellar evolution have, until now, focused mainly on the effects of massloss and mass transfer between the components in late evolutionary stages. In the case of massive stars, binary interactions during the main sequence have been incorporated only indirectly in the sense that tides are invoked to maintain shortperiod binary stars in rapid rotation, allowing them to be treated as rapid rotators throughout their main sequence lifetime. The possible modification in the stellar structure that results from the tidal interactions, however, is generally neglected. These effects include tidal shear energy dissipation and turbulent viscosity, which, in turn, have an impact on the internal energy budget as well as the rates of angular momentum and chemical transport.
The presence of tidal perturbations does not depend on the age of a star, but it may be most easily detected in very close binaries or those in which the radius of one of the stars is significant compared to the orbital separation – specifically, stars at the end of the main sequence. We find that evolutionary tracks of our tidally heated stars extend further to the blue during the end stages of the main sequence than does the track for the standard model. This effect is reminiscent of the extended main sequence turnoff (eMSTO) phenomenon that was first detected by Bertelli et al. (2003) and Mackey & Broby Nielsen (2007) in the Large Magellanic Cloud clusters NGC 2173 and NGC 1846; this is now considered an ubiquitous feature of Magellanic Cloud Clusters with ages between ∼20 Myr and ∼2 Gyr (Milone et al. 2018; Goudfrooij et al. 2014). The eMSTO has been interpreted as the result of a prolonged star formation (Mackey et al. 2008; Glatt et al. 2008; Goudfrooij et al. 2011; Keller et al. 2011) or as a result of stellar rotation (Bastian & de Mink 2009). We suggest that tidal heating may be an additional potential explanation. In globular clusters, the presence of blue straggler stars not showing evidence of masstransfer (Ferraro et al. 2006) could also potentially be a manifestation of the effectiveness of this process.
Another interesting application of our model refers to the stability of a binary star’s outer layers as it leaves the main sequence and heads up the giant branch. As the stellar radius increases, so do the tidal amplitudes (Appendix A). We can speculate that the growing surface velocities could attain the escape velocity before the Roche Lobe radius is reached. Because the tidal amplitudes are largest around the equator, this is where the star would become most bloated and where mass loss might be expected to occur first and take the form of an excretion disk. Observational evidence for asymmetrical massloss episodes is found in planetary nebulae, many of which display bipolar morphologies (Soker & Rappaport 2001; De Marco 2009; Jones 2011; Jones & Boffin 2017; Boffin & Jones 2019), as well as the presence of structures such as jets, rings, and halos (Harpaz et al. 1997; Kwok 2002; Phillips et al. 2009). Therefore, taking into account the perturbations caused by tidal forces on the progenitor star during its expansion stages may contribute to improving our understanding of the processes that give rise to the wide variety of morphologies of such interstellar structures. Finally, we note that although our focus is the potential effects on the stellar structure and evolution due to tidal perturbations, any nonstandard process that can heat subsurface stellar layers would produce the same effects we describe in this paper.
4.4. Caveats
There are several important caveats to our results. The first is that our simplified prescription for the turbulent viscosity depends on a λ parameter, which can be associated with the fraction of kinetic energy in the shearing flows that is transformed into turbulent eddies. Our calculations were performed for λ = 0.1 and 1, and we find that the values of Ė scale approximately linearly with λ. Significant differences between the heated and the standard MESA models appear mainly for the larger λ value, which is likely to be unrealistic. However, Ė values large enough to produce significant differences can also be obtained by increasing the synchronicity parameter, as we showed for the cases with β_{0} = 1.05 and 1.1, or decreasing it as illustrated for β_{0} = 0.2. This means that even if λ is small, stars having significant departures from synchronous rotation could be observed to display the tidal heating effects that we have described.
An underlying assumption in our treatment is that the criteria for triggering shear instabilities are met, an issue that depends on the hydrodynamical properties and the microphysics of the fluid (Garaud & Kulenthirarajah 2016, and references therein) processes that our approach cannot compute. In the case of solartype stars, turbulent viscosity is associated with the convective eddies, so in principle, the criteria for triggering the shear instability are met. However, the nature of the interaction between convective motions and the tidal flows remains to be resolved (Terquem 2021; Vidal & Barker 2020; Duguid et al. 2020).
Finally, it is important to keep in mind that Ė is highly nonisotropic in 3D space, attaining maximum values near the equator and decreasing toward the poles, as well as having a dependence on azimuth. Although the TIDES model captures the 3D tidal perturbation structure, the MESA models are currently only 1D. The Ė radial profile that was injected into the h1 and h10 MESA models corresponds to the total energy dissipation in each layer and, hence, the heating is strongly concentrated around the equator. This would lead to equatorial bloating unless the horizontal energy transport processes are sufficiently efficient to redistribute the added energy.
5. Conclusions
Stars in a binary system can interact in different ways, depending on both their evolutionary stage and their orbital parameters (Soker 1997), and these interactions have an impact on the stellar evolution of the components. In this paper, we explore the potential role of tidal shear energy dissipation in altering not only the evolutionary path of a star in its postmain sequence stages, but also its internal structure during the main sequence. Tidal heating offers a possible alternative for describing discrepancies between observations and the standard stellar structure models. Examples include phenomena such as the eMSTO in clusters, bloated or overluminous binary components, age discrepancies, and aspherical mass ejection. However, establishing the actual role of tidal heating requires incorporating the nonspherically symmetric properties of the tidal perturbations into stellar structure models. It also requires a hydrodynamical approach to determining the turbulent viscosity of the stellar fluid.
The doubleprime notation indicates that it is measured with respect to the reference frame S″, which rotates at the same constant rate as the primary star core, as opposed to the S′ reference frame that rotates with the orbital motion of the companion which for eccentric orbits is not constant (see Koenigsberger et al. 2021).
These diagrams were made with the mkipp.py software written by Pablo Marchant https://github.com/orlox/mkipp/
Acknowledgments
We acknowledge support from CONACYT project 252499 and DGAPA/PAPIIT projects IN103619 and IN105723. We thank an anonymous referee for the very insightful and helpful comments. G.K. thanks Catherine Pilachowski and Constantine Deliyannis for enlightening discussions.
References
 Bastian, N., & de Mink, S. E. 2009, MNRAS, 398, L11 [NASA ADS] [CrossRef] [Google Scholar]
 Bertelli, G., Nasi, E., Girardi, L., et al. 2003, AJ, 125, 770 [NASA ADS] [CrossRef] [Google Scholar]
 Boffin, H. M. J., & Jones, D. 2019, The Importance of Binaries in the Formation and Evolution of Planetary Nebulae (Switzerland: Springer Nature) [Google Scholar]
 Chabrier, G., Gallardo, J., & Baraffe, I. 2007, A&A, 472, L17 [CrossRef] [EDP Sciences] [Google Scholar]
 Clausen, J. V., Baraffe, I., Claret, A., & Vandenberg, D. A. 1999, in Stellar Structure: Theory and Test of Connective Energy Transport, eds. A. Gimenez, E. F. Guinan, & B. Montesinos, ASP Conf. Ser., 173, 265 [NASA ADS] [Google Scholar]
 Clausen, J. V., Bruntt, H., Claret, A., et al. 2009, A&A, 502, 253 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Currie, L. K., & Browning, M. K. 2017, ApJ, 845, L17 [NASA ADS] [CrossRef] [Google Scholar]
 De Marco, O. 2009, PASP, 121, 316 [NASA ADS] [CrossRef] [Google Scholar]
 Duguid, C. D., Barker, A. J., & Jones, C. A. 2020, MNRAS, 497, 3400 [Google Scholar]
 Feiden, G. A. 2016, A&A, 593, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ferraro, F. R., Sabbi, E., Gratton, R., et al. 2006, ApJ, 647, L53 [NASA ADS] [CrossRef] [Google Scholar]
 Garaud, P., & Kulenthirarajah, L. 2016, ApJ, 821, 49 [Google Scholar]
 Glatt, K., Grebel, E. K., Sabbi, E., et al. 2008, AJ, 136, 1703 [Google Scholar]
 Goldreich, P., & Keeley, D. A. 1977, ApJ, 212, 243 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Nicholson, P. D. 1977, Icarus, 30, 301 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Nicholson, P. D. 1989, ApJ, 342, 1079 [Google Scholar]
 Goudfrooij, P., Puzia, T. H., Chandar, R., & KozhurinaPlatais, V. 2011, ApJ, 737, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Goudfrooij, P., Girardi, L., KozhurinaPlatais, V., et al. 2014, ApJ, 797, 35 [CrossRef] [Google Scholar]
 Greenberg, R., Geissler, P., Hoppa, G., et al. 1998, Icarus, 135, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Harpaz, A., Rappaport, S., & Soker, N. 1997, ApJ, 487, 809 [NASA ADS] [CrossRef] [Google Scholar]
 Harrington, D., Koenigsberger, G., Moreno, E., & Kuhn, J. 2009, ApJ, 704, 813 [NASA ADS] [CrossRef] [Google Scholar]
 Herrero, A., Kudritzki, R. P., Vilchez, J. M., et al. 1992, in The Atmospheres of EarlyType Stars, eds. U. Heber, & C. S. Jeffery, 401, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Hoxie, D. T. 1973, A&A, 26, 437 [NASA ADS] [Google Scholar]
 Hussmann, H., Spohn, T., & Wieczerkowski, K. 2002, Icarus, 156, 143 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, D. 2011, PhD Thesis, University of Manchester, UK [Google Scholar]
 Jones, D., & Boffin, H. M. J. 2017, Nat. Astron., 1, 0117 [Google Scholar]
 Keller, S. C., Mackey, A. D., & Da Costa, G. S. 2011, ApJ, 731, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Koenigsberger, G., & Moreno, E. 2016, Rev. Mex. A&A, 52, 113 [Google Scholar]
 Koenigsberger, G., Moreno, E., & Langer, N. 2021, A&A, 653, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kopal, Z. 1968, Ap&SS, 1, 411 [NASA ADS] [CrossRef] [Google Scholar]
 Kundu, P. K. 1990, Fluid Mechanics (San Diego: Academic Press) [Google Scholar]
 Kwok, S. 2002, Am. Astron. Soc. Meet. Abstr., 201, 22.06 [NASA ADS] [Google Scholar]
 Lacy, C. H. 1977, ApJS, 34, 479 [Google Scholar]
 Landau, L. D., & Lifshitz, E. M. 1987, Fluid Mechanics (Oxford: Pergamon Press) [Google Scholar]
 LópezMorales, M. 2007, ApJ, 660, 732 [Google Scholar]
 Mackey, A. D., & Broby Nielsen, P. 2007, MNRAS, 379, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Mackey, A. D., Broby Nielsen, P., Ferguson, A. M. N., & Richardson, J. C. 2008, ApJ, 681, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Massey, P., Morrell, N. I., Neugent, K. F., et al. 2012, ApJ, 748, 96 [NASA ADS] [CrossRef] [Google Scholar]
 Mathis, S., Palacios, A., & Zahn, J. P. 2004, A&A, 425, 243 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., AuclairDesrotour, P., Guenel, M., Gallet, F., & Le PoncinLafitte, C. 2016, A&A, 592, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., Prat, V., Amard, L., et al. 2018, A&A, 620, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Milone, A. P., Marino, A. F., Di Criscienzo, M., et al. 2018, MNRAS, 477, 2640 [Google Scholar]
 Morales, J. C., Ribas, I., & Jordi, C. 2008, A&A, 478, 507 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Moreno, E., & Koenigsberger, G. 1999, Rev. Mex. A&A, 35, 157 [NASA ADS] [Google Scholar]
 Moreno, E., Koenigsberger, G., & Harrington, D. M. 2011, A&A, 528, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ogilvie, G. I. 2014, ARA&A, 52, 171 [Google Scholar]
 Pappalardo, R. T., Belton, M. J. S., Breneman, H. H., et al. 1999, J. Geophys. Res., 104, 24015 [NASA ADS] [CrossRef] [Google Scholar]
 Patel, R., & Penev, K. 2022, MNRAS, 512, 3651 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
 Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [Google Scholar]
 Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
 Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [Google Scholar]
 Penev, K., Sasselov, D., Robinson, F., & Demarque, P. 2007, ApJ, 655, 1166 [Google Scholar]
 Phillips, J. P., RamosLarios, G., Schröder, K. P., & Contreras, J. L. V. 2009, MNRAS, 399, 1126 [NASA ADS] [CrossRef] [Google Scholar]
 Popper, D. M. 1997, AJ, 114, 1195 [Google Scholar]
 Richard, D., & Zahn, J.P. 1999, A&A, 347, 734 [Google Scholar]
 Scharlemann, E. T. 1981, ApJ, 246, 292 [NASA ADS] [CrossRef] [Google Scholar]
 Segatz, M., Spohn, T., Ross, M. N., & Schubert, G. 1988, Icarus, 75, 187 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1988, Adv. Space Res., 8, 135 [NASA ADS] [CrossRef] [Google Scholar]
 Soker, N. 1997, ApJS, 112, 487 [NASA ADS] [CrossRef] [Google Scholar]
 Soker, N., & Rappaport, S. 2001, ApJ, 557, 256 [NASA ADS] [CrossRef] [Google Scholar]
 Spiegel, E. A. 1971, Comments Astrophys. Space Phys., 3, 53 [Google Scholar]
 Terquem, C. 2021, MNRAS, 503, 5789 [NASA ADS] [CrossRef] [Google Scholar]
 Toledano, O., Moreno, E., Koenigsberger, G., Detmers, R., & Langer, N. 2007, A&A, 461, 1057 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Torres, G., Lacy, C. H., Marschall, L. A., Sheets, H. A., & Mader, J. A. 2006, ApJ, 640, 1018 [NASA ADS] [CrossRef] [Google Scholar]
 Vidal, J., & Barker, A. J. 2020, ApJ, 888, L31 [NASA ADS] [CrossRef] [Google Scholar]
 Zahn, J. P. 1966, Ann. Astrophys., 29, 313 [NASA ADS] [Google Scholar]
 Zahn, J. P. 1989, A&A, 220, 112 [Google Scholar]
Appendix A: Polytropic structures, angular velocity, and viscosity for different values of TIDES input parameters
We present the behavior of angular velocity as a function of azimuth angle for the models of the block 1 in the Table 2, both with constant viscosity and the cases with variable viscosity for different radii (Figs. A.1, A.2, A.3, and A.4). The aim is to observe more clearly how the variation of the radius as well as of λ modifies the behavior of ω″. Analogous graphs of the viscosity as a function of azimuth angle calculated for λ = 0.1 and 1 are illustrated in Figs. A.5, A.6, A.7, and A.8. The shape of these functions reveals four maxima per 360° in azimuth angle, each maximum corresponding to the extrema in the angular velocity curves. The maximum value of these computed viscosities^{5} is listed in column 7 of Table 2 for comparison with the constant viscosity values used in KM2016 and shows that they are of similar order of magnitude. However, there is a major difference for the inner layers as mentioned above, where the smaller velocity gradients lead to significantly smaller viscosity values.
Fig. A.1. Angular velocity at the equator as a function of azimuth angle in five layers in the rest frame of the star. Left: Model 2, constant ν; Middle: Model 4, λ = 1; Right: Model 6, λ = 0.1. 
Fig. A.2. Angular velocity at the equator as a function of azimuth angle in five layers in the rest frame of the star. Left: Model 7, constant ν; Middle: Model 8, λ = 1; Right: Model 10, λ = 0.1. 
Fig. A.3. Angular velocity at the equator as a function of azimuth angle in five layers in the rest frame of the star. Left: Model 11, constant ν; Middle: Model 12, λ = 1; Right: Model 14, λ = 0.1. 
Fig. A.4. Angular velocity at the equator as a function of azimuth angle in five layers in the rest frame of the star. Left: Model 23, constant ν; Middle: Model 24, λ = 1; Right: Model 26, λ = 0.1. 
Fig. A.5. Viscosity at the equator as a function of azimuth angle in five layers in the rest frame of the star. Left: Model 4, λ = 1; Right: Model 6, λ = 0.1. 
Fig. A.6. Viscosity at the equator as a function of azimuth angle in five layers in the rest frame of the star. Left: Model 8, λ = 1; Right: Model 10, λ = 0.1. 
Fig. A.7. Viscosity at the equator as a function of azimuth angle in five layers in the rest frame of the star. Left: Model 12, λ = 1; Right: Model 14, λ = 0.1. 
Fig. A.8. Viscosity at the equator as a function of azimuth angle in five layers in the rest frame of the star. Left: Model 24, λ = 1; Right: Model 26, λ = 0.1. 
Table B.1 shows the comparison of the value of the energy dissipation rate in different orbital cycles for the M4 and M72 models whose parameters are the same, but with five and ten layers, respectively. The values shown correspond to layers 1 to 5 for M4 and layers six to ten (the outermost) for M72. Table B.2 presents the values of the rate of energy dissipation in different orbital cycles for all the layers of M72.
The TIDES calculations of Moreno & Koenigsberger (2016) were performed under the assumption of a n = 1.5 polytropic structure which resulted in densities that are significantly larger than those obtained in this paper and which are obtained under the assumption that n = 3 or larger. In Fig. A.9, we show a comparison of the polytropic structures that we explored with the density structure given by the MESA models.
Fig. A.9. Density structure of the MESA models (blue) compared to a polytropic structure (red dash). The MESA models are: Standard model at 0.13 Gyr (h03, top left) and at 8.1 Gyr (h011, top right) as well as the set 2 10x heated model at 10.5 Gyr (h1022, bottom). The polytropes are labeled with their corresponding index, and their central density was chosen to coincide with that of the MESA models for illustration purposes. The density is given in units of g cm^{−3}. 
Appendix B: Comparison of model runs with five and with ten layers
In Table B.1, we compare the energy dissipation rates obtained in the fivelayer computation (model 4) with those obtained in the top five layers of a tenlayer computation (model 72), demonstrating that they are comparable. In Table B.2, we show that the energy dissipation rates in the inner layers have larger cycletocycle variation than those of outer layers. The orbital cycles listed are 150, 200, and 250.
Comparison of a fivelayer and a tenlayer computation.
Comparison of different orbital cycles for model 72.
All Tables
All Figures
Fig. 1. Angular velocity, viscosity and dissipated energy as a function of azimuth angle in the ten layers of model 63. The azimuth angle is measured in the direction of the companion’s orbital motion and φ = 0 corresponds to the subbinary longitude. Top: angular velocity ω″ at the equator in the perturbed star rest frame. The surface layers (largest amplitude curves) are more strongly coupled to the tidal field than the inner layers and thus are forced to lag behind the supersynchronously rotating core. Middle: turbulent viscosity values computed by the model. Bottom: tidal shear energy dissipation rate. 

In the text 
Fig. 2. Properties of the MESA models of set 2, in which the heating profile is introduced continuously throughout the evolutionary trajectory. Top: stellar radius (left) and luminosity (right). The blue curve corresponds to the nonheated models. The orange and green curves correspond, respectively, to the heating profile of model 63 and a heating profile that is ten times larger. Bottom: protonproton reaction rate power (left) and CNO reaction rate power (right). The blue curve corresponds to the nonheated models and its ordinate is on the right side of the plot and given in units of log(L_{⊙}). The orange and green curves correspond to the ratio of rates in the heated and nonheated models, with colors the same as in the top panels. 

In the text 
Fig. 3. Luminosity structure of MESA set 2 models at ages of (top left) 2 Gyr and ages between (top right) 9.257, (bottom left) 9.653 and (bottom right) 11.644 Gyr. The blue curves correspond to the models without added heating; orange to those with heating as given by model 63; green to heating by ten times that given by model 63. The tidally heated models have surface luminosities that are larger than the nonheated models, but with lower luminosities in deeper layers. 

In the text 
Fig. 4. Kippenhahn diagram for the standard model (top) and the set2h10 model (bottom). The panels from left to right represent: Left: evolution from ZAMS to TAMS (defined as abundance of H first drops below 10^{−4} in the stellar core). Center: evolution from the TAMS to 90% of the remaining time. Right: final 10% of the postTAMS evolution. The ordinate is the mass coordinate that runs from the center of the star to the surface. The nuclear energy generation rates are indicated by the different shades of blue. The green hatches indicate zones in which the energy transport is dominated by convection. 

In the text 
Fig. 5. Properties of the MESA models of set 3, in which the heating profile is introduced gradually, starting at the ZAMS and reaching maximum at the TAMS. Top: stellar radius (left) and luminosity (right). The blue curve corresponds to the nonheated models. The orange and green curves correspond, respectively, to the heating profile of model 63 and a heating profile that is ten times larger. Bottom: protonproton reaction rate power (left) and CNO reaction rate power (right). The blue curve corresponds to the nonheated models and its ordinate is on the right side of the plot and given in units of log(L_{⊙}). The orange and green curves correspond to the ratio of rates in the heated and nonheated models, with colors the same as in the top panels. 

In the text 
Fig. 6. Luminosity structure of MESA set 3 models at ages 9.257, 9.653, and 11.644 Gyr. Blue curves correspond to the models without added heating; red to those with heating as given by model 63; green by heating ten times that given by model 63. The tidally heated models have a surface luminosity that is larger than the nonheated models, but lower luminosity in deeper layers. 

In the text 
Fig. 7. Evolutionary tracks on the HRD of the standard model (light blue), the set3h1 model (orange), and the set3h10 model (green) showing the result of introducing a gradual heating near the end of the main sequence. The dark blue curve corresponds to the set2h10 model in which there was constant heating throughout the main sequence. 

In the text 
Fig. A.1. Angular velocity at the equator as a function of azimuth angle in five layers in the rest frame of the star. Left: Model 2, constant ν; Middle: Model 4, λ = 1; Right: Model 6, λ = 0.1. 

In the text 
Fig. A.2. Angular velocity at the equator as a function of azimuth angle in five layers in the rest frame of the star. Left: Model 7, constant ν; Middle: Model 8, λ = 1; Right: Model 10, λ = 0.1. 

In the text 
Fig. A.3. Angular velocity at the equator as a function of azimuth angle in five layers in the rest frame of the star. Left: Model 11, constant ν; Middle: Model 12, λ = 1; Right: Model 14, λ = 0.1. 

In the text 
Fig. A.4. Angular velocity at the equator as a function of azimuth angle in five layers in the rest frame of the star. Left: Model 23, constant ν; Middle: Model 24, λ = 1; Right: Model 26, λ = 0.1. 

In the text 
Fig. A.5. Viscosity at the equator as a function of azimuth angle in five layers in the rest frame of the star. Left: Model 4, λ = 1; Right: Model 6, λ = 0.1. 

In the text 
Fig. A.6. Viscosity at the equator as a function of azimuth angle in five layers in the rest frame of the star. Left: Model 8, λ = 1; Right: Model 10, λ = 0.1. 

In the text 
Fig. A.7. Viscosity at the equator as a function of azimuth angle in five layers in the rest frame of the star. Left: Model 12, λ = 1; Right: Model 14, λ = 0.1. 

In the text 
Fig. A.8. Viscosity at the equator as a function of azimuth angle in five layers in the rest frame of the star. Left: Model 24, λ = 1; Right: Model 26, λ = 0.1. 

In the text 
Fig. A.9. Density structure of the MESA models (blue) compared to a polytropic structure (red dash). The MESA models are: Standard model at 0.13 Gyr (h03, top left) and at 8.1 Gyr (h011, top right) as well as the set 2 10x heated model at 10.5 Gyr (h1022, bottom). The polytropes are labeled with their corresponding index, and their central density was chosen to coincide with that of the MESA models for illustration purposes. The density is given in units of g cm^{−3}. 

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