Structure and evolution of a tidally heated star

The shearing motion of tidal flows that are excited in non-equilibrium binary stars transform kinetic energy into heat via a process referred to as tidal heating. In this paper we aim to explore the way tidal heating affects the stellar structure. We used the TIDES code, which solves the equations of motion of the three-dimensional (3D) grid of volume elements that conform multiple layers of a rotating binary star to obtain an instantaneous value for the angular velocity, $\omega''$, as a function of position in the presence of gravitational, centrifugal, Coriolis, gas pressure, and viscous forces. The released energy, $\dot{E,}$ was computed using a prescription for turbulent viscosity that depends on the instantaneous velocity gradients. The $\dot{E}$ values for each radius were injected into a MESA stellar structure calculation. The method is illustrated for a 1.0+0.8 M$_\odot$ binary system, with an orbital period of $P$=1.44d and departures from synchronous rotation of 5% and 10%. We find that 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. 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.


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 so-called "turbulent viscosity", ν turb , is present, then the internal stellar structure could be affected.
Visiting scholar, Astronomy Department, Indiana University.
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. 2004Mathis et al. , 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 mixing-length 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(Zahn , 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. 2004Mathis et al. , 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 turbulence-induced 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 single-star counterpart.Focusing on the post-main 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 time-variable.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 n-layer calculation of the tidal perturbations, including a prescription for the turbulent viscosity allowing it to be computed as a function of the time-varying and location-dependent 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 mass-transfer 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.

Method
The response of the stellar layers to the external gravitational field of a companion is computed with the n-layer TIDES code1 (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 quasi-hydrodynamic 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 ( 0 r , 0 ϕ , 0 θ ) 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 Runge-Kutta integrator.The secondary is considered to be a point-mass 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 time-marching 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.
A44, page 2 of 16 For the new version of TIDES used in this paper, we implemented a self-consistent 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 n-layer 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 v i = ω i r sin θ its linear velocity.Here, ω i 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 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.

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.
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 correspond- Tolerance for the Runge-Kutta integration ing surface equatorial velocities V rot in the unperturbed star are listed in Cols. 2 and 6, respectively, of Tables 2 and 3.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 low-mass 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 super-synchronous.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   A44, page 3 of 16  Table 2. TIDES five-layer models.
Model 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 five-layer and a ten-layer 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 long-term 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 phase-dependent 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.

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 star2 as a function of the azimuth angle.The sinusoidal shape corresponds to the equilibrium tide.The peak-to-peak 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 peak-topeak 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 one-layer calculations that were performed by KM2016.
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 ten-layer runs for which Ėk=1 correspond to a layer that lies at ∼0.4 R .A layer-by-layer comparison down to this depth clearly illustrates the density-dependence of Ėr , as shown in Table 4 by comparing models 62 and 63.
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 lay- ers.Thus, increasing departures from synchronicity while holding other parameters constant does lead to increasing energy dissipation rates.
Finally, the near-zero-values 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 base-level viscosity associated with the convective eddies are expected to always be present.

Stellar structure models
The stellar structure models were computed with version 15140 of the MESA open-source 1D stellar evolution code3  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: proton-proton 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.(Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 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 time-dependent factor of which increases sharply to unity as t → t TAMS , where it becomes saturated.Table 6 lists properties of the three sets.
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 set2-h10 and set3-h10 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.

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 set2-h1 and set2-h10, which correspond to time-constant 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.
The bottom panels show the total power generated by the proton-proton 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.
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.

Results for gradual injected heat
An overview of the results for gradual time-increasing 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 4 These diagrams were made with the mkipp.pysoftware written by Pablo Marchant https://github.com/orlox/mkipp/A44, page 7 of 16 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.is approached.However, at the TAMS, the set 3 models very rapidly approach the properties of the set 2 models at similar ages.
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 Hertzsprung-Russell Diagram (HRD), as shown in Fig. 7, where we see that the set2-h10 and the set3-h10 tracks are identical after the MS turnoff, which is when both sets are being equally heated.

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, north-south 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 co-rotation.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.

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 R 2 /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 A44, page 8 of 16 .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: proton-proton 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.11.644 Gyr 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.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 5000 5250 5500 5750 6000 Fig. 7. Evolutionary tracks on the HRD of the standard model (light blue), the set3-h1 model (orange), and the set3-h10 model (green) showing the result of introducing a gradual heating near the end of the main sequence.The dark blue curve corresponds to the set2-h10 model in which there was constant heating throughout the main sequence.
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 one-layer 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 tidally-induced 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.

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 A44, page 9 of 16 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 ν −1 turb 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 recently-formed 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 low-mass binary stars in short-period orbits have been determined to be as much as 10% larger than their counterparts in long-period 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 low-mass 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ópez-Morales 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 mass-discrepancy problem may reside in the phenomena we discuss in this paper.

Implications for stellar evolution and population synthesis
The effects of binary interactions on stellar evolution have, until now, focused mainly on the effects of mass-loss 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 short-period 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 mass-transfer (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 mass-loss 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 sub-surface stellar layers would produce the same effects we describe in this paper.

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.
A44, page 10 of 16 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 re-distribute the added energy.

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 post-main 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.(h10-22, 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 . A44

Fig. 1 .
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 sub-binary 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.

Fig. 2 .
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: proton-proton 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.

Fig. 3 .
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.

Fig. 4 .
Fig. 4. Kippenhahn diagram for the standard model (top) and the set2-h10 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 post-TAMS 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.
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: proton-proton 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.
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.

Table 1 .
Description of TIDES input parameters.

Table 4 .
Energy dissipation rate radial profiles from the TIDES ten-layer models.

Table 5 .
Notes.Column 1 lists the layer number; Cols.2, 5, and 7 list the distance of the layer from the star's center; Cols.3, 4, 6, and 8 list the energy dissipation rate for each layer in units of ergs s −1 .Models 62 and 63 have the same radial grid.Heat injection in MESA models.
Notes.The name of each MESA model is listed in Col. 1; Col. 2 indicates whether it includes injected heat.The injected heating profile is indicated in Col. 3 and the manner in which it was injected in Col. 4.

Table 6 .
Ages and events in MESA models.= 2.25 R end of the run h1 11.728 R 1 = 2.25 R end of the run h10 Table B.2.Comparison of different orbital cycles for model 72.In TableB.1,wecompare the energy dissipation rates obtained in the five-layer computation (model 4) with those obtained in the top five layers of a ten-layer computation (model 72), demonstrating that they are comparable.In TableB.2,weshow that the energy dissipation rates in the inner layers have larger cycleto-cycle variation than those of outer layers.The orbital cycles listed are150, 200, and 250.