11institutetext:  Instituto de Ciencias Físicas, Universidad Nacional Autónoma
de México, Ave. Universidad S/N, Chamilpa, Cuernavaca, México
11email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Eccentric binaries: Periastron events and tidal heating

G. Koenigsberger Visiting scholar, Astronomy Department, Indiana University11
D. Estrella-Trujillo 11

(August 2023)

###### Abstract

Context. Binary stars cannot be fully understood without assessing the
interaction effects between the two components and the impact of these effects
on observational diagnostics. Periastron brightening events, also known as the
heartbeat phenomenon, are a clear manifestation of this type of interaction.

Aims. We aim to explore the role of tidal shear energy dissipation in stars
undergoing periastron brightening events.

Methods. We performed a computation from first principles that uses a quasi-
hydrodynamic Lagrangian scheme to simultaneously solve the orbital motion and
the equations of motion of a 3D grid of volume elements covering the inner,
rigidly rotating ”core” of a rotating, the tidally perturbed star. The
equations of motion include the gravitational acceleration of both stars, the
centrifugal, Coriolis, gas pressure accelerations, and viscous coupling
between volume elements. The method is illustrated for a grid of model binary
systems with a 10 M⊙ primary that is perturbed by a 6.97 M⊙ companion in
eccentric orbits ($e$=0 - 0.7). The model was then applied to the heartbeat
star MACHO 80.7443.1718.

Results. We find an increase by factors $\sim$10-6 -10-3 in tidal shear energy
dissipation at periastron, consistent with the majority of observed heartbeat
stars. The magnitude of the periastron effect correlates with the degree of
departure from synchronicity: stars rotating much faster or much slower than
the synchronous rate at periastron present the strongest effect. We confirm
that for eccentricities $\leq$0.3, pseudo-synchronization occurs for
0.8$<\omega/\Omega_{ave}<$1, where $\Omega_{ave}$ is the average orbital
angular velocity. The minimum energy rotation rate (pseudo-synchronism) for
$e$=0.5 and 0.7 occurs for 1.0$<\omega/\Omega_{ave}<$1.15. The tidal shear
energy dissipation model reproduces from first principles the $\sim$23%
maximum brightness enhancement at periastron of MACHO 80.7443.1718.

Conclusions. Our results suggest that the magnitude and shape of the heartbeat
signal may serve as diagnostics for the internal stellar rotation and
turbulent viscosity values.

##  1 Introduction

A growing number of binary systems have been reported to undergo periastron
events, that is, a change in the observable properties around the time of
periastron, such as a brightness enhancement, an increase in surface activity,
and/or a triggering of oscillations. Early examples include Hutchings (1979),
who cited a significant increase in mass-loss rates around periastron in the
two B-type systems, HD 187399 and AZ Cas, and in HD 163181 (B0Ia+0), and van
Genderen & Sterken (2007) at periastron in six systems, including $\eta$
Carinae.

More recently, photometric observations obtained with the Kepler space
telescope have led to the discovery of objects that not only display
periastron-related brightness enhancements, but also low-amplitude sinusoidal
oscillations that tend to have frequencies that are harmonics of the orbital
frequency (Thompson et al., 2012; Fuller, 2017; Guo, 2020). These objects have
come to be known as heartbeat stars, and more than 180 are known from Kepler
alone (Hambleton et al., 2013; Kirk et al., 2016) and many more from BRITE
(Pigulski, 2018), ASAS-SN (Jayasinghe et al., 2020), and TESS (Kołaczek-
Szymański et al., 2021). These stars are all considered to be eccentric orbit
binary systems with typical photometric variations $\Delta$L/L$\leq$10-3.

Although a large majority of the heartbeat stars have masses in the range of
A-type and F-type stars, similar properties are also reported in high-mass OB-
type objects such as $\iota$ Orionis (Pablo et al., 2017) and $\eta$ Carinae
(Richardson et al., 2018) in the Milky Way, MACHO 80.7443.1718 (Jayasinghe et
al., 2019) in the Large Magellanic Cloud, and the massive multiple system HD
5980 (Sterken & Breysacher, 1997; Kołaczek-Szymański et al., 2021) in the
Small Magellanic Cloud. This large diversity points to a common physical
mechanism which is now widely accepted to be the variable gravitational
perturbation caused by the binary companion.111It must be noted that
periastron brightening in massive stars is often interpreted in terms of
colliding winds, a phenomenon that is however hard-pressed to explain low-
amplitude high-frequency post-periastron oscillations.

A companion distorts a star’s intrinsic spherical shape producing two tidal
lobes, one facing the companion and one in the opposite hemisphere. The
extension of the tidal lobes is not necessarily the same on both sides and
their orientation is not necessarily parallel to the line joining the stars’
centers. These depend on the orbital separation as well as the star’s radius,
its rotation velocity, and the viscosity of its material.

In binary systems with circular orbits, the shape of the distorted stars
remains approximately constant throughout the orbital cycle. However, because
the viewing angle to the deformed star changes with orbital phase, the
changing projected size of the disk in the sky leads to brightness variations.
This so-called ellipsoidal effect is a purely geometric effect, and can be
modeled using an appropriate stellar atmosphere model (Orosz & Hauschildt,
2000).

Eccentric binaries differ from circular orbit systems in that the orientation
of the tidal bulges with respect to the line joining the two stellar centers
varies over the orbital cycle, as does their extension. These variations occur
most rapidly around periastron, triggering high frequency dynamical effects
that show up as low-amplitude photometric variations.

Press & Teukolsky (1977) presented a mathematical framework to model the tidal
distortions based on a linear superposition of the non-radial pressure and
gravity mode oscillations, a framework further developed by Kumar et al.
(1995) and Lai (1997), predicting the brightening around periastron and the
appearance of tidally excited oscillations (TEOs) that are harmonics of the
orbital period. Recent advances were achieved by Fuller (2017) to account for
rotation, incorporate nonadiabatic effects, and allow for spin-orbit
misalignment; Schenk et al. (2001) who studied the nonlinear interaction of
modes in rotating Newtonian stars; and by Dewberry & Wu (2024) who proposed an
alternative approach to the treatment of tidally driven oscillations damping
that accounts for dissipative mode coupling.

As with any theoretical framework, the linear mode-decomposition models have
limitations. As pointed out by Press & Teukolsky (1977), a linear mode
calculation is valid if nonlinear effects are very small. Specifically, with
$R_{1}$ the stellar radius and $r_{per}$ the periastron distance, they point
out that the validity of a linear mode calculation is to be doubted if
$(R_{1}/r_{per})^{2}\geq 0.10$. Indeed, using a 3D non-perturbative modeling
method, Fellay et al. (2024) show that the linear approach underestimates
tidal interactions by up to 40% for binary stars with small separations. In
addition, the tidal interaction imposes a large-scale differential-rotation
structure as well as induced oscillatory shearing motions on orbital
timescales (Koenigsberger et al., 2021). An important consequence of the
differential rotation is that the shearing motions can be damped through the
presence of turbulent viscosity, thus converting kinetic energy into heat in a
process known as tidal shear energy dissipation. However, although
observational evidence exists for deviations from an adiabatic approximation
(Guo et al., 2020), with a few exceptions (see, for example, Baruteau &
Rieutord (2013); Astoul et al. (2021)) the current analytical models assume a
rigidly rotating star and thus neglect tidal shear energy heating. Injection
of this added heat into a star’s subphotospheric layers causes it to expand
and become more luminous, at least in the case of a solar-type star (Estrella-
Trujillo et al., 2023).

In the presence of turbulent viscosity, kinetic energy in the tidal flows is
converted into heat at all orbital phases, not only at periastron. However, as
the star approaches and recedes from periastron, the shearing layers undergo
rapid velocity variations leading to a more abrupt release of energy. Moreno
et al. (2011) investigated the potential role of tidal heating as a source of
periastron brightening events, concluding that it could play a significant
role in this phenomenon. However, they were unable to provide quantitative
constraints for three reasons. The first is that the study was limited to a
single surface layer, thus neglecting the effects of underlying layers. The
second is that the model needed as input a turbulent viscosity, values of
which are highly uncertain. The third reason is that the viscosity remained
constant throughout the calculation, whereas in reality turbulent viscosity is
non-isotropic and time variable. These three limitations have now been
overcome (Koenigsberger et al., 2021; Estrella-Trujillo et al., 2023) enabling
this study.

The objective of this paper is to show the general features associated with
tidal shear energy dissipation in eccentric binary systems, particularly the
increased brightness predicted to occur around the time of periastron. The
method is summarized in Section 2. In Section 3 we present the results from a
grid of model calculations for a grid of binary systems having a range of
orbital eccentricities and rotation velocities. In Section 4 we apply the
model to one of the most extreme heartbeat stars found to date, MACHO
80.7443.1718, and illustrate the manner in which our calculation from first
principles directly predicts the observed variability amplitude. A discussion
and the conclusions are presented in Section 5.

##  2 Method

Binary stars are defined to be in equilibrium when the orbit is circular, the
rotation axis coincides with the orbital angular velocity vector, and the
rotation period is equal to the orbital period. A departure from any of these
conditions induces deformations of the stellar surface and excites
oscillations.

Classical tidal theory divides the analytical treatment of the perturbations
into the equilibrium tide model (Darwin 1880; Alexander 1973; Hut 1981, 1982)
and the dynamical tide model (Zahn, 1975; McMillan et al., 1987; Goldreich &
Nicholson, 1989; Lai, 1997; Savonije et al., 1995). The first analyzes the
large-scale distortion imposed by the external gravitational perturber. It
assumes that the viscosity of the material leads to a lag between the
principal axis of the deformed spheroid and the line of centers of the two
stars. This lag is responsible for the torque that over time synchronizes the
rotation and circularizes the orbit, and it is caused by the fact that the
stellar material is not inviscid. The second model treats the star as an
oscillator with a number of modes that are excited by the gravitational
potential of the companion, and assumes that the dissipation of these
oscillations is responsible for the eventual synchronization and
circularization of the orbit (Ogilvie 2014 and references therein). As noted
by Eggleton et al. (1998), rather than competing, each model provides a useful
description of the interaction under different regimes; that is, for different
masses and radii, orbital separations and eccentricities, and rotation rates.
However, both models suffer from important uncertainties in the treatment of
dissipation (Ogilvie, 2014; Mathis et al., 2016; Barker & Astoul, 2021;
Terquem, 2021; Townsend & Sun, 2023; Sun et al., 2023). Both models are
intended for studying the long term evolution of the stellar rotation and
orbital properties of binary stars, with a focus on determining the timescales
over which an equilibrium configuration is reached. However, the uncertainty
in the dissipation limits the possibility of quantitative results.

The approach used in this paper differs significantly from the two classical
models. It makes no a priori assumptions other than that the damping mechanism
can be associated with the viscous stresses suffered by fluid parcels that are
in relative motion with respect to each other. Such relative motions naturally
occur throughout the star given the non-isotropic dependence of the tidal
forcing in the radial and horizontal directions. Thus, our method relies on a
viscosity parameter as it appears in the shear part of the kinematic stress
tensor (e.g., Symon (1971)). The aim is to analyze the short term (orbital
timescales) variations induced by the tidal interaction in order to allow a
direct comparison with observational data.

The method is a computation from first principles which is named TIDES222The
Tidal Interactions with Dissipation of Energy due to Shear code in all its
versions is available upon request to G. Koenigsberger and is easily
implemented in any operating system running a Fortran or GNU Fortran compiler.
It uses a quasi-hydrodynamic Lagrangian scheme which simultaneously solves the
orbital motion of the companion and the equations of motion of a 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 numerical scheme was first applied to the B-type system $\iota$ Orionis
(Moreno & Koenigsberger, 1999). It predicted a rapid increase in the stellar
radius around periastron as well as small-amplitude, short-timescale
oscillations throughout the orbit. It was subsequently used to produce
synthetic absorption line profiles as a function of the orbital phase in
$\epsilon$ Persei (Moreno et al., 2005) and $\alpha$ Virginis (Harrington et
al. 2009; 2014; Palate et al. 2013b) and other systems (Palate et al., 2013b).
The simulated line-profile variability captures many of the characteristics of
the observations, namely, the spectroscopic equivalent of the ellipsoidal
effects and the traveling “bumps” that are characteristic of tidally induced
non-radial pulsators. The method was also used to derive the synchronization
timescales for eccentric binaries (Toledano et al., 2007; Moreno et al.,
2011), finding good agreement with the expression given by Zahn (2008) (his
Eq. 2.5) for viscous dissipation for circular and $e$=0.1 orbits, the cited
range of applicability.

The above calculations were performed for a single surface layer contiguous to
the rigidly rotating internal region. The output includes the tidal
perturbations in the radial and azimuthal directions and it was shown that the
tidal velocity in the radial direction, ${\mathrm{v}}_{r}$, is generally an
order of magnitude smaller than that in the azimuthal direction,
${\mathrm{v}}_{\varphi}$, consistent with the analytical results for the
equatorial latitudes333The inviscid approximation leads to poorer results for
the tidal amplitude at latitudes closer to the poles because it neglects
angular momentum to these regions from the strongly perturbed equatorial
latitudes, see Harrington et al. (2009), Fig. 7. of Scharlemann (1981),
derived by assuming small departures from corotation. The dominance of the
azimuthal component of the tidal flows was also found by Lee (2017). This
original version of the numerical method is named TIDES-1. It requires a
viscosity value as input and its output includes the radial and azimuthal
velocities for the surface layer. These velocities on the hemisphere facing
the observer are then projected onto the line of sight allowing the
computation of the orbital phase-dependent profile variability of photospheric
absorption line profiles (Moreno et al., 2005).

The extension of the TIDES method to $n$ layers (TIDES-n) was introduced in
Koenigsberger et al. (2021) where, for circular-orbits, it was shown that the
rotation rate tends toward a uniform rotation structure when synchronization
is attained. However, the synchronization timescale is different for each
layer, because it goes as the viscous timescale, $\sim\nu$/r2, where $\nu$ is
the viscosity and $r$ is the radius of the particular layer. The calculation
in Koenigsberger et al. (2021) showed how each layer procedes toward
synchronization and that the average angular velocity of the surface layer (at
the equator) is the first to reach synchronicity, in agreement with the
conclusion reached analytically by Goldreich & Nicholson (1989). TIDES-n
requires a turbulent viscosity value as input. It is isotropic and remains
constant throughout the calculation.

The primary results presented in this paper use the latest version of the code
(TIDES-nvv). It is based on the same algorithm as TIDES-n but has the
additional capability of computing a turbulent viscosity at each grid point
(Estrella-Trujillo et al., 2023). It outputs the azimuthal velocity, turbulent
viscosity and corresponding energy dissipation rate. These variables are all
non-isotropic and time-variable.

The strongest approximation in our calculation concerns the mathematical
formulation for the turbulent viscosity. Turbulent viscosity arises in
unstable flows that are characterized as turbulent. Mathis et al. (2004) noted
that turbulence might be expected even in stellar radiative zones due to
nonlinear processes that appear at large Reynolds numbers. Understanding
turbulent viscosity is a complex problem that falls in the realm of fluid
dynamics and turbulent transport (see, for example, Le Bars et al. (2022);
Zahn (1992); Richard & Zahn (1999); Zahn (2008); Mathis et al. (2004, 2016);
Garaud & Kulenthirarajah (2016)). We use the simplest functional
representation, as given by (Landau & Lifshitz, 1987), in which its value is
directly proportional to the instantaneous gradient in velocity of neighboring
fluid parcels. This representation yields a non-isotropic and time-dependent
turbulent viscosity value, which is a significantly less arbitrary approach
than simply adopting a constant value. Furthermore, the viscosity thus derived
includes both the shear due to the global differential rotation and the shear
due to the tidally excited oscillations.

The need for a turbulent viscosity in massive star models follows from the
fact that rotation induced shear is required to account for the enhanced
abundances of nuclear-processed elements on the surface of main sequence
massive stars (Maeder & Meynet, 1996; Talon et al., 1997; Maeder & Meynet,
2000). The fluid velocity used in these formulations is associated with
meridional circulation currents and the differential effect of the Coriolis
force (Hirschi & Maeder, 2010). In binary stars with asynchronous rotation,
the tidal perturbations can induce significantly larger fluid velocities and
velocity gradients (Koenigsberger et al. 2021). Also, observational data of
binary systems in clusters indicates that systems with orbital periods shorter
than $\sim$6-8 days have mostly circular orbits, which means that the
circularization process must occur on much shorter timescales than possible in
the presence of only molecular viscosity. The circularization timescale that
is usually applied in this case is as derived by Zahn (1977), with typical
turbulent viscosities in the range 1010cm2 s-1-1013cm2 s-1 (Zahn, 2008). It is
also found that turbulent viscosity as estimated using arguments from the
mixing length theory should have values on the order of 1012-1013 cm2/s
(Käpylä et al., 2020).

It is interesting to note that viscosity values that are orders of magnitude
larger than the molecular (kinematical) viscosities (1-103 cm2/s) appear in
many additional contexts. For example, numerical simulations of the solar
convective region use a viscosity value 1010-1012 cm2/s (Elliott et al., 2000)
based on the idea that the dynamics are expected to be governed principally by
the largest (and most energetic) eddies, and not to be sensitive to small-
scale flow features. 444Elliott et al. (2000) note, however, that the large
viscosity used in their models should be regarded as an artificial viscosity
whose values are chosen to suitably truncate the nonlinear energy cascade and
thereby prevent nonlinear instability. This is in contrast to the value of
kinematic viscosity $\sim$27 cm2/s at the solar center discussed in Fig. 3 of
Caleo et al. (2016), where the plasma conditions are significantly different
from those in the convective zone. Also, Pérez-de-Tejada (1999) estimated
$\nu$=6. 1013 cm2/s for this parameter in boundary layers around solar system
planetary ionospheres; Efroimsky (2018) used the outgoing energy flux due to
the vapor plumes in Saturn’s moon Enceladus to estimate that its mean tidal
viscosity $\eta$=$\nu\rho\sim$0.24 1014 Pa s (which translates into $\nu\sim$2
1014 cm2/s assuming a mass density $\rho\sim$0.9 for ice); and Shi et al.
(2006) found $\eta\sim$ 1019 Pa-s for Earth’s crustal medium at a depth of 32
km.

###  2.1 Definition of parameters

The star whose tidal perturbations are being considered is called the
”primary”, and its mass and unperturbed radius are denoted, respectively,
$m_{1}$ and $R_{1}$. The companion is assumed to be a point source and its
mass is labeled $m_{2}$. The primary’s initial rotation is uniform and the
rotation axis is perpendicular to the orbital plane. A number $N_{r}$ of
stellar layers are allowed to respond to the forces in the system. This leads
to a departure from the initial uniform rotation structure as it evolves over
time.

We characterize the initial rotation velocity in terms of the synchronicity
parameter $\beta_{0}$=$\omega_{0}/\Omega_{0}$, where $\omega_{0}$ is the
initial rotation angular velocity and $\Omega_{0}$ is the orbital angular
velocity at a reference point. For eccentric orbits, periastron is a
convenient reference point, so that $\beta_{0}$ can be cast in the form:
$\beta_{0}$=0.02P$\frac{V_{rot}(1-e)^{3/2}}{R_{1}(1+e)^{1/2}}$. Here, $P$ is
the orbital period in days, $V_{rot}$ is the initial equatorial rotation
velocity of the uniformly rotating star (initial condition) in km/s, and
$R_{1}$ is the stellar radius in solar units.

Table 1: Description of input parameters.

Parameter | Value | Description
---|---|---
Porb | (a) | Orbital period (days)
e | (a) | Orbital eccentricity
m1 | 10 | Primary mass (perturbed star) (M⊙)
m2 | 6.97 | Secondary mass (point source) (M⊙)
R1 | 5.48 | Undisturbed primary radius (R⊙)
$\beta_{0}$ | (a) | Initial synchronicity parameter
$\lambda$ | 1 | Turbulent viscosity coefficient
IND | 3 | Polytropic index
$\Delta$R/R1 | 0.06 | Layer thickness
Nr | 10 | Number of layers
Nφ | 200 | Number of partitions in longitude
Nθ | 20 | Number of partitions in latitude
Tol | 10-7 | Tolerance for Runge-Kutta integration
555(a) See Tables 2, 4, and 3.

Table 2: Overview of the model grid.

Set | $e$ | P | $a$ | Ncy | Nphases
---|---|---|---|---|---
E0.0 | 0.0 | 6.00 | 35.7 | 30, 50 | 10
E0.1 | 0.1 | 6.00 | 35.7 | 30, 50 | 10
E0.3 | 0.3 | 8.75 | 45.9 | 30, 50 | 130
E0.5 | 0.5 | 14.49 | 64.3 | 12, 20 | 130
E0.7 | 0.7 | 31.18 | 107.2 | 5, 10 | 169
666$e$ is the orbital eccentricity, the orbital period P is given in days, and
the orbital semimajor axis $a$ is given in R⊙. Ncy is the number of orbital
periods at which the results were sampled and Nphases is the number of orbital
phases at each of the sampled cycles for which output was written.

The rate of energy dissipation per unit volume is as given in Moreno et al.
(2011):

$\dot{E}_{V}\simeq{\nu\rho}\left\\{\frac{4}{3}\left(\frac{\partial\omega}{\partial{\varphi}}\right)^{2}+\left[r^{2}\left(\frac{\partial\omega}{\partial
r}\right)^{2}+\left(\frac{\partial\omega}{\partial\theta}\right)^{2}\right]{\sin}^{2}{\theta}\right\\},$
(1)

with $\nu$ the kinematical viscosity, $\rho$ the mass density, $\omega$ the
angular velocity, and $r$, $\theta$, $\varphi$, respectively, the radius,
latitude, and longitude coordinates. The units of $\dot{E}_{V}$ are energy per
unit time per unit volume.

The viscosity is assumed to be a turbulent viscosity and is calculated at each
timestep and volume element interface using the expression (Landau & Lifshitz,
1987):

$\nu_{\mathrm{turb}}=\lambda\ell_{\mathrm{t}}\Delta u_{\mathrm{t}},$ (2)

where $\ell_{\mathrm{t}}$ is the characteristic length of the largest eddies
that are associated with the turbulence, here assumed to be on the order of
the layer thickness, and $\Delta u_{\mathrm{t}}$ is the typical average
velocity variation of the flow over the length $\ell_{\mathrm{t}}$. The
variable $\lambda\in$[0,1] is a proportionality parameter that describes the
fraction of kinetic energy transformed into thermal energy. We used
$\lambda$=1. Inspection of Eqs. 1 and 2 and numerical experiments (Estrella-
Trujillo et al. 2024, in preparation) indicate that the energy dissipation
rate for different values scales approximately with the value of $\lambda$ and
is thus straightforward to estimate for smaller values. The above prescription
is applied to the interface between each volume element and its neighboring
elements. Thus, the non-isotropic and time-variable nature of the viscosity is
incorporated in the calculation.

The density is obtained assuming a polytropic stellar structure. Each volume
element of the grid contains a fixed mass of gas that obeys the ideal gas
equation of state, which provides the corresponding gas pressure.

The total energy dissipation rate at particular times during the orbital cycle
is:

$\dot{E}_{tot}(t)=\int_{V}\dot{E}_{V}dV.$ (3)

Another quantity that is of interest, especially when comparing with
investigations of tidal theory, is the energy dissipation rate averaged over
the orbital cycle:

$\left<\dot{E}_{tot}\right>=\frac{\int_{0}^{P}\dot{E}_{tot}(t)dt}{P}.$ (4)

Table 3: Eccentric models.

Num | $\beta_{0}$ | Vrot | $<\dot{E}_{tot}>$ | $\dot{E}_{tot}^{per}$ | $\dot{E}_{tot}^{ap}$ |  $\Delta\dot{E}$/L∗
---|---|---|---|---|---|---
P=6.0 $e$=0.1 $r_{per}$=32.15 |
201 | 0.20 | 11 | 7.60 10-1 | 1.65 | 4.10 10-1 | 3.10 10-3
202 | 0.40 | 22 | 2.30 10-1 | 5.40 10-1 | 1.20 10-1 | 1.05 10-3
203 | 0.60 | 34 | 4.80 10-2 | 1.40 10-1 | 2.00 10-2 | 3.00 10-4
204 | 0.80 | 45 | 4.20 10-3 | 1.20 10-2 | 1.20 10-3 | 2.70 10-5
205 | 0.88 | 49 | 2.50 10-3 | 4.80 10-3 | 2.00 10-3 | 7.00 10-6
206 | 0.95 | 53 | 4.20 10-3 | 2.40 10-3 | 5.50 10-3 | -7.75 10-6
207 | 1.00 | 56 | 8.80 10-3 | 3.40 10-3 | 1.10 10-2 | -1.90 10-5
208 | 1.05 | 59 | 1.50 10-2 | 7.10 10-3 | 2.00 10-2 | -3.22 10-5
209 | 1.20 | 67 | 6.30 10-2 | 4.90 10-2 | 7.60 10-2 | -6.75 10-5
210 | 1.40 | 79 | 2.70 10-1 | 2.70 10-1 | 3.00 10-1 | -7.50 10-5
211 | 1.60 | 90 | 9.30 10-1 | 9.30 10-1 | 8.80 10-1 | 1.25 10-4
212 | 1.88 | 105 | 3.01 | 3.48 | 3.12 | 9.00 10-4
213 | 1.95 | 109 | 4.20 | 4.84 | 4.32 | 1.30 10-3
P=8.73 $e$=0.3 $r_{per}$=32.1 |
301 | 0.20 | 12 | 6.50 10-1 | 2.47 | 1.80 10-1 | 5.73 10-3
302 | 0.40 | 24 | 1.80 10-1 | 6.90 10-1 | 3.50 10-2 | 1.64 10-3
303 | 0.60 | 37 | 4.90 10-2 | 2.30 10-1 | 1.00 10-2 | 5.50 10-4
304 | 0.80 | 49 | 2.40 10-2 | 7.90 10-2 | 1.10 10-2 | 1.70 10-4
305 | 0.88 | 54 | 1.60 10-2 | 1.70 10-2 | 9.60 10-3 | 1.85 10-5
306 | 0.95 | 58 | 1.30 10-4 | 7.60 10-3 | 1.30 10-2 | -1.35 10-5
307 | 1.00 | 61 | 2.20 10-2 | 6.70 10-3 | 1.60 10-2 | -2.33 10-5
308 | 1.05 | 64 | 2.40 10-2 | 5.60 10-3 | 1.70 10-2 | -2.85 10-5
309 | 1.20 | 73 | 5.30 10-2 | 3.00 10-2 | 3.60 10-2 | -1.50 10-5
310 | 1.40 | 86 | 1.70 10-1 | 2.10 10-1 | 1.20 10-1 | 2.25 10-4
311 | 1.60 | 98 | 6.00 10-1 | 9.50 10-1 | 4.80 10-1 | 1.17 10-3
312 | 1.88 | 115 | 4.42 | 5.60 | 3.64 | 4.90 10-3
313 | 1.95 | 119 | 6.45 | 8.11 | 5.53 | 6.43 10-3
P=14.49 $e$=0.5 $r_{per}$=32.1 |
501 | 0.20 | 13 | 5.40 10-1 | 2.83 | 3.62 | -1.97 10-3
502 | 0.40 | 26 | 2.20 10-1 | 1.04 | 3.00 10-2 | 2.52 10-3
503 | 0.60 | 39 | 1.30 10-1 | 4.30 10-1 | 6.00 10-2 | 9.25 10-4
504 | 0.80 | 52 | 4.30 10-2 | 9.00 10-2 | 2.00 10-2 | 1.75 10-4
507 | 1.00 | 66 | 2.40 10-2 | 1.10 10-2 | 4.72 10-3 | 1.57 10-5
508 | 1.05 | 69 | 2.50 10-2 | 8.93 10-3 | 5.53 10-3 | 8.50 10-6
509 | 1.10 | 72 | 2.80 10-2 | 1.03 10-2 | 6.76 10-3 | 8.85 10-6
510 | 1.20 | 78 | 3.90 10-2 | 3.01 10-2 | 8.38 10-3 | 5.43 10-5
511 | 1.40 | 92 | 1.10 10-1 | 2.40 10-1 | 3.09 10-2 | 5.23 10-4
512 | 1.60 | 105 | 4.30 10-1 | 1.13 | 1.60 10-1 | 2.43 10-3
513 | 1.88 | 123 | 4.65 | 8.18 | 3.55 | 1.16 10-2
514 | 1.95 | 128 | 8.80 | 13.2 | 7.42 | 1.44 10-2
P=31.17 $e$=0.7 $r_{per}$=32.15 |
701 | 0.20 | 14 | 3.60 10-1 | 3.07 | 4.96 10-3 | 7.66 10-3
702 | 0.40 | 28 | 2.00 10-1 | 1.40 | 3.86 10-2 | 3.40 10-3
703 | 0.60 | 42 | 1.20 10-1 | 5.10 10-1 | 5.89 10-2 | 1.13 10-3
704 | 0.80 | 56 | 4.70 10-2 | 1.20 10-1 | 1.12 10-2 | 2.72 10-4
707 | 1.00 | 70 | 2.20 10-2 | 2.43 10-2 | 3.02 10-3 | 5.32 10-5
708 | 1.05 | 73 | 2.06 10-2 | 2.09 10-2 | 2.87 10-3 | 4.51 10-5
709 | 1.10 | 77 | 2.05 10-2 | 2.40 10-2 | 2.69 10-3 | 5.33 10-5
711 | 1.40 | 98 | 6.90 10-2 | 4.20 10-1 | 1.11 10-2 | 1.02 10-3
712 | 1.60 | 112 | 3.00 10-1 | 2.21 | 8.33 10-2 | 5.32 10-3
713 | 1.88 | 131 | 3.71 | 16.2 | 2.24 | 3.49 10-2
714 | 1.95 | 136 | 5.66 | 22.6 | 3.70 | 4.72 10-2
777Surface equatorial rotation velocity Vrot is in km/s; total energy
dissipation rate averaged over the orbital cycle, total energy dissipation
rate at periastron and at apastron are in units of 1035 ergs/s. The last
column gives ($\dot{E}_{tot}^{per}$-$\dot{E}_{tot}^{ap}$)/L∗ where
L∗=4.$\times$1037 ergs/s is the adopted total luminosity for a 10⊙ star such
as the one we model. A negative value occurs when the brightness at apastron
is greater than at periastron. The heading for each set models lists the
corresponding orbital period in days, eccentricity, and orbital separation at
periastron $R_{per}$ in solar units.

Table 4: Circular orbit set.

Model | $\beta_{0}$ | Vrot | $<\dot{E}_{tot}>$ |  $\Delta\dot{E}$/L∗
---|---|---|---|---
P=6 d, e=0
102 | 0.40 | 18 | 0.24 | 6.0$\times$10-4
103 | 0.60 | 27 | 4.2$\times$10-2 | 1.0$\times$10-4
104 | 0.80 | 36 | 5.8$\times$10-3 | 1.4$\times$10-5
106 | 0.95 | 43 | 7.9$\times$10-5 | 2.0$\times$10-7
108 | 1.05 | 48 | 6.8 $\times$10-5 | 1.7$\times$10-7
109 | 1.20 | 55 | 4.3 $\times$10-3 | 1.1$\times$10-5
111 | 1.60 | 73 | 0.17 | 4.3$\times$10-4
113 | 1.95 | 89 | 0.92 | 2.3$\times$10-3
888Surface equatorial rotation velocity Vrot is given in km/s; total energy
dissipation rate averaged over the orbital cycle in units of 1035 ergs/s.

###  2.2 Input parameters

The TIDES input parameters are described in Table 1, where we also list the
values that were kept constant. The masses and the primary star’s radius were
chosen to be similar to those used in Harrington et al. (2009), Harrington et
al. (2016), and Koenigsberger et al. (2021). The eccentricities ($e$=0, 0.1,
0.3, 0.5, 0.7) were chosen to sample the range observed in heartbeat stars
reported by Kolaczek-Szymánski et al. (2021). The orbital periods were chosen
so as to keep the periastron distance at an approximately constant value,
$r_{per}$=32.1 R⊙. The synchronicity parameters $\beta_{0}$=0.2 - 1.95 were
chosen to sample the range in values of the eclipsing binaries studied by
Lurie et al. (2017) (their Fig. 7)999It is important to note that the rotation
velocities quoted in Lurie et al. (2017) were determined based on the
assumption that the light variations are caused by starspots.

Table 2 provides an overview of the parameters that were held constant for
each eccentricity set: Orbital period P, orbital semi-major axis $a$, number
of orbital periods and phases that were sampled, Ncy and Nphases,
respectively. Detailed descriptions for each model are provided in Table 3 for
the eccentric orbit models and Table 4 for the circular orbit models. Column 1
lists the model number, column 2 the value of $\beta_{0}$ and column 3 the
surface equatorial rotation velocity. This velocity corresponds to the initial
equatorial rotation velocity of the uniformly rotating star.

##  3 Results

An initial uniform rotation state transitions into one in which the average
rotation rate of each layer is either faster (if $\beta_{0}<1$) or slower (if
$\beta_{0}>1$) than the initial rate (Koenigsberger et al., 2021). In all
cases, the surface layer changes its velocity faster than any of the other
layers. As a result, the average angular rotation is a function of distance
from the stellar center. Furthermore, because the equatorial latitude is the
one most strongly affected by the tidal force, it slows down or speeds up
(depending on the $\beta_{0}$ value) faster than the polar regions. Thus,
there is also an average differential rotation structure in the polar
direction. Both of these differential rotation structures contribute to the
total energy dissipation rate when a turbulent viscosity is present. In
addition, the models rapidly develop the tidal bulges with the corresponding
inflow and outflow of material. This relative motion of different mass parcels
also contribute to the total energy dissipation rate.

###  3.1 Pseudosynchronicity

The energy dissipation rate averaged over an orbital period
$\left<\dot{E}_{tot}\right>$ is listed in Table 3 (column 3), as are the
instantaneous energy dissipation rates $\dot{E}_{tot}$ at periastron (column
4) and apastron (column 5). This table clearly shows that larger departures
from synchronicity lead to larger $\left<\dot{E}_{tot}\right>$ values. As
shown in Fig. 1, minimum dissipation rates occur near the synchronization
rotation rate $\beta_{0}\simeq$1 but, except for $e$=0, the actual minimum
does not occur at $\beta_{0}$=1 but at the pseudosynchronous rate,
$\beta_{ps}$, as shown analytically by Hut (1981). However, contrary to the
analytical result which predicts a relatively constant $\beta_{ps}$ value for
$e\geq$0.2, our $\beta_{ps}$ values depend on eccentricity. Another difference
is that although the analytical result predicts a unique $\beta_{ps}<$1 value
for the minimum energy rotation state, our $e$=0.5 and 0.7 models show a flat
and extended minimum at $\beta_{0}>$1\. A similar result was found by Moreno
et al. (2011) based on denser $\beta_{0}$ grids. The flat minimum is also
consistent with the recent finding of a blurred pseudosynchronization in
numerical calculations of the high eccentricy binary KOI-54 (Townsend & Sun,
2023).

Figure 1: Energy dissipation rate as a function of the degree of synchronicity
$\beta_{0}$. Colors correspond to the eccentricity of the model, as indicated
in the legend. Left: Average over an orbital cycle,
$\left<\dot{E}_{tot}\right>$. Right: Value at periastron,
$\dot{E}_{tot}^{per}$. The vertical line indicates $\beta_{0}$=1,
corresponding to the synchronous rotation rate for a circular orbit system.

A bifurcation between subsynchronous and supersynchronous rotation is present
in the perturbed rotation velocity. This is quantified by finding the
maximum/minimum angular velocity at each latitude and layer depth, and is
defined as $\pm|\omega_{max}|$. That is, the value of the maximum departure
from the rest frame velocity, with the sign indicating whether it is in the
direction of rotation (positive sign) or in the opposite direction. In the
subsynchronous rotators, $|\omega_{max}|>$0, as expected since the action of
the tidal force is to accelerate the stellar rotation toward corotation. The
opposite is true for the supersynchronous rotators. Examples of this behavior
for the surface layer at the equator are shown in Figs. 9 and 10. These show
that when rotation is close to the pseudosynchronous rate, $|\omega_{max}|$
displays high frequency oscillations in addition to the rapid excursions to
either large positive or negative values around periastron. For systems with
circular orbits, the energy dissipation rate decreases symmetrically around
$\beta_{0}$=1, tending toward a null value at synchronicity as expected (Table
4).

###  3.2 Orbital phase-dependence and periastron events

Fig. 2 illustrates the orbital phase-dependent behavior of $\dot{E}_{tot}$ in
the $e$=0.3 models. Several trends are clearly visible: a) there is always an
increase in $\dot{E}_{tot}$ around the time of periastron; b) maximum
intensity tends to occur shortly after periastron; c) the durations of the
periastron event is shorter in the subsynchronous models than in the
supersynchronous models; d) models with $|\beta_{0}-1|\simeq$1 present a dip
at periastron; and e) models with $|\beta_{0}-1|\simeq$1 display high
frequency oscillations. The orbital phase-dependent behavior for $e$=0.5 and
$e$=0.7 is similar to that of $e$=0.3, as shown in Fig. 3, but the amplitude
of the rise in brightness around periastron increases with increasing
eccentricity.

Assuming that the primary star in our models is similar to the B-type star in
the Spica system, we can adopt $L_{*}$=4. 1037 ergs/s (Tkachenko et al., 2016)
and compare the brightness enhancements at periastron due to tidal shear
energy dissipation to this luminosity. We list in Column 6 of Table 3 the
difference between the maximum around periastron and the mean value at
apastron normalized to $L_{*}$:
$\Delta\dot{E}_{tot}$/L∗=($\dot{E}_{tot}^{per}-\dot{E}_{tot}^{ap}$)/L∗. We
find values mostly in the range 10-5 to 10-3, with a few as low as 10-6 and
others as high as 10-2, the latter in models with the highest rotation rates
and largest eccentricity. Values in the range 10-6-10-3 are typically found in
the heartbeat stars.

A curious feature in Column 6 of Table 3 is that some of the values of
$\Delta\dot{E}_{tot}$/L∗ are negative. This is because of the dip at
periastron which, in models with $|\beta_{0}-1|\simeq$1 can have a lower
energy around apastron than the value of the mean.

The value of viscosity $\nu_{turb}$ plays a crucial role in determining the
energy dissipation rate. Our results confirm the idea that $\nu_{turb}$ has a
time-dependence in addition to a dependence on spatial coordinates. For
example, we find $\nu$ in the range 1012 cm2/s - 1015 cm2/s for the inner and
surface layers, respectively (see Table 6) at periastron and smaller values
around apastron.

The turbulent viscosities we derive for inner layers are on the order of those
found by Penev et al. (2007) and others (Zahn, 2008) 1010 \- 1013 cm2/s, but
the surface values are likely to be overestimates, an issue that could be
resolved with a variable $\lambda$ value if a straightforward theoretical
formulation existed for predicting the rate of conversion from the tidal
kinetic energy into heat. This remains an open question.

Figure 2: Total energy dissipation rate (in units of 1035 ergs/s) as a
function of the orbital phase for the $e$=0.3 set of models. Models with
$\beta_{0}<$0.80 (blue) show a relatively narrow periastron effect with a
single peak centered in phase slightly after periastron passage. Models with
$\beta_{0}\geq$1.20 show a significantly broader periastron effect with a
structured peak. Models with 1.20$\leq\beta_{0}>$0.80 tend to have a dominant
peak after periastron and multiple weaker peaks at other phases. Left: All
models of the set; Right: Models with 0.60$\geq\beta_{0}\leq$1.40, which have
the weakest energy dissipation rates and are plotted on a scale that shows the
variations associated with the tidally excited oscillations.

Figure 3: Total energy dissipation rate (in units of 1035 ergs/s) as a
function of the orbital phase for the $e$=0.5 (top) and $e$=0.7 (bottom)
models. All $\beta_{0}$ values modeled are plotted on the left while on the
right we zoom in on the 0.80$\leq\beta_{0}\leq$1.40 cases. Although the
orbital elements are all the same in these models, the perturbation duration
around periastron differs, being longer for faster rotation rates (larger
$\beta_{0}$ values). The lowest-lying curves are those for $\beta_{0}$=1.00
(green), 1.05 (red dash) and 1.10 (red dot), all of which display several
discrete peaks between orbital phase 0.9-1.3.

###  3.3 Conclusions from this section

We find that the pseudosynchronous rotation rate $\beta_{ps}$ for $e$=0.1 and
0.3 is slower than the synchronous rate at periastron while that for $e$=0.5
and 0.8 is faster. Contrary to the results of Hut (1981), $\beta_{ps}$ is not
approximately constant for eccentricities $e>$0.2 (see his Fig. 3, the curve
labeled $\Omega_{ps}$/$n_{p}$, which is equivalent to our $\beta_{ps}$). A
possible source for this discrepancy may be the assumption of uniform rotation
in Hut (1981) while our results are obtained for a star that has differential
rotation. Our results are consistent with those of Moreno et al. (2011) who
found that $\beta_{ps}$ does not have a unique value for $e$=0.5 and 0.7 but
that the minimum energy state occurs for a range of rotation rates near the
synchronous rate at periastron. This broad $\beta_{ps}$ region is also
consistent with the blurred pseudosynchronization found in the numerical
calculations of Townsend & Sun (2023).

There are clear trends in the orbital phase-dependent behavior of the energy
dissipation rates visible in all eccentricities that we modeled: a) maximum
brightness tends to occur slightly after periastron, not at periastron; b) the
duration over which an increased brightness may be observed is shorter in the
subsynchronous models than in the supersynchronous models; c) for fixed values
of orbital period and eccentricity, the energy dissipation rate increases with
increasing departure from synchronicity, with the largest values found for
$\beta_{0}>>$1; d) for fixed values of $\beta_{0}$ and periastron distance,
the energy dissipation rate increases with increasing eccentricity; e) models
with $|\beta_{0}-1|\simeq$1 present a dip in brightness at periastron as well
as prominent high frequency variations that we identify with tidally excited
oscillations.

Table 5: Input parameters and results MACHO 80.7443.1718.

Parameter | Value | Notes
---|---|---
Data from observations and input parameters
Porb | 32.83627 d | (a)
e | 0.507 | (a)
m1 | 34.5 M⊙ | (b)
m2 | 15.7 M⊙ | (b)
R1 | 23.7 R⊙ | (b)
Vrot | 250 km/s | (c)
L∗ | 1.57$\times$1039 ergs/ | (b)
$\Delta$Fperi/F |  $\sim$23% | (d)
$\beta_{0}$ | 1.95 |
$\lambda$ | 1 |
Ncy | 100 |
Nr | 5 |
$a$ | 159 R⊙ |
$r_{per}$ | 78.5 R⊙ |
Output values and results
$\dot{E}_{peri}$ | 3.7$\times$1038 ergs/s | (e)
$\dot{E}_{apo}$ | 6.6$\times$1036 ergs/s | (e)
$\left<\dot{E}_{tot}\right>$ | 6.5$\times$1037 ergs/s | (e)
$\Delta\dot{E}_{peri}$/L∗ | 23.5% | (f)
$R_{max-peri}$ | 23.83 R⊙ | (g)
$R_{min-per}$ | 23.63 R⊙ | (g)
$R_{max-ap}$ | 23.76 R⊙ | (g)
$R_{min-ap}$ | 23.66 R⊙ | (g)
$\left<R_{per}\right>$ | 23.71 R⊙ | (h)
$\left<R_{ap}\right>$ | 23.70 R⊙ | (h)
101010(a) Jayasinghe et al. (2019); (b) Jayasinghe et al. (2021); (c) Adopting
an orbital inclination $i$=44∘ and ${\mathrm{v}}\sin(i)$=174 km/s from
Jayasinghe et al. (2021); (d) estimated from Fig. 9 of Jayasinghe et al.
(2021) which shows maximum light at levels of $\sim$1.2 to $\sim$1.26
(depending on the TESS sector) above the average value $\sim$1 between
subsequent periastron passages. ; (e) Output from TIDES-nvv; (f)
($\dot{E}_{tot}^{peri}$-$\dot{E}_{tot}^{apo}$)/L∗; (g) Maximum and minimum
equatorial radius at periastron and apastron, including high frequency
oscillations, from the TIDES-1 version of the code with $\nu$=0.028
R${}^{2}_{\odot}$/d. The corresponding values from the TIDES1-vv model with
$\lambda$=0.01 are 23.76 and 23.67 R⊙. (h) average polar radius at periastron
and apastron from the TIDES-1 version of the code.

##  4 Modeling the heartbeat star MACHO 80.7443.1718

MACHO 80.7443.1718 is located in the Large Magellanic Cloud Lucke-Hodge 58
association, which contains stars formed within the past few million years
(Garmany et al., 1994). It was classified as an eclipsing binary system
(Alcock et al., 1997) and has recently been found to display a $\sim$40% peak-
to-peak variability amplitude with morphology that places it among the
heartbeat-type stars (Jayasinghe et al., 2019). An analysis of the light curve
and radial velocity curve of photospheric absorption lines was performed by
Jayasinghe et al. (2021), leading to the conclusion that minimum light
coincides with the time of periastron passage, defined to occur at orbital
phase $\varphi$=0, and maximum light at $\varphi\sim$0.04. The times of minima
in the light curve were found to undergo a systematic change by 11.1 s/yr
(Kołaczek-Szymański et al., 2022). If interpreted to represent times of
central eclipse, this implies an extraordinarily large rate of change in the
orbital period.

The observed brightness increase $\Delta$F/F$\sim$23% at periastron far
exceeds the typical values in the range 10-6 \- 10-3 usually found in
heartbeat stars (Fuller, 2017). MacLeod & Loeb (2023) pointed out that the
stellar radius of 24 R⊙ deduced from the spectroscopic observations
(Jayasinghe et al., 2021) was inconsistent with the significantly larger tidal
distortion needed by standard theories to explain such strong light curve
variations. In order to address this problem, MacLeod & Loeb (2023) performed
a hydrodynamic simulation that predicts a tidal wave raised at periastron
having an extension of several solar radii, and concluded that the breaking
wave could lead to shock-heating on the surface of the star. These results,
however, are based on several assumptions, the strongest of which is that the
core of the primary star is excised from the calculation and treated as a
point mass leaving the envelope detached from the core. In addition, both the
radius and orbital eccentricity are overestimated to obtain the desired
result.

In the previous sections, we showed that tidal shear energy dissipation can
play an important role in describing the heartbeat-star phenomenon. In this
section we present results from applying the model to MACHO 80.7443.1718.

###  4.1 Brightness variations from tidal heating

We computed a TIDES-nvv model for MACHO 80.7443.1718 using as input the values
that were obtained by Jayasinghe et al. (2021) from their analysis of the
observational data. Table 5 lists these input parameters. The projected
rotation velocity deduced from the widths of the photospheric absorption lines
is 174 km/s which, with an assumed orbital inclination of 44∘, implies an
equatorial rotation velocity Vrot=250 km/s. This yields $\beta_{0}$=1.95. The
range in Vrot values implied by the measurement uncertainties in V$\sin$(i)
quoted by Jayasinghe et al. (2021) corresponds to a range
$\beta_{0}$=1.57-2.34. Hence, the primary star of the system is in
supersynchronous rotation for all possible Vrot values. We adopt
$\beta_{0}$=1.95 for the calculation, which means that at apastron it is even
more highly supersynchronous, $\beta_{ap}$=18.2. We used the same
computational parameters as in Table 1 except for the number of layers that we
set to $N_{r}$=5 since the largest tidal shear energy dissipation rates occur
in the surface layers. The model was evolved out to 100 orbital periods (9
years), with output retained at 5, 10, 20 and 50 orbital periods.

The resulting total tidal energy dissipation rate $\dot{E}_{tot}$ at
periastron and apastron is, respectively, 3.7 1038 ergs/s and 6.6 1036 ergs/s.
The normalized light curve, obtained by adding the phase-dependent
$\dot{E}_{tot}$ to the stellar luminosity L∗=1.57 1039 ergs/s cited in
Jayasinghe et al. (2021) and dividing by this number, is illustrated in Fig.
4. It shows a 23% increase at maximum, with maximum occurring at orbital phase
$\sim$0.06 after periastron passage. The brightness increase is not
symmetrical around periastron but initiates at $\varphi\sim$0.9 and ends at
$\varphi\sim$0.3. These results are all consistent with the observations.

We emphasize that our results follow directly from a calculation performed
from first principles, with no fitting to the data, no prior assumptions and
no adjustable parameters other than the conversion rate of kinetic energy into
heat ($\lambda$). However, several comments are in order. The n-layer,
variable-viscosity version of our model provides only the tidally induced
perturbations in the azimuthal direction due to computational
challenges111111With n-layers, the absence of a lower boundary for each layer
allows vertical displacements of the volume elements that exceed their
vertical extension; that is, they can become detached from each other, a
condition that halts the numerical integration. A full SPH computation is
required to keep track of the 3D location of each volume element and allow for
buoyancy effects. Thus, it is likely to yield an underestimate of the total
energy dissipation rate. On the other hand, the value of $\lambda$ we have
used is the upper limit of this parameter and is certainly an overestimate.
Adding in the energy dissipation rate associated with the radial perturbation
components might counterbalance the use of a smaller, more realistic $\lambda$
value.

A final comment is in order regarding the average tidal energy dissipation
rate that we find in our model which, at first sight, might seem very large in
the context of the tidal evolution of the system. Using Eq. (9) of
Koenigsberger & Moreno (2016) and the energy dissipation rate average over an
orbit, we find $\dot{P}/P\simeq$1.4$\times$10-3 s/yr. Kołaczek-Szymański et
al. (2022) reported variations in the light curve that they interpret in terms
of a decreasing orbital period at a rate of 11 s/yr, which would imply an
energy loss from the orbit that is orders of magnitude larger than that
associated with the tidal shear energy dissipation rate that we have derived.

Figure 4: Predicted MACHO 80.7443.1718 light curve from the tidal shear energy
dissipation calculation added to the stellar luminosity given by J21, L∗=1.63
1039 ergs/s, and normalized to this value. The vertical dotted line indicates
the time of periastron and the horizontal line indicates the normalized
luminosity level. Tidal shear contributes $\sim$23% of the luminosity peak
around the time of periastron.  Figure 5: Equatorial radius of MACHO
80.7443.1718 at periastron as a function of the longitude obtained with the
TIDES-1 models using $\nu$=0.028 R${}^{2}_{\odot}$/d (green) and $\nu$=2
R${}^{2}_{\odot}$/d (red), and the TIDES-1vv with $\lambda$=0.01 (blue) and
$\lambda$=1 (magenta). The sub-binary longitude is at $\varphi$=0∘. Each dot
in the plot correspond to a grid point. The smaller viscosity models display
28 wave-like oscillations in the radius over the 2$\pi$ radians of the stellar
circumference.

Figure 6: Velocity in the radial direction (left) and azimuthal direction
(right) along the equator of the MACHO 80.7443.1718 surface at periastron as a
function of the longitude. The colors indicate results from the calculations
using the code TIDES-1vv with $\lambda$=0.01 (blue); and with $\lambda$=1
(magenta). The sub-binary longitude is at $\varphi$=0∘. Each dot in the plot
correspond to a grid point. The smaller viscosity models ($\lambda$=0.1)
display wave-like velocity structures along the stellar circumference.

###  4.2 Radius and perturbed radial velocity component

In order to explore the effects of the radial component of the tidal
perturbation we performed two tests with a one-layer model based on the TIDES-
nvv algorithm but that allows radial motions. As with the original version, in
this new version (TIDES-1vv) the surface layer remains attached to the
underlying rotating core at all times and is coupled to it through the viscous
force. The first test was performed with $\lambda$=1, as was used with the
n-layer calculation and the second with $\lambda$=0.01. For comparison, we
also computed two models using the original one-layer model discussed in
Moreno et al. (2005) and Moreno et al. (2011). All input parameters are the
same as for the TIDES-1vv calculation except that in this case a fixed
viscosity value is needed as input. We explored values $\nu$=0.028
R${}_{\odot}^{2}$/d (1.6 1015 cm2/s) and $\nu$=2 R${}_{\odot}^{2}$/d (2.3 1017
cm2/s) which are on the order of magnitude of those that were computed by
TIDES-nvv with $\lambda$=1. All the calculations were performed with 500
partitions in the azimuthal direction.

We find that the primary tidal bulge at periastron extends no more than
$\sim$0.13 R⊙ above the equilibrium radius. The four calculations yield a
similar average shape (Fig. 5). The differences can be traced to the viscosity
values, with smaller viscosities allowing for the appearance of high frequency
oscillations in both radius and velocity.

The size of the main tidal bulge found in our calculation can be compared with
the estimate (see, for example, Zahn (2008):

$\frac{\delta
R}{R_{1}}\simeq\frac{m_{1}}{m_{2}}\left(\frac{R_{1}}{a}\right)^{3},$ (5)

where $R_{1}$ is the equilibrium radius, $m_{1}$ and $m_{2}$ the masses of the
primary star and the perturber, respectively, $a$ is the semimajor axis, and
$\delta R=R_{bulge}-R_{1}$. This relation gives $\delta R\simeq$0.17 R⊙, which
is comparable to our maximum value 0.13 R⊙.

In general, the tidal velocity field has components in the three spatial
directions. Our one-layer computations calculate the values in the radial and
azimuthal directions. The two TIDES-1vv tests yield the behavior of the radial
component that is shown in Fig. 6 (left). The $\lambda$=0.01 calculation
displays high frequency oscillations with maximum amplitudes $\sim$2 km/s at
longitudes 300∘ \- 360∘. Thus, our model indicates that the rise to the tidal
bulge is through a series of oscillations of increasing amplitude. The
azimuthal velocity component also displays high frequency oscillations which,
however, have a small amplitude compared to the large-scale equilibrium tide
component which reaches speeds as high as $\sim$100 km/s as it rises into the
primary tidal bulge (Fig. 6 right).

The most significant difference between the one-layer calculations is the
presence of high frequency oscillations in the models having smaller viscosity
values. This suggests that although large $\lambda$ values are needed for our
results to be consistent with the observed maximum brightness, the value of
this parameter needs to decrease at the surface in order to allow for TEOs, a
reasonable expectation given the significantly lower density of the surface
compared to inner layers, as well as its optically thinner properties allowing
energy to escape. Allowing for a variable $\lambda$ parameter is, however,
beyond the scope of the current paper.121212We note that the percentage
contribution from each layer of radius $r$/R⊙=[16.6,18.0,19.4,20.9,22.3] to
the total energy dissipation rate at periastron is respectively [$<$0.01, 0.2,
0.55, 0.2, 0.05]. Hence, a very small $\lambda$ value in the surface layer
will not affect the conclusion that tidal shear energy dissipation can explain
the brightness variations of MACHO 80.7443.1718 and other heartbeat stars.

Figure 7: Line profiles over an orbital cycle computed with the one-layer
model of Moreno et al. (2005) for MACHO 80.7443.1718 assuming $\beta_{0}$=1.95
and the masses, radius, and orbital elements given by Jayasinghe et al.
(2021). The unperturbed and rotationally broadened profiles (red dots) are
superposed on the predicted profiles (black) which show significant
distortions only in the orbital phase between just before periastron
($\varphi\sim$0.97) until $\sim$0.3 in phase after periastron, consistent with
the duration of the luminosity increase.

###  4.3 Photospheric line profile variations

The tidal flows and radial motions excited by the tidal perturbation alter the
smoothly varying rotation velocity across the stellar disk, projected onto the
line of sight to the observer. Hence, instead of the typical rotationally
broadened photospheric absorption line profile that would be observed if the
surface were stationary in the rest frame, it appears distorted and varies
over time as different portions of the stellar surface come into view. Moreno
et al. (2005) incorporated the calculation of line profiles into the TIDES-1
calculation. The tidally perturbed velocities on the visible portion of the
stellar disk are projected onto the line of sight to the observer. The
corresponding Doppler shifts are applied to the local (on each surface
element) line profile and all line profiles are added to produce the
photospheric line profile that is predicted for each particular orbital phase.
The method has successfully modeled the line profile variations that are
observed in the eccentric B-type binary Spica (Harrington et al., 2009; Palate
et al., 2013a; Harrington et al., 2016).

Fig. 7 is a montage of line profiles predicted by the TIDES-1 calculation
using the input parameters for MACHO 80.7443.1718 in Table 5. They are stacked
from bottom to top in order of increasing orbital phase and show that the
strongest perturbations are expected to be observed just before periastron and
until approximately 0.3 in phase after periastron. It would be interesting to
determine whether similar variations are indeed present in lines that are not
affected by a stellar wind, such as He 1 $\lambda$ 4921 shown in Jayasinghe et
al. (2021) but, unfortunately, the high quality spectra presented by these
authors only cover the orbital phase range 0.27-0.81.

Figure 8: Distribution of the tidal shear energy dissipation in MACHO
80.7443.1718 from the TIDES-nvv model at periastron (left) and apastron
(right). Top: Maps of the surface, with the sub-binary longitude at
$\varphi$=0∘ and 360∘ (the rims of the maps),and the poles at the top and
bottom. Dark brown represents maximum energy dissipation rate and minimum is
color coded white. Bottom: Plots of the corresponding $\dot{E}$ per unit
volume (in units of 1035 ergs/s/unit volume) at various latitudes between the
equator and the pole. Each curve is labeled with its corresponding colatitude
angle. Each latitude is shifted vertically by 0.03 units for display purposes.

###  4.4 Non-isotropic energy dissipation structure

The distribution of the tidal shear energy dissipation rate over the surface
is non-isotropic, and changes between periastron and apastron as shown in Fig.
8. At periastron, maximum energy is concentrated around the sub-binary
longitude and extends to colatitudes as high as $\sim$60∘. Three additional
locations of large energy dissipation rates can be seen at approximately 90∘,
180∘ and 270∘ from the sub-binary longitude. At apastron, maxima occur near
90∘ and 270∘ at the equator, however at higher latitudes maxima are spread out
in long arch-like structures spanning broad ranges in longitude.

###  4.5 Conclusions of this section

Our n-layer calculation of the tidal shear energy dissipation rate yields a
light curve with a $\sim$23% maximum which occurs slightly after periastron
passage, as is observed. These results follow directly from a calculation
performed from first principles with no fitting to the data and no adjustable
parameters other than the conversion rate of kinetic energy into heat
($\lambda$). The only basic assumption in our treatment is that relative
motions excited by the tidal force gives rise to a turbulent viscosity.

We find that the maximum radial deformation of the star at periastron is
$\leq$0.2 R⊙, contrary to the several R⊙ that was found by MacLeod & Loeb
(2023), and that the radial velocity of the material rising into the primary
bulge is subsonic, $\sim$2 km/s, contrary to the 20 km/s found by these
authors. Our results indicate that it is the horizontal velocity component
that reaches supersonic velocities, attaining speeds as high as $\sim$50 km/s.

The light curve does not predict the minimum that is observed at periastron.
This minimum could result from a physical eclipse as suggested by MacLeod &
Loeb (2023), but according to the parameters derived by Jayasinghe et al.
(2021) such an eclipse would be at best a grazing eclipse. The tidal shear
energy dissipation models discussed in Section 3 do predict a minimum at
periastron but only in stars that are rotating close to synchronicity at this
orbital phase. MACHO 80.7443.1718 is rotating at approximately twice the
synchronous rate. An alternative scenario for a light curve minimum is one in
which a puff of dense material is ejected at periastron. This material would
dim the star until it has expanded enough to become optically thin, at which
time it would become the source of emission lines such as H$\alpha$ and
H$\beta$. MacLeod & Loeb (2023) also suggested that the activity at periastron
could lead to material leaving the star.

Jayasinghe et al. (2021) detected spectral variations in MACHO 80.7443.1718
that suggest that, indeed, enhanced mass-loss episodes occur. Hence, the
prominent dip in the light curve, not predicted by our tidal shear energy
dissipation model, could be caused by the ejected material. We speculate that
the mechanism for mass ejection may be found in the large velocity gradients
around periastron combined with the rapid rise in the stellar bulge.
Furthermore, similar proceses are likely also present in the companion, whose
contribution to the observed phenomenon should in the future be assessed.

##  5 Summary and conclusions

We examine the role of the often neglected process of tidal shear energy
dissipation $\dot{E}$ in shaping the observational properties of stars in
which periastron brightening events occur. Our model consists of a direct
calculation from first principles of the motions induced on a rotating star by
a companion in an eccentric orbit. The value of $\dot{E}$ is computed using a
non-isotropic and time-dependent viscosity. We explore the impact of stellar
rotation on $\dot{E}$ in a test binary consisting of a 10 M⊙ primary that is
perturbed by a 7 M⊙ companion with orbital periods in the range 6 d-31 d and
eccentricities $e$ in the range 0$-$0.7.

The model grid includes rotation rates at periastron that are subsynchronous
($\beta_{0}<$1), synchronous ($\beta_{0}$=1), and supersynchronous
($\beta_{0}>$1), where $\beta_{0}$=$\omega_{0}/\Omega_{0}$. Here, $\omega_{0}$
is the rotation angular velocity of the inner, rigidly rotating region that
includes the stellar core and $\Omega_{0}$ is the orbital angular velocity at
periastron.

We find that the shear energy dissipation rate increases by factors $\sim$10-6
-10-3 at periastron, consistent with the majority of observed heartbeat stars,
and we find the following trends in $\dot{E}$: a) maximum $\dot{E}$ occurs at
orbital phases slightly after periastron, not at periastron; b) the timescale
over which $\dot{E}$ increases and then decreases is shorter in the
subsynchronous models than in the supersynchronous models; c) for fixed values
of orbital period and eccentricity, $\dot{E}$ grows with increasing departure
from $\beta_{0}\sim$1, with the largest values appearing for $\beta_{0}>>$1;
d) for fixed values of $\beta_{0}$ and periastron distance, models with larger
eccentricity have a stronger periastron effect; e) models with
$|\beta_{0}-1|\simeq$1 have the lowest dissipation rates, present a dip at
periastron, and display the relatively prominent high frequency variations
that may be identified with tidally excited oscillations.

Pseudosyncronism is defined as the minimum energy rotation state in an
eccentric binary. In this state, $\beta_{ps}$=$\omega_{ps}/\Omega_{ave}$,
where $\omega_{ps}$ is the rotation angular velocity and $\Omega_{ave}$ is the
average orbital angular velocity. We confirm the classical result that
0.8$<\beta_{ps}<$1.0 for the cases with e$<$0.3. However, in the $e$=0.5 and
0.7 cases, we find 1.0$<\beta_{ps}<$1.15, contrary to the results of Hut
(1981) who found an approximately constant $\beta_{ps}\simeq$0.82 for
eccentricities $e>$0.2.

In order to compare with observational resuls, we apply the model to the
massive (M1+M2=34.5 + 15.7 M⊙, $P$=33 d, $e$=0.5) heartbeat star MACHO
80.7443.1718 which displays an unusually large ($\sim$23%) brightness increase
at periastron. Our model directly predicts a $\sim$23% increase in tidal shear
energy dissipation rates at periastron. It thus accounts in a straightforward
manner for the extraordinarily large brightness increase that is observed.
However, this assumes that all of the dissipated energy contributes to the
luminosity. If part of the energy remains within the star, the radiated energy
would be smaller. On the other hand, however, the energy that may be trapped
in sub-photospheric layers could trigger a temporary increase in the stellar
radius which would then lead to even larger tidal perturbations. Such
processes would be interesting to analyze in the future.

Based on the results of our grid of models discussed in the first part of the
paper, we are led to conclude that MACHO 80.7443.1718’s unusually strong
periastron effect is due to its highly super-synchronous rotation rate
($\beta_{0}\sim$1.95), a state that is not expected to be very common in
binary systems if they have had time to pseudosynchronize. MACHO 80.7443.1718
belongs to the Lucke-Hodge 58 association which has an estimated age of only a
couple of million years (Garmany et al., 1994), which could explain why it
still is in such rapid rotation. It is also interesting to note that its
$\beta_{0}$ value exceeds the limit $\beta_{crit}$=18/11 at which the
eccentricity tends to increase over time instead of decrease (Zahn, 2008).

Our model for MACHO 80.7443.1718 does not predict a dip in the brightness,
whereas it’s observed light curve does show such a dip. MacLeod & Loeb (2023)
suggested that the dip could be due to obscuration caused by ejected material.
It is tempting to speculate that the steep surface angular velocity gradients
and fast bulge radius increase that take place around periastron could trigger
brief mass-shedding events. In this respect, a determination of the mass-loss
rates at different orbital phases would be most useful.

The effect of tidal shear energy dissipation is often neglected when
considering the structure of binary stars and their observational properties.
This neglect is due to gaps in our understanding of the dissipation processes
that are active and the assumption that viscous shear energy dissipation is
negligible in stars.131313This process, however, plays a major role in
celestial objects such as Jupiter’s moons Io and Europa and most likely play a
role in close-in exoplanets. However, if it is non-negligible, as we assume in
this paper, it can have a significant effect on the stellar properties which,
in eccentric systems, become strongly orbital phase-dependent.

###### Acknowledgements.

We thank Werner Schmutz for a critical reading of this paper and the anonymous
referee who provided thoughtful comments and suggestions. GK is grateful to
the Astronomy Department of the Indiana University for hosting a sabbatical
visit during which the research in this paper was concluded. This work was
supported by UNAM-PAPIIT 105723.

## References

  * Alcock et al. (1997) Alcock, C., Allsman, R. A., Alves, D., et al. 1997, ApJ, 486, 697
  * Astoul et al. (2021) Astoul, A., Park, J., Mathis, S., Baruteau, C., & Gallet, F. 2021, A&A, 647, A144
  * Barker & Astoul (2021) Barker, A. J. & Astoul, A. A. V. 2021, MNRAS, 506, L69
  * Baruteau & Rieutord (2013) Baruteau, C. & Rieutord, M. 2013, Journal of Fluid Mechanics, 719, 47
  * Caleo et al. (2016) Caleo, A., Balbus, S. A., & Tognelli, E. 2016, MNRAS, 460, 338
  * Dewberry & Wu (2024) Dewberry, J. W. & Wu, S. C. 2024, MNRAS, 527, 2288
  * Efroimsky (2018) Efroimsky, M. 2018, Icarus, 300, 223
  * Eggleton et al. (1998) Eggleton, P. P., Kiseleva, L. G., & Hut, P. 1998, ApJ, 499, 853
  * Elliott et al. (2000) Elliott, J. R., Miesch, M. S., & Toomre, J. 2000, ApJ, 533, 546
  * Estrella-Trujillo et al. (2023) Estrella-Trujillo, D., Arthur, S. J., Koenigsberger, G., & Moreno, E. 2023, A&A, 670, A44
  * Fellay et al. (2024) Fellay, L., Dupret, M. A., & Rosu, S. 2024, arXiv e-prints, arXiv:2401.02573
  * Fuller (2017) Fuller, J. 2017, MNRAS, 472, 1538
  * Garaud & Kulenthirarajah (2016) Garaud, P. & Kulenthirarajah, L. 2016, ApJ, 821, 49
  * Garmany et al. (1994) Garmany, C. D., Massey, P., & Parker, J. W. 1994, AJ, 108, 1256
  * Goldreich & Nicholson (1989) Goldreich, P. & Nicholson, P. D. 1989, ApJ, 342, 1079
  * Guo (2020) Guo, Z. 2020, in Stars and their Variability Observed from Space, ed. C. Neiner, W. W. Weiss, D. Baade, R. E. Griffin, C. C. Lovekin, & A. F. J. Moffat, 203–207
  * Guo et al. (2020) Guo, Z., Shporer, A., Hambleton, K., & Isaacson, H. 2020, ApJ, 888, 95
  * Hambleton et al. (2013) Hambleton, K., Degroote, P., Conroy, K., et al. 2013, in EAS Publications Series, Vol. 64, EAS Publications Series, ed. K. Pavlovski, A. Tkachenko, & G. Torres, 285–294
  * Harrington et al. (2009) Harrington, D., Koenigsberger, G., Moreno, E., & Kuhn, J. 2009, ApJ, 704, 813
  * Harrington et al. (2016) Harrington, D., Koenigsberger, G., Olguín, E., et al. 2016, A&A, 590, A54
  * Hirschi & Maeder (2010) Hirschi, R. & Maeder, A. 2010, A&A, 519, A16
  * Hut (1981) Hut, P. 1981, A&A, 99, 126
  * Hutchings (1979) Hutchings, J. B. 1979, in Mass Loss and Evolution of O-Type Stars, ed. P. S. Conti & C. W. H. De Loore, Vol. 83, 3–20
  * Jayasinghe et al. (2021) Jayasinghe, T., Kochanek, C. S., Strader, J., et al. 2021, MNRAS, 506, 4083
  * Jayasinghe et al. (2020) Jayasinghe, T., Stanek, K. Z., Kochanek, C. S., et al. 2020, MNRAS, 491, 13
  * Jayasinghe et al. (2019) Jayasinghe, T., Stanek, K. Z., Kochanek, C. S., et al. 2019, MNRAS, 489, 4705
  * Käpylä et al. (2020) Käpylä, P. J., Rheinhardt, M., Brandenburg, A., & Käpylä, M. J. 2020, A&A, 636, A93
  * Kirk et al. (2016) Kirk, B., Conroy, K., Prša, A., et al. 2016, AJ, 151, 68
  * Koenigsberger & Moreno (2016) Koenigsberger, G. & Moreno, E. 2016, Rev. Mexicana Astron. Astrofis., 52, 113
  * Koenigsberger et al. (2021) Koenigsberger, G., Moreno, E., & Langer, N. 2021, A&A, 653, A127
  * Kołaczek-Szymański et al. (2021) Kołaczek-Szymański, P. A., Pigulski, A., Michalska, G., Moździerski, D., & Różański, T. 2021, A&A, 647, A12
  * Kołaczek-Szymański et al. (2022) Kołaczek-Szymański, P. A., Pigulski, A., Wrona, M., Ratajczak, M., & Udalski, A. 2022, A&A, 659, A47
  * Kumar et al. (1995) Kumar, P., Ao, C. O., & Quataert, E. J. 1995, ApJ, 449, 294
  * Lai (1997) Lai, D. 1997, ApJ, 490, 847
  * Landau & Lifshitz (1987) Landau, L. D. & Lifshitz, E. M. 1987, Fluid Mechanics
  * Le Bars et al. (2022) Le Bars, M., Barik, A., Burmann, F., et al. 2022, Surveys in Geophysics, 43, 229
  * Lee (2017) Lee, U. 2017, MNRAS, 468, 1864
  * Lurie et al. (2017) Lurie, J. C., Vyhmeister, K., Hawley, S. L., et al. 2017, AJ, 154, 250
  * MacLeod & Loeb (2023) MacLeod, M. & Loeb, A. 2023, Nature Astronomy, 7, 1218
  * Maeder & Meynet (1996) Maeder, A. & Meynet, G. 1996, A&A, 313, 140
  * Maeder & Meynet (2000) Maeder, A. & Meynet, G. 2000, ARA&A, 38, 143
  * Mathis et al. (2016) Mathis, S., Auclair-Desrotour, P., Guenel, M., Gallet, F., & Le Poncin-Lafitte, C. 2016, A&A, 592, A33
  * Mathis et al. (2004) Mathis, S., Palacios, A., & Zahn, J. P. 2004, A&A, 425, 243
  * McMillan et al. (1987) McMillan, S. L. W., McDermott, P. N., & Taam, R. E. 1987, ApJ, 318, 261
  * Moreno & Koenigsberger (1999) Moreno, E. & Koenigsberger, G. 1999, Rev. Mexicana Astron. Astrofis., 35, 157
  * Moreno et al. (2011) Moreno, E., Koenigsberger, G., & Harrington, D. M. 2011, A&A, 528, 48
  * Moreno et al. (2005) Moreno, E., Koenigsberger, G., & Toledano, O. 2005, A&A, 437, 641
  * Ogilvie (2014) Ogilvie, G. I. 2014, ARA&A, 52, 171
  * Orosz & Hauschildt (2000) Orosz, J. A. & Hauschildt, P. H. 2000, A&A, 364, 265
  * Pablo et al. (2017) Pablo, H., Richardson, N. D., Fuller, J., et al. 2017, MNRAS, 467, 2494
  * Palate et al. (2013a) Palate, M., Koenigsberger, G., Rauw, G., Harrington, D., & Moreno, E. 2013a, A&A, 556, A49
  * Palate et al. (2013b) Palate, M., Rauw, G., Koenigsberger, G., & Moreno, E. 2013b, A&A, 552, A39
  * Penev et al. (2007) Penev, K., Sasselov, D., Robinson, F., & Demarque, P. 2007, ApJ, 655, 1166
  * Pérez-de-Tejada (1999) Pérez-de-Tejada, H. 1999, ApJ, 525, L65
  * Pigulski (2018) Pigulski, A. 2018, in XXXVIII Polish Astronomical Society Meeting, ed. A. Ró&żańska, Vol. 7, 151–156
  * Press & Teukolsky (1977) Press, W. H. & Teukolsky, S. A. 1977, ApJ, 213, 183
  * Richard & Zahn (1999) Richard, D. & Zahn, J.-P. 1999, A&A, 347, 734
  * Richardson et al. (2018) Richardson, N. D., Pablo, H., Sterken, C., et al. 2018, MNRAS, 475, 5417
  * Savonije et al. (1995) Savonije, G. J., Papaloizou, J. C. B., & Alberts, F. 1995, MNRAS, 277, 471
  * Scharlemann (1981) Scharlemann, E. T. 1981, ApJ, 246, 292
  * Schenk et al. (2001) Schenk, A. K., Arras, P., Flanagan, É. É., Teukolsky, S. A., & Wasserman, I. 2001, Phys. Rev. D, 65, 024001
  * Shi et al. (2006) Shi, X., Wen, D., Bao, X., Li, C., & Huang, Y. 2006, Science China Earth Sciences, 49, 492
  * Sterken & Breysacher (1997) Sterken, C. & Breysacher, J. 1997, A&A, 328, 269
  * Sun et al. (2023) Sun, M., Townsend, R. H. D., & Guo, Z. 2023, ApJ, 945, 43
  * Symon (1971) Symon, K. R. 1971, Mechanics (Addison-Wesley)
  * Talon et al. (1997) Talon, S., Zahn, J. P., Maeder, A., & Meynet, G. 1997, A&A, 322, 209
  * Terquem (2021) Terquem, C. 2021, MNRAS, 503, 5789
  * Thompson et al. (2012) Thompson, S. E., Everett, M., Mullally, F., et al. 2012, ApJ, 753, 86
  * Tkachenko et al. (2016) Tkachenko, A., Matthews, J. M., Aerts, C., et al. 2016, MNRAS, 458, 1964
  * Toledano et al. (2007) Toledano, O., Moreno, E., Koenigsberger, G., Detmers, R., & Langer, N. 2007, A&A, 461, 1057
  * Townsend & Sun (2023) Townsend, R. H. D. & Sun, M. 2023, ApJ, 953, 48
  * van Genderen & Sterken (2007) van Genderen, A. M. & Sterken, C. 2007, Information Bulletin on Variable Stars, 5782, 1
  * Zahn (1975) Zahn, J. P. 1975, A&A, 41, 329
  * Zahn (1977) Zahn, J. P. 1977, A&A, 57, 383
  * Zahn (1992) Zahn, J. P. 1992, A&A, 265, 115
  * Zahn (2008) Zahn, J. P. 2008, in EAS Publications Series, Vol. 29, EAS Publications Series, ed. M. J. Goupil & J. P. Zahn, 67–90

##  Appendix A Extremum value of the angular velocity at each timestep

The eccentric binary models all display a variability pattern in the tidal
perturbations on the orbital timescale, the largest amplitudes occurring
around the time of periastron. We quantify the time-dependent perturbations at
each computational timestep by finding the largest deviation of the angular
velocity from the underlying rigidly rotating core at each computational time
step. We define this quantity as max(—$\omega(r,\theta,\varphi,t|$), where
$\omega$ is the instantaneous angular velocity of the volume element located
at ($r,\theta,\varphi$) at time $t$ as measured in the reference frame that
rotates with the primary star. Once this quantity is located at each timestep
of the computation, it is assigned a positive or negative sign depending on
whether this maximum velocity is in the direction of the orbital motion or in
the opposite direction. The first case occurs mostly in the subsynchronous
models ($\beta_{0}<$1) where the tidal force acts to increase the velocity
toward corotation. The second case occurs in the supersynchronous models where
the tidal force acts to slow the velocity down. This general behavior is
illustrated in Fig. 9 for $e$=0.1 and 0.3 models and Fig. 10 for $e$=0.5.

Around the time of periastron, the perturbation becomes strongest and the
angular velocity changes rapidly creating either a maximum (subsynchronous
systems) or a minimum peak (highly supersynchronous systems). Systems with
rotation speeds that are not very far from corotation at periastron tend to
present high frequency oscillations throughout the orbital cycle and a variety
of behaviors at periastron.

Some of the models show a maximum peak at periastron and negative values at
other phases (for example, $e$=0.3, $\beta_{0}$=0.60 and $e$=0.5,
$\beta_{0}$=0.40) but in general have either positive or negative values
throughout the orbit.

Figure 9: Maximum perturbation of the surface angular velocity at the equator
plotted as a function of time for $e$=0.1 (left) and $e$=0.3 (right). The
velocity is measured in the rest reference frame of the rotating star.
Positive values indicate motion in the direction of the orbital motion,
negative values to motion opposite to the orbital motion. Vertical lines
indicate times of periastron. Each curve is labeled with its value of
$\beta_{0}$. The smallest tidal perturbations for $e$=0.1 occur for
$\beta_{0}$=0.88. For $e$=0.3, the $\beta_{0}$=0.80, 0.88, 1.00 and 1.05
curves show low amplitude oscillations and are very similar. The smallest
perturbations are observed for $\beta_{0}$=0.95. Observational perspective:
Surface layer is rotating faster than the inner core for $\beta_{0}<$1 and
rotating slower in $\beta_{0}>$1 models. Also, these stars are likely to
display variations in the shape of the photospheric line profiles as discussed
in Moreno et al. (2005) and Harrington et al. (2009, 2016). TEOs are evident
mostly for $\beta_{0}\sim$1 models.

Figure 10: Maximum perturbation of the surface angular velocity at the equator
for $e$=0.5 plotted as a function of time. The velocity is measured in the
reference frame at rest in the rotating star. Positive values indicate motion
in the direction of the orbital motion, negative values motion opposite to the
orbital motion. Vertical lines indicate times of periastron. Left: All models
with $\beta_{0}$ increasing from top to bottom. Models close to
pseudosynchronism are not labeled for clarity in the figure. Right: Models
with 0.60$\leq\beta_{0}\leq$1.00 are plotted in blue and cyan and those with
1.05$\leq\beta_{0}\leq$1.20 in red.  Table 6: Turbulent viscosity ten layer
models.

Model | $\beta_{0}$ | Vrot | $<\dot{E}_{tot}>$ | $\log\left<\nu/R_{\odot}^{2}/d\right>$
---|---|---|---|---
|  |  |  | Periastron | Apastron
| k=2 | k=4 | k=6 | k=8 | k=10 | k=2 | k=4 | k=6 | k=8 | k=10
$e$=0.1 P=6 d
201 | 0.20 | 11 | 7.61e-1 | -3.20 | -2.86 | -2.72 | -2.29 | -1.17 | -3.20 | -2.89 | -2.72 | -2.78 | -1.96
202 | 0.40 | 22 | 2.26e-1 | -3.60 | -3.03 | -2.86 | -2.43 | -1.36 | -3.59 | -3.03 | -2.84 | -2.97 | -1.94
203 | 0.60 | 34 | 4.78e-2 | -3.81 | -3.33 | -3.02 | -2.62 | -1.50 | -4.18 | -3.32 | -3.02 | -3.12 | -2.28
204 | 0.80 | 45 | 4.19e-3 | -3.80 | -3.95 | -3.64 | -2.92 | -1.81 | -4.59 | -4.20 | -3.64 | -3.37 | -1.97
205 | 0.88 | 49 | 2.48e-3 | -3.87 | -4.35 | -3.66 | -3.12 | -1.91 | -4.76 | -4.50 | -3.81 | -3.21 | -1.87
206 | 0.95 | 53 | 4.19e-3 | -3.84 | -4.04 | -3.36 | -3.33 | -2.10 | -4.74 | -4.06 | -3.36 | -3.09 | -1.83
207 | 1.00 | 56 | 8.78e-3 | -3.85 | -3.78 | -3.21 | -3.26 | -2.21 | -4.85 | -3.78 | -3.21 | -2.99 | -1.78
208 | 1.05 | 59 | 1.48e-2 | -3.86 | -3.56 | -3.12 | -3.17 | -2.35 | -4.61 | -3.55 | -3.12 | -2.93 | -1.73
209 | 1.20 | 67 | 6.32e-2 | -3.86 | -3.19 | -2.94 | -2.96 | -1.96 | -4.10 | -3.19 | -2.94 | -2.79 | .1.61
210 | 1.40 | 79 | 2.69e-1 | -3.44 | -2.95 | -2.78 | -2.67 | -1.61 | -3.43 | -2.95 | -2.78 | -2.64 | -1.47
211 | 1.60 | 90 | 8.33e-1 | -3.12 | -2.81 | -2.68 | -2.49 | -1.40 | -3.11 | -2.81 | -2.68 | -2.52 | -1.36
212 | 1.88 | 105 | 3.01 | -2.88 | -2.66 | -2.56 | -2.29 | -1.18 | -2.88 | -2.66 | -2.57 | -2.36 | -1.18
213 | 1.95 | 109 | 4.20 | -2.83 | -2.63 | -2.54 | -2.25 | -1.12 | -2.83 | -2.63 | -2.54 | -2.32 | -1.13
$e$=0.3 P=8.73 d
301 | 0.20 | 13 | 6.61e-1 | -3.55 | -3.02 | -2.83 | -2.14 | -1.09 | -3.54 | -3.02 | -2.85 | -2.77 | -1.61
302 | 0.40 | 27 | 1.81e-1 | -4.16 | -3.42 | -3.12 | -2.29 | -1.27 | -4.15 | -3.42 | -3.12 | -3.14 | -1.70
303 | 0.60 | 40 | 5.28e-2 | -4.22 | -3.85 | -3.35 | -2.46 | -1.35 | -4.40 | -4.14 | -3.73 | -3.24 | -1.63
304 | 0.80 | 53 | 2.47e-2 | -4.29 | -3.88 | -3.26 | -2.66 | -1.51 | -4.40 | -3.88 | -3.26 | -3.29 | -1.49
305 | 0.88 | 59 | 1.58e-2 | -4.46 | -3.88 | -3.30 | -2.90 | -1.75 | -4.54 | -3.86 | -3.30 | -3.23 | -1.40
306 | 0.95 | 63 | 4.21e-4 | -4.37 | -4.42 | -4.00 | -3.92 | -2.40 | -4.55 | -4.21 | -3.73 | -3.81 | -1.30
307 | 1.00 | 67 | 2.34e-2 | -4.47 | -3.61 | -3.18 | -3.32 | -2.02 | -4.44 | -3.60 | -3.18 | -3.03 | -1.32
308 | 1.05 | 70 | 2.66e-2 | -4.48 | -3.56 | -3.16 | -3.43 | -2.13 | -4.27 | -3.55 | -3.16 | -2.98 | -1.29
309 | 1.20 | 80 | 5.80e-2 | -4.25 | -3.32 | -3.03 | -2.97 | -1.86 | -4.26 | -3.32 | -3.03 | -2.96 | -1.22
310 | 1.40 | 93 | 2.02e-1 | -3.74 | -3.06 | -2.86 | -2.62 | -1.53 | -3.72 | -3.06 | -2.86 | -2.84 | -1.14
311 | 1.60 | 107 | 7.18e-1 | -3.21 | -2.85 | -2.70 | -2.41 | -1.31 | -3.20 | -2.85 | -2.70 | -2.69 | -1.11
312 | 1.88 | 125 | 4.36 | -2.80 | -2.60 | -2.49 | -2.21 | -1.09 | -2.80 | -2.60 | -2.49 | -2.52 | -1.14
313 | 1.95 | 130 | 6.45 | -2.75 | -2.56 | -2.46 | -2.17 | -1.03 | -2.75 | -2.56 | -2.46 | -2.44 | -1.10
$e$=0.0 P=6 d |  |
102 | 0.40 | 18 | 1.83e-1 | -3.23 | -3.04 | -2.96 | -2.74 | -1.62 |  |  |  |  |
103 | 0.60 | 27 | 3.79e-2 | -4.20 | -3.29 | -3.00 | -2.88 | -1.79 |  |  |  |  |
104 | 0.80 | 36 | 6.01e-3 | -4.09 | -3.53 | -3.34 | -3.22 | -2.10 |  |  |  |  |
106 | 0.95 | 43 | 2.09e-4 | -4.28 | -4.31 | -3.96 | -3.76 | -2.70 |  |  |  |  |
108 | 1.05 | 48 | 1.41e-3 | -4.65 | -4.43 | -3.93 | -3.80 | -2.70 |  |  |  |  |
109 | 1.20 | 55 | 4.13e-3 | -4.22 | -3.86 | -3.30 | -3.16 | -2.10 |  |  |  |  |
111 | 1.60 | 73 | 1.52e-1 | -3.67 | -3.04 | -2.84 | -2.73 | -1.61 |  |  |  |  |
113 | 1.95 | 89 | 9.27e-1 | -3.12 | -2.81 | -2.68 | -2.51 | -1.38 |  |  |  |  |

141414Vrot is the initial equatorial rotation velocity in units of km/s;
$\dot{E}_{tot}$ is the total energy dissipation rate after the numerical
simulation has evolved to 50 orbital periods in units of 1035 ergs/s; $\nu$ is
the viscosity as computed in the model using Eq. 2 with $\lambda$=1 in units
of R${}_{\odot}^{2}$/d. $k$ indicates the layer with $k$=1 the layer adjoining
the core and $k$=10 the surface layer, each layer having a depth of
$\Delta$R=0.329 R⊙. The header of each set of model grids indicates the
corresponding orbital period and eccentricity.