Photometric Analysis of the OGLE Heartbeat Stars

We present an analysis of 991 heartbeat stars (HBSs) from the OGLE Collection of Variable Stars (OCVS). The sample consists of 512 objects located toward the Galactic bulge (GB), 439 in the Large Magellanic Cloud (LMC) and 40 in the Small Magellanic Cloud (SMC). We model the $I$-band OGLE light curves using an analytical model of flux variations, reflecting tidal deformations between stars. We present distributions of the model parameters that include the eccentricity, orbital inclination, and argument of the periastron, but also the period-amplitude diagrams. On the Hertzsprung-Russell (HR) diagram, our HBS sample forms two separate groups of different evolutionary status. The first group of about 90 systems, with short orbital periods ($P\lesssim50$~days), consists of an early-type primary star lying on (or close to) the main sequence (MS). The second group of about 900 systems, with long orbital periods ($P\gtrsim100$~days), contains a red giant (RG). The position of RG HBSs on the period-luminosity diagram strongly indicates their binary nature. They appear to be a natural extension of confirmed binary systems that include the OGLE ellipsoidal and Long Secondary Period (LSP) variables. We also present a time-series analysis leading to detection of tidally-excited oscillations (TEOs). We identify such pulsations in about 5\% of stars in the sample with a total number of 78 different modes. This first relatively large homogeneous sample of TEOs allowed us to construct a diagram revealing the correlation between the TEO's orbital harmonic number and the eccentricity of the host binary system.


INTRODUCTION
Heartbeat stars (HBSs) are a subclass of ellipsoidal variable stars with eccentric orbits. The ellipsoidal variable stars are close binary systems whose components are deformed due to gravitational tidal forces. For a "classical" ellipsoidal variable, which typically has a circular orbit, the light curve is characterized by a sinusoidal shape with two maxima and two minima per orbital period. This is a result of observing different sides of the stellar surface, which is distorted into a rotational ellipsoid (Morris 1985).
If the orbit of the system is eccentric, the deformation strength depends on the orbital phase, and the strongest effect appears during the periastron passage. This variation of the tidal force affects the shape of the light curve, which starts to deviate from the typical sinusoidal shape. The response of the stellar surface and volume to the presence of varying tidal potential can be expressed as a sum of two effects (Zahn 1975).
The first one is an equilibrium tide, which refers to the instantaneous deformation of a star due to tides. It is the only component of the tidal response provided that the system is a circular and synchronized binary. The equilibrium tide is responsible for the presence of the "heartbeat" feature in the light curves of HBSs. This prominent characteristic resembles a single electrocardiogram pulse, hence the name of this group of variables.
The second type of stellar tidal response is the dynamical tide, which may manifest itself as tidally excited oscillations (TEOs; e.g., Zahn 1970, Kumar et al. 1995, Fuller 2017. The periodically changing tidal potential may act as a driving force and induce even naturally damped pulsations, which would not be visible in the absence of a nearby companion. Most of the TEOs observed in main-sequence (MS) stars are high-radial order gravity modes (Guo 2021). Kumar et al. (1995) provided a convenient analytical model of the flux variations driven by tidal interactions (we will refer to it as the K95 model). The K95 model allows for determination of the orbital parameters of the system based only on the shape of the single-passband light curve. For instance, this model was successfully applied to derive orbital parameters for a sample of HBSs discovered by Thompson et al. (2012). The K95 model accounts for the ellipsoidal variability only. It neglects other proximity effects, such as the irradiation/reflection effect and Doppler beaming/boosting. While the latter should not be pronounced in our sample of HBSs because of their long orbital periods and hence low radial velocities, the former effect may play a potentially significant role (see, e.g., Welsh et al. 2011, their Figure 7).
The prototype of the entire HBS class is KOI-54 (HD 187091), studied in detail by many authors, e.g., Welsh et al. (2011), Burkart et al. (2012), and Fuller & Lai (2012). For the first time, the HBS was reported as a separate class of variable stars in the work of Thompson et al. (2012), but such objects were known in the literature much earlier (e.g., Handler et al. 2002, Maceroni et al. 2009). A group of more than 100 ellipsoidal variables with an eccentric orbit was highlighted in the work of Soszyński et al. (2004). A subsample of them were analyzed later by Nicholls et al. (2010), , and Nie et al. (2017). Their studies showed that HBSs containing a red giant (RG) typically have longer periods than classical ellipsoidal variables. The vast majority of HBSs in our sample belong to this group; therefore, we could verify that HBSs are indeed an extension of the classical ellipsoidal systems to longer periods for a given brightness.
The majority of papers on HBSs published so far have focused on individual objects, e.g., KOI-54 (Welsh et al. 2011), KIC 3749404 (Hambleton et al. 2016), KIC 8164262 (Hambleton et al. 2018), and MACHO 80.7443.1717 (OGLE-LMC-HB-0254; Jayasinghe et al. 2019, 2021. Due to the limited number of known HBSs, a detailed quantitative analysis has been challenging to carry out. A sketchy analysis of HBSs containing RG stars was introduced in the work of Soszyński et al. (2004), where the authors explained the reason for the unusual shape of the light curves of the ellipsoidal variables. Later, a subset of those systems were analyzed by  and then by Nie & Wood (2014). In both works, the authors measured changes in radial velocities confirming the binary nature of those stars and showing that the strongest brightness changes occur near the periastron passage. Afterward, Nie et al. (2017) presented an analysis of 81 ellipsoidal RG stars, including 22 HBSs. They mainly focused on the evolutionary status of those stars and the primary/secondary mass distributions. In turn, Beck et al. (2014) studied the asteroseismic properties of 18 RG HBSs found in the data obtained by NASA's Kepler Space Telescope (Borucki et al. 2010).
Another set of HBSs, but consisting of stars located on or close to the MS, was found in the Kepler database (17 objects were the subject of an analysis by Thompson et al. 2012). The catalog, including HBSs, among other variable stars, was released by Kirk et al. (2016). Those HBSs are mainly low-and intermediate-mass A-F-type stars. They are characterized by a short orbital period (days to tens of days) and very small amplitudes of brightness variations (a few millimagnitudes) contrary, to the HBSs with an RG star, which have much longer periods (hundreds of days) and an order of magnitude higher amplitudes.
Recently, using the 9th Catalogue of Spectroscopic Binary Orbits (Pourbaix et al. 2004), Ko laczek-Szymański et al. (2021) selected HBS candidates and examined their light curves delivered by the Transiting Exoplanet Survey Satellite (TESS) mission from sectors 1 -16 (mostly the southern ecliptic hemisphere). The authors discovered 20 massive and intermediate-mass HBSs, seven of which exhibit several TEOs lying at low harmonics of the orbital frequency.
In this work, we conduct a general analysis of 991 HBSs cataloged in the OGLE Collection of Variable Stars (OCVS; Wrona et al. 2022, hereafter Paper I).
In parallel to the analysis of the heartbeat phenomenon itself, we search for TEOs in the presented collection of OGLE HBSs. The derived sample of binaries exhibiting TEOs may be a valuable test bed for future studies on the influence of dynamical tides on the orbital evolution in binary systems, including those with evolved companions.
The structure of this paper is as follows. In Section 2, we describe the origin of the photometric data of the HBSs and the way we prepare them for the examina-tion. In Section 3, we describe the modeling process of the light curves using the K95 model. In Section 4, we study the impact of irradiation and the reflection effect on the light curve of an HBS. In Section 5, we introduce the method used to search for TEOs and present the results of the search. The core of our work is presented in Section 6, where we discuss the results of the analysis. In Section 7, we summarize and conclude our work.

PHOTOMETRIC DATA
In this work, we analyze the sample of HBSs found in the OGLE project database. The detailed specifications of the data are presented in Section 2 of Paper I.
The time-series data used in the analysis were obtained using the 1.3 m Warsaw Telescope located at Las Campanas Observatory, Chile. The majority of the data come from the fourth phase of the OGLE project, in operation since 2010 (Udalski et al. 2015). The observations were conducted using Cousins I-band and Johnson V -band filters. The photometry was obtained using differential image analysis (DIA; e.g., Alard & Lupton 1998, Woźniak 2000. All of the data for each filter were calibrated separately to the standard photometric system using the scheme presented by Udalski et al. (2015). We also corrected the photometric uncertainties based on the work of Skowron et al. (2016).
In the analysis and modeling of the light curves, we used mainly the I-band data. The V -band photometry was used to determine the (V − I) color information for all HBSs from our sample. To prepare the data for modeling and analysis, we have taken the following steps.
First, we cleared the light curves from outliers by rejecting all of the data that lay outside the 3σ level from the average flux, where σ is the standard deviation of the flux in the entire range of time. This step was taken separately for all of the available data obtained by different OGLE phases (OGLE-II, OGLE-III, and OGLE-IV) because in some cases, there was a significant shift in the mean flux between those data. During this step, we also removed sets of data consisting of clearly improperly determined photometry, which could be caused, for instance, by bad weather conditions during the observing night or some failure in the DIA pipeline (this is a very common case for sources with high proper motions).
The second step was to remove trends in the data sets. This procedure was done separately for data obtained during OGLE-II, OGLE-III, and OGLE-IV. To each part of the light curve, we fitted cubic splines. Then, the obtained sets of splines were subtracted from the data.
In the third step, we shifted the detrended data in the fluxes from OGLE-II and OGLE-III to OGLE-IV (or   to OGLE-III if there were no data from OGLE-IV). For this purpose, we calculated the mean flux for each phase and shifted the data to the latest phase by the difference between these averages.
Finally, we used additional cleaning procedures. Using the combined data from all phases, we prepared light curves phase-folded with the orbital period and divided them into bins, the width of which was set on 0.1 of the orbital phase. Then we calculated a standard deviation in each bin and removed points lying more than A · σ from the mean magnitude. The A parameter was chosen individually for each star and mainly depended on the number and positions of outlying points. The A parameter was usually about 3.
To assess photometric temperatures for the HBSs containing a hot primary, we used U BV photometry obtained by Massey (2002). Data were taken with the Curtis Schmidt telescope at Cerro Tololo Inter-American Observatory (CTIO), Chile, using a Tektronix 2048 × 2048 CCD with a 2.32 /pixel scale. Observations were made at the beginning of 1999 and in 2001 March/April.
For some objects, the aforementioned U BV data were unavailable; therefore we decided to use U BV photometry obtained by Zaritsky et al. (2004). The authors used the 1 m Swope Telescope, located at Las Campanas Observatory, right next to the Warsaw Telescope. Observations were taken with the Great Circle Camera with a 2K CCD with a 0.7 pixel scale between 1995 October and 1999 December, with additional observations in 2001 December.
In this work, we also present, among others, periodluminosity (PL) diagrams using the W JK Wesenheit index. To calculate this quantity, we used JHK s photometry collected by Kato et al. (2007) using the InfraRed Survey Facility (IRSF) 1.4 m telescope at Sutherland, the South African Astronomical Observatory, with the SIRIUS camera, which is equipped with three 1024 × 1024 HgCdTe arrays. The pixel scale for this camera is 0.45 .
As complementary data, we utilized photometry collected during the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006Skrutskie et al. , 2019 In this work, we decided to use the analytic model of tidally induced stellar deformation presented in Kumar et al. (1995), where the normalized flux variations are described by their Equation (44). The K95 model was derived under certain assumptions and simplifications, such as spin-orbit alignment, inclusion of dominant modes only (spherical harmonics with l = 2, m = 0, ±2), and alignment of the tidal bulge with the line connecting the mass centers of stars and without the irradiation and Doppler beaming effects. Therefore, even a model perfectly fitted to the light curve only gives an assessment of the orbital parameters. The slightly modified version of the K95 formula was used by many authors (e.g., Thompson et al. 2012;Jayasinghe et al. 2019;Ko laczek-Szymański et al. 2021) to assess the basic orbital parameters of the system based on the light curve. According to Thompson et al. (2012), the relative change of the flux caused by tidal deformation as a function of time, t, can be expressed as where S is the scaling factor of the amplitude, C is the zero-point offset, i is the inclination of the orbit, ω is the argument of the periastron, ϕ(t) represents the true anomaly as a function of time, R(t) describes the distance between the components of the system as a function of time, and a is the semimajor axis. Both time-dependent variables R(t) and ϕ(t) can be rewritten as functions of the eccentricity, e, and eccentric anomaly, E: The quantities E and t are connected by Kepler's equation, where T 0 is the time of periastron passage, and P is the orbital period. However, we believe that there is a mistake in Equation (1). In the top panel of Figure 1, we present the light curve of the HBS OGLE-BLG-HB-0081. In the modeling process, we used Equation (1) and the Markov Chain Monte Carlo (MCMC) fitting procedure (described in more detail in Section 3.2). We obtained the following orbital parameters: e K = 0.803, i K = 64 • , and ω K = 17 • . The fitted model is represented by F F (t) = S 1 3sin 2 isin 2 ( (t) ) (R(t)/a) 3  (1), while in the bottom panel, we used the corrected version shown in Equation (5). We obtained the orbital parameters of the system using two independent methods: the K95 model (red line, index K) and the PHOEBE Legacy program, which is based on the WD code (blue line, index WD). the red line in Figure 1. We also performed modeling of that light curve using the PHysics Of Eclipsing BinariEs (PHOEBE) Legacy program 1 (Prša & Zwitter 2005), which is based on the Wilson-Devinney (WD) code (Wilson & Devinney 1971), and thus independent of the K95 model. The resulting values of e WD = 0.78 and i WD = 55 • are similar to the ones obtained with the K95 model (within 3σ), but the argument of periastron is far from that: ω WD = 152 • . Setting these values to Equation (1) and plotting the results in Figure 1 we got the blue line. It is clearly seen that this model is far from being proper, but one can notice that the blue line is a symmetric reflection of the red one about the vertical axis.
We also compared light curves generated on the basis of Equation (1) with synthetic light curves presented in Figure 5 of Thompson et al. (2012). The results turned out to be very similar to the ones described above. Our light curves were symmetric reflections of the ones presented by Thompson et al. (2012). We are certain that the fault is connected with the sign of the ω in Equation (1) because by inverting this sign, we get proper models. We were not able to track down the source of this discrepancy. Nevertheless, the version of Equation (1) with the plus sign before ω makes the model fit properly to the data and is in agreement with the synthetic light curves generated with alternative methods; thus, we decided to use the following equation instead of Equation (1): The default unit of a star's brightness in our catalog is magnitude, while during the modeling, we operated with a fractional change of the flux. We can express the change of the magnitude using Pogson's equation, where we used the median magnitude as a zero-point, m 0 . Thus, we can calculate the fractional change of the flux using following formula: σ f = 2.5(δF/F )σ m ln 10 , where σ f is an uncertainty of the flux change with a given uncertainty of the magnitude change, σ m .

Fitting the K95 Model to the Light Curve
In the fitting procedure, we decided to use the MCMC method, first to search for a proper model and investigate plausible degeneracies and then to estimate the uncertainties of the model's parameters. We used Python's emcee v3.0.2 package, described in detail by Foreman-Mackey et al. (2013).
The time span of the OGLE HBS light curves often exceeds a dozen years; therefore, in the case of a high apsidal motion rate, the shape of the light curve can change significantly (e.g., Hambleton et al. 2016). That may degrade the quality of the K95 model fit. To estimate the role of apsidal motion, we visually compared phasefolded light curves from the first and last three observational seasons. After the eye inspection, we detected clearly visible changes in the shape of the light curve in only six objects: OGLE-BLG-HB-0240, OGLE-BLG-HB-0317, OGLE-BLG-HB-0358, OGLE-LMC-HB-0072, OGLE-LMC-HB-0146, and OGLE-SMC-HB-0001. A light curve of the last object and a brief description of the possible causes of the light curve's shape changes are presented in Section 4.3.4 in Paper I. We did not detect any significant changes in the rest of the HBSs; thus, we decided to neglect the apsidal motion in the fitting procedure.
3.2.1. Uniqueness of the Solution Found Using the K95

Model
To assess the risk of degeneracy of the K95 model, we created 20,000 synthetic light curves using the K95 model with randomly selected parameters from the entire hyperspace. We also took into account the noise of the brightness and the uneven time sampling for the typical OGLE light curve. We assumed the normal distribution of the noise with the standard deviation depending on the mean magnitude, which was also randomly selected. The time sampling included the mid-season gaps and typical cadence of the observations for the central regions of the Magellanic Clouds (MCs) during the OGLE-IV project.
In the fitting process, we used the MCMC method. In the MCMC run, we used the K95 model, according to the formulae introduced in Section 3.1. We adopted the flat prior distribution for e, i, and ω. We limited their final values to the physically reasonable range but with a small margin for angular variables to avoid sharp cuts at the final distribution: 0 < e < 1, 0 • < i < 92 • , and −10 • < ω < 190 • . The range of the argument of periastron is only a half of the full angle because the stellar distortion due to tidal deformation is symmetric along the elongation axis. For the remaining parameters, we used the normal prior distribution. The range limits for P and T 0 were connected to their estimated initial values P e , T 0,e : 0.95P e < P < 1.05P e , T 0,e − 0.5P e < T 0 < T 0,e + 0.5P e . In the modeling process, we used 50 walkers and 20,000 steps for a single chain. As a result, we used parameters from the model with the maximum likelihood.
We found that systems with low inclination angles (i 20 • ) are difficult to model, and the obtained orbital parameters are usually unreliable. This behavior is caused by the sin 2 i term in the K95 model (Equation 5). In our sample, however, there is a very low number of HBSs with such a low inclination angle (see Figure 8); thus, to assess the level of the degeneracy of the K95 model, we will focus on the systems with i > 20 • .
For 90 % of simulated light curves with i > 20 • (≈ 16, 000 systems), the MCMC fitting procedure returned the correct values of the parameters within the 3σ region. For the remaining 10 % of the simulated light curves, we did not recover the initial parameters of the K95 model. In one of the most common cases (4 % of the sample), the program found two equivalent solutions for ω, one about 0 • and the second about 180 • . These two solutions are mathematically indistinguishable due to the 180 • ambiguity, which is the result of the periodicity of the sin 2 (ϕ(t) + ω) term of the K95 model. We also noticed an excess in the number of systems (about 4.5 % of the sample) for which the program found i ≈ 90 • , despite the true value of inclination often being far from 90 • . The rest of the parameters obtained for these systems (except for P and C) often differ significantly from the original values. We found that such behavior is exhibited by systems in which the mean scatter of the light curve is comparable with the amplitude of the heartbeat. The third most common case (1.5 % of the sample) was observed for systems with low inclination angles (i 30 • ) and the argument of the periastron about 0 • or 180 • . If we take a smaller i of about 10 • and a larger e of about 0.3 for e ≈ 0.2 or about 0.1 for e ≈ 0.7 (the value of the correction is inversely proportional to the eccentricity), and if we change ω to the right angle and lower amplitude scaling factor by a factor 3.5, the obtained model will be similar to the initial one.

Modeling of the OGLE HBS Light Curves
We conducted fitting with the MCMC method for all 991 HBSs from our sample. We used light curves cleaned from outliers, detrended, and shifted in average flux to the latest OGLE phase, as described in Section 2. Stars with additional brightness variations, such as eclipses or spots, were manually cleaned by removing the affected part of the light curve if it could be easily separated from the heartbeat modulation. Nevertheless, in some cases, this procedure could affect the heartbeat shape; therefore, the final parameters may not be reliable. Each OGLE HBS for which we did not find a proper K95 model or the fitted model is unreliable has been appropriately flagged (see the description of Table A2 in Paper I).
We ran the fitting process with similar prior distributions of parameters and their limits as described in Section 3.2.1. We calculated initial orbital periods, P e , using the FNPEAKS 2 program (created by Z. Ko laczkowski, W. Hebisch, and G. Kopacki), which is based on Fourier frequency spectra. The initial time of periastron passage, T 0,e , was estimated based on the local extrema in the light curves.
In the MCMC fitting process, for each star, we applied 50 walkers with 20,000 steps in a single chain. Only about 6% of the HBSs showed signs of degeneracies. In most cases, the problems were ω bouncing between 0 • 2 http://helas.astro.uni.wroc.pl/deliverables.php?active=fnpeaks and 180 • or the program finding a different solution for the period. For those problematic stars, we used another MCMC run with a larger number of walkers and steps and narrower prior distributions. This approach led to proper models for most of the troublesome light curves. We could not find the satisfactory model for only about 2% of our HBSs.
Finally, we once more used MCMC fitting to derive the uncertainties of the K95 model's parameters. The priors were randomly selected using a normal distribution centered on the value from the previous step, with the standard deviation being one-hundredth of this value (for T 0 , we used P/100, and for C, we used C/1000). We used 100 walkers and 25,000 steps in a single chain.
As a final value of the given parameter, we used the maximum value of the distribution, but we also considered using the median value. In most cases, these two approaches give similar values, but for the non-Gaussian distribution (e.g., when the i value is about 90 • ), we found that the first one results in better-fitted models. As uncertainties σ − and σ + , we assigned distances between the 50th and the 16th and 84th percentiles of the distribution, respectively. In Figure 2, we present a typical corner plot (top panel) and the corresponding phase-folded light curve including the model (bottom panel). Both corner plots and light curves with the fitted model for each OGLE HBS are available on the OGLE websites 3 . In most cases, there are no clearly visible correlations between parameters, except for the following pairs: i-S, i-C, ω-T 0 , and less often P -T 0 .
The final values of the K95 model parameters for each star with their uncertainties are presented in Paper I (their Table A3).

THE IMPACT OF THE IRRADIATION/REFLECTION EFFECT
In order to quantitatively estimate the impact of the irradiation/reflection effect on the shape of the OGLE I-band light curves of HBSs located in the Large Magellanic Cloud (LMC), we ran a set of dedicated simulations. We used PHOEBE 2 modeling software (version 2.3; Prša et al. 2016a;Horvat et al. 2018;Jones et al. 2020;Conroy et al. 2020) for this purpose and performed two groups of simulations. The first one refers to HBSs with the primary component being an RG. The majority of our HBSs fall into this group. In the second group, we considered HBSs with a massive primary component located on the MS or passing through the Hertzsprung gap. We will refer to these two groups of simulations as S-RG and S-MS, respectively. The S-RG simulations were done as follows. We considered three masses of RGs being a primary component, 1, 2, and 6 M with corresponding radii of 50, 75, and 100 R (see right panel of Figure 10). For each mass of the primary component, we generated 10,000 binaries with mass ratios q, e, i, and ω, and periastron distances relative to the sum of the radii of components r peri /(R 1 + R 2 ) drawn from uniform distributions U [α,β] on the interval [α, β]. We comment on them below.
• q ∼ U [0.2,0.9] -We did not consider q > 0.9 because we realistically assumed that the secondary companions are MS stars. Here q ≈ 1 would suggest that the secondary is also an RG, which in turn would result in very long orbital periods, unobserved by us. On the other hand, the overall strength of the tides is proportional to q; therefore, we omitted systems with q < 0.2. The typical peak-to-peak amplitude of the heartbeat in our sample of RGs is relatively high, ∼ 0.05 mag; hence, it rather excludes the possibility of low-q companions.
• e ∼ U [0.15,0.70] -We adopted a representative range of eccentricities observed in the analyzed collection of OGLE HBSs.
• r peri /(R 1 + R 2 ) ∼ U [1,4] -The strength of the tidal forces at periastron falls off rapidly as r −3 peri . Therefore, we set a maximum limit of r peri /(R 1 + R 2 ) to 4 in order to simulate only the orbits with a chance to reproduce detectable heartbeat signals.
The values of q, i, ω, and r peri /(R 1 + R 2 ) were generated independently from each other. The drawn values of r peri /(R 1 + R 2 ) were subsequently transformed to P , using Kepler's third law. During the simulations, only those systems were accepted that did not overflow at the periastron; i.e., the detached geometry was always preserved. The remaining physical parameters of the components, namely, their effective temperatures and the radii of the secondaries, were taken from the MIST isochrones of the age compatible with the parameters of the primary RG star. We assumed [Fe/H] = −0.7 for systems with 1 and 2 M primaries, while for 6 M we adopted a higher value of [Fe/H] = −0.4 in order to account for the metallicity gradient observed in the LMC. The bolometric albedos were set to 0.6 and 1.0 for components with convective and radiative envelopes, respectively. The surfaces of both components were simulated within PHOEBE 2 with 4000 -6000 triangular elements, depending on the drawn orbital parameters. The radiative properties of stars with effective temperatures above 4000 K were obtained from ATLAS9 model atmospheres (Castelli & Kurucz 2003). Otherwise, the PHOENIX models (Hauschildt et al. 1997, Husser et al. 2013 were used for cooler components. Both grids of model atmospheres are incorporated into PHOEBE 2 with the accompanying limb-darkening tables. Finally, for each system, two I C -band light curves were calculated in the full range of orbital phases, φ; one with the irradiation/reflection effect being switched off, F irr−off (φ), and the second one with the aforementioned effect treated in the formalism developed by Wilson (1990), F irr−on (φ). Since our calculations in PHOEBE 2 were performed in the absolute scale, the following condition was satisfied for every φ, F irr−on (φ) > F irr−off (φ). This is due to the fact that the irradiation/reflection effect can only add flux to the beam. Next, we were searching for the largest difference between these two synthetic light curves, which we denote as max(∆m irr ). We defined this quantity in the following way: The upper part of Figure 3 summarizes the effects of the S-RG simulations. The first thing that draws attention is the clear correlation between q, r peri , and max(∆m irr ) (second column in Figure 3). The larger the q and lower the r peri , the more pronounced the impact of the irradiation/reflection effect. Nevertheless, for systems with the primary's mass 2 M , the max(∆m irr ) 1 mmag. This fact allows us to state that the K95 model should return trustworthy orbital parameters for the majority of the OGLE HBSs containing an RG. This is because most of them are low-mass stars, as indicated by their positions on the Hertzsprung-Russell (H-R) diagram ( Figure 10). The situation is different for massive RGs, companions of which may have high effective temperatures and luminosities; hence, the irradiation/reflection effect is significant. In such a scenario, the 95th percentile of max(∆m irr ) is equal to circa 10 mmag, but systems with max(∆m irr ) up to 40 mmag can occur. Depending on the amplitude of the heartbeat observed in HBSs that contain the massive RG, one should be aware that the irradiation/reflection effect may become comparable to the ellipsoidal variability. Therefore, the orbital parameters obtained for such systems via fitting the K95 model to their light curves should be treated as estimates (see Section 4.1).
The S-MS simulations were performed analogously to the S-RG, but here we considered different masses of the primary components, 4, 10, and 25 M (see left panel of Figure 10). To realistically reflect the evolutionary phase of these components visible in Figure 10, we assumed the 4 M primary to be in the middle of its Hertzsprung gap, while the 10 and 25 M primaries were located halfway between the zero-age MS (ZAMS) and terminal-age MS (TAMS). All models were assumed to have [Fe/H] = −0.4. The rest of the parameters and methods necessary to run S-MS simulations were identical to the S-RG setup described above.
Similarly to the upper part of Figure 3, the lower part shows the results of the S-MS simulations. We found that for massive HBSs that are still on the MS or near it, the irradiation/reflection effect in the I C band is generally nonnegligible. The 95th percentile of max(∆m irr ) is equal to around 9 mmag for all variants in these simulations. However, still numerous systems can have an amplitude of the irradiation/reflection effect in the range of 10 -30 mmag. Considering these facts, we would like to emphasize that orbital parameters derived by us for massive MS HBSs may deviate from actual values and should be treated with caution (see Section 4.1).

K95 Model Fitted to the Synthetic Light Curves
In the previous subsection, we were interested in the overall amplitudes of the irradiation/reflection effect among the HBSs from our sample. However, how this phenomenon will affect the derived orbital parameters is another question. Recalling that the K95 model neglects the aforementioned effect, one can expect that some systematic biases in the derived parameters might occur. Moreover, in some circumstances, the morphology of the model may not be able to effectively reproduce the actual brightness changes. To explore both this issue and the properties of the model, we performed two types of tests, which were done on the synthetic irradiated I C -band light curves from the S-RG and S-MS simulations.
First, we examined if the K95 model can successfully fit the synthetic light curves provided that the orbital parameters are fixed on the values injected into PHOEBE 2. In other words, the only free parameters during the fitting procedure were S and C (see Equation 5). In order to quantitatively describe the global effectiveness of the best fit in the entire range of orbital phases, we constructed the following metric: where K best−fit (φ) denotes the K95 model that best fits the F irr−on (φ), which stands for (F irr−on (φ)/ F irr−on − 1). The quantity η is a ratio between the integral of the residuals from the fit and the integral of the normalized synthetic light curve; therefore, it is independent of the amplitude of the heartbeat. The leftmost column in Figure 4 (titled "K95 fixed") presents the distribution of the mean value of η, η in the hex-binned i -e plane, obtained for the S-RG and S-MS sets of simulations. For both types of simulations, regardless of the primary's mass, it can be seen that the K95 model with fixed orbital parameters is able to successfully fit the simulated variability when i 15 • . For greater values of i, the results are always worse because the contribution of irradiation/reflection to the total light curve can no longer mimic the ellipsoidal variability. This is especially pronounced in the S-MS simulations (lower part of Figure 4) for high-e binaries and i 40 • . Another important feature is the significant difference in η between the 4 and 10/25 M S-MS binaries. The former are generally characterized by a much better quality of the fit. The situation is notably different for the S-RG simulations (upper part of Figure 4). The worst matches lie in a vertical strip with inclinations in the ∼20 • -50 • range, almost independent of the eccentricity of the system. This is in contrast to the effects of the S-MS simulations discussed above. In an obvious way, the models with the primary's mass of 6 M exhibit far larger departures from the K95 model when compared to the less massive S-RG systems. This is due to much greater amplitudes of the irradiation/reflection effect for this kind of HBS (see Figure 3, upper part). Nevertheless, the values of η are statistically significantly smaller for binaries containing an RG as a primary component than those that harbor two MS stars.
Our next test was analogous to the previous one, except for one major difference. This time, the orbital parameters were treated as free parameters in a leastsquares fit together with S and C. Thanks to this approach, we were able to track down any systematic biases in the determination of the orbital parameters. The results of our second test are summarized in Figure 4. Similarly to the leftmost column in Figure 4, the column denoted as "K95 free" shows the distributions of η across the i -e plane. It can be easily seen that the shape of these distributions remained nearly the same as in the previous experiment, but now the values of η are significantly smaller. This is because the K95 model is trying to fit both the ellipsoidal and irradiation/reflection variability. The price to pay for "forcing" this fit is some systematic errors in the optimized orbital parameters. Let us denote this error for a given orbital parameter, The series of hex-binned panels on the left shows the distributions of the average value of η, η (color-coded), for simulated orbits with different eccentricities and inclinations. The column of panels denoted as "K95 fixed" presents the results for the K95 model with fixed orbital parameters during the fitting procedure. In turn, the column denoted as "K95 free" corresponds to the fits with orbital parameters being free during the optimization. On the right, we present the absolute errors of the orbital parameters, ∆X , obtained from the "K95 free" simulations. The value presented in the corner of each plot is equal to the median of the |∆X |. See the discussion in Section 4.1 for more details. Note the different range of η for S-RG and S-MS. X , as ∆X = X PHOEBE − X K95 , where X PHOEBE is an orbital parameter used to generate a PHOEBE 2 synthetic light curve, and X K95 corresponds to the parameter obtained from fitting the K95 model to the synthetic data. We note that ∆X > 0 means that X is being underestimated by the K95 model, while ∆X < 0 suggests the opposite. The right-hand part of Figure 4 presents quadriads of panels showing ∆e, ∆i, ∆ω, and ∆T 0 /P for each primary's mass in the S-RG and S-MS simulations. They are plotted against i, which is the parameter most strongly correlated with these ∆X . The values presented in the corners correspond to the medians of |∆X |. There are several conclusions about the behavior of the K95 model that can be drawn from distributions of ∆X . We discuss them briefly below.
• The systematic errors in determining all four orbital parameters for the S-RG simulations are much smaller than those for the S-MS simulations.
Since the majority of our HBSs are systems with the primary being an RG, we infer that the K95 model returns reliable parameters for them.
• The most accurately estimated parameters, regardless of the strength of the irradiation/reflection effect, are e and T 0 . In the set of S-RG simulations, |∆e| < 0.03 and |∆T 0 /P | < 0.03 in almost all cases. For the S-MS simulations, the corresponding ranges are −0.04 < ∆e < 0.1 and |∆T 0 /P | < 0.03.
• In general, i is estimated with the lowest accuracy among all four orbital parameters, and its systematic error, ∆i, depends in a nontrivial way on the actual i of the system. For i 15 • , the K95 model has a tendency to overestimate i. In turn, ∆i for the orbits with relatively high values of i can behave in two ways. If the orbit is only slightly eccentric (e 0.2), it is very likely that the K95 model will return i ≈ 90 • . This can be seen from the diagonal line of points at the bottom of each ∆i distribution. On the contrary, if the orbit is more eccentric, the K95 model will certainly underestimate the actual value of i. This is because the irradiation/reflection effect fills the "dips" in the heartbeat signal, so the light curve seems to originate from a binary with a fictitious smaller i.
• As expected, ∆ω shows a typical dependence on i.
The argument of periastron is well constrained for orbits with relatively high i's and e's but becomes poorly constrained for small values of i. Eventually, it is undefined for i = 0 • . Let us also empha-size at this point that the K95 model allows for determining ω with a 180 • ambiguity.
Although the above analysis was performed for HBSs located in the LMC, its results should also remain valid for systems from the Small Magellanic Cloud (SMC) and Galactic bulge (GB). Thanks to the S-MS and S-RG simulations that we performed, we are able to estimate the average accuracy of the orbital parameters obtained by us via modeling the I C -band light curves of OGLE HBSs. Combining the results of all simulations, the median values of |∆e|, |∆i|, |∆ω|, and |∆T 0 /P | are equal to about 0.01, 3.0, 2.5, and 0.002, respectively. However, one should be aware that some individual situations can still be characterized by relatively large ∆X , especially when it comes to the determined i. In particular, the massive MS HBSs are vulnerable to such effects in our modeling.

General Properties of TEOs
The majority of TEOs known so far are forced damped gravity (g) modes that may dissipate the total orbital energy and make the system tighter with time. Therefore, TEOs may play a significant role in the dynamical evolution of the binary system. In general, TEOs come in two "varieties." The first one occurs when the frequency of the eigenmode temporarily coincides with the harmonic, n, of the orbital frequency, f orb . In such circumstances, which are called a "chance resonance," the amplitude of a TEO does not exceed the parts per thousand level and is most often significantly smaller. The other variety is associated with the so-called "resonantly locked modes" (see Fuller 2017;Hambleton et al. 2018). These are forced normal modes with frequencies evolving due to the stellar evolution at the same rate and direction as the nearest harmonic of the orbital frequency. Therefore, resonantly locked TEOs have enough time to gain relatively high amplitudes and hence effectively dissipate the total orbital energy. High-amplitude TEOs are expected to be quadrupole (l = 2) modes with a mainly azimuthal order m = 0 or 2. Although the history of theoretical studies of the dynamical tides dates back to the 1970s (Zahn 1970), the actual effectiveness of energy dissipation due to TEOs and their impact on the binary's evolution are still a matter of large uncertainties. Studies of large-amplitude TEOs offer a unique opportunity to make progress in this subject. The main observational difference between self-excited pulsations and TEOs is that the latter have frequencies exactly equal to the integer multiples 4 of f orb ; therefore, they phase well with the orbital period, P (see the top left panel in Figure 5). Nevertheless, because of the resonant nonlinear mode coupling (NLMC, described extensively by, for instance, Dziembowski 1982;Dziembowski & Królikowska 1985;Dziembowski et al. 1988), it is also possible to observe nonharmonic TEOs that are "daughter" modes of the resonant (harmonic) TEO, which is the "mother" mode (e.g., Guo 2020). The simplest manifestation of this mechanism is the decay of the resonant TEO into two modes, the sum of the frequencies of which equal n·f orb . However, one can also expect the quintuplets, septuplets, etc. of frequencies formed via higher-order NLMC and/or non-harmonic TEOs being e.g. a "granddaughter" modes (i.e. the "daughter" modes of some prior "daughter" mode).
There are several other features that make TEOs unique when compared to the self-excited oscillations. Therefore, the whole branch of research on this subject is called "tidal asteroseismology." One of the most extensively studied HBSs in terms of tidal asteroseismology is the aforementioned KOI-54, for which this type of analysis was performed by Burkart et al. (2012). There are also other papers dedicated solely to this topic, e.g. Guo et al. (2020).
While the TEOs observed in the MS stars are mostly g modes, the situation is diametrically different for the RGs, which make up the vast majority of our sample of HBSs. Both the Brunt-Väisälä and Lamb frequencies are very large inside the dense core of an RG; therefore, a nonradial oscillation mode has a dual nature in its interior. It propagates as an acoustic mode in the convective envelope of an RG but as a g mode in its radiative core. This is the reason why these kinds of modes are called mixed modes (Aerts et al. 2010). They were detected in many RGs, mainly thanks to the ultraprecise photometry delivered by the Kepler mission (see, e.g., Beck et al. 2011;Hekker & Christensen-Dalsgaard 2017, and references therein). The mixed modes are generally problematic when it comes to their analysis and modeling, for several reasons. Their wavelengths are extremely short in the core of an RG; hence, the corresponding eigenfunctions are characterized by a large number of nodes in the radial direction (typically of the order of 10 3 ). Next, they are expected to undergo severe radiative damping when traveling through the dense core (Dziembowski 1971) and be subject to significant nonlinear effects in the upper part of the red giant branch (RGB; Wein-4 Assuming the linear theory of tides. In turn, nonlinear tidal effects can lead to a minor offset between the observed frequency of a TEO and its corresponding harmonic of the orbital frequency. berg & Arras 2019). Dupret et al. (2009) performed the theoretical investigation of amplitudes and lifetimes of the self-excited radial and nonradial (l = 1, 2) oscillations in the RGs (their "model C" is the closest to the RG HBSs analyzed in our work). The authors found that the aforementioned properties of mixed modes are a sensitive function of the density contrast between the core and the envelope (i.e., the position in the RGB). The amplitudes and lifetimes of nonradial modes also depend on which part of an RG the mode is trapped in. As expected, the modes trapped in the core should have small amplitudes on the surface, in contrast to the modes trapped in the extended envelope. Considering these facts, one does not expect to observe many nonradial modes in the RG stars, especially located high in the RGB, provided that they are stochastically driven by the turbulent convection. An interesting question arises, however: how will this picture change for the mixed modes that are excited by periodic and coherent variation in tidal forces (see, e.g., Fuller et al. 2013)?
More importantly, what is the impact of tidally excited mixed modes on the dynamical evolution of the HBSs with an evolved component(s)? In order to help answer the questions raised above, we provide the community with the compilation of high-amplitude TEOs detected in the OGLE HBSs (see Appendix A, Tables A1 and A2).

Methodology of Search
We search for TEOs in our sample of HBSs by means of the Fourier analysis. We performed a standard iterative prewhitening procedure on the residual light curves obtained after the subtraction of the best-fitting K95 model. The vast majority of the residual light curves reveal the presence of long-term variability that, regardless of whether it is physical or not, significantly enhances the signal in the frequency spectra at low frequencies. In order to get rid of these long-term brightness changes, prior to the prewhitening, it was modeled with Akima cubic splines (Akima 1970) and subtracted from each residual light curve. After calculating the error-weighted Fourier frequency spectra, only peaks with a signal-tonoise ratio (S/N)≥ 4 were considered statistically significant. The mean noise level, N , was derived as the mean signal in the frequency spectrum in the frequency range 0 -10 day −1 . The final parameters characterizing the coherent variability were obtained by means of the error-weighted nonlinear least-squares fitting of the truncated Fourier series to the corrected residual light curves. The formal errors of the extracted frequencies and their amplitudes were estimated using the covariance matrix. The frequency f was considered a TEO if it was sufficiently close to the nearest harmonic of the orbital frequency, i.e., if the following condition was satisfied: |f /f orb − n| ≤ 3σ f /f orb . We estimated the error of f /f orb , σ f /f orb using a standard error-propagation formula, σ f /f orb = (P 2 σ 2 f + f 2 σ 2 P ) 1/2 , where σ f and σ P stand for the error of f and P , respectively. Recalling that TEOs can also have a nonharmonic nature (provided that they formed via NLMC), we examined if any sum of two nonresonant frequencies, f 1 and f 2 , fulfills the inequality |(f 1 +f 2 )/f orb −n| ≤ 3σ (f1+f2)/f orb , where σ (f1+f2)/f orb was calculated analogously to σ f /f orb . We did not look for the presence of nonharmonic TEOs resulting from higher-order NLMC due to the insufficient precision of the fitted frequencies.

Results of Search
The results of our search are presented in Tables A1 and A2. In total, we were able to find 52 systems (five on the MS and 47 containing a post-MS star) out of 991 (∼ 5 %) that exhibit at least a single TEO, while the total number of detected TEOs amounts to 78. Figure 5 shows the compilation of frequency spectra of four sample OGLE HBSs with detected TEOs. The lefthand side of Figure 5 refers to one of the most prominent TEOs detected by us in OGLE-BLG-HB-0451. The outof-periastron variability of this object reveals the clear beating pattern between the dominant n = 26 and 30 TEOs. We also found evidence that 12 nonresonant frequencies present in four systems, namely, OGLE-BLG-HB-0066, OGLE-BLG-HB-0157, OGLE-BLG-HB-0208, and OGLE-BLG-HB-0362, are possible nonharmonic TEOs formed via the NLMC. We did not detect their corresponding "mother" modes; however, this can be explained with their amplitudes below the detection limit (typically between 0.2 and 3.0 mmag in the analyzed OGLE light curves).
Additionally, we have marked the systems in Tables A1 and A2 that, in parallel to the TEOs, also show a pronounced intrinsic periodic variability 5 . This may help in future research about the interaction between tides and intrinsic pulsations. Do tidal interactions suppress self-excited pulsations, or do they not have a major impact on them (e.g., Springer & Shaviv 2013, Fuller et al. 2020? The cases in which we did not detect any significant periodic variability are also interesting because of the question raised above. The presented sample of TEOs is the largest homogeneous sample of this kind known so far, which allows us to statistically investigate the dependence between parameters like n, e, and the amplitude of a TEO. From a theoretical point of view, the amplitude of a TEO depends on several factors that were described in detail by Fuller (2017) (his Equation (2) and related equations). Nevertheless, one can estimate the range of n, which favors excitation of a TEO, considering only the product of the so-called overlap integral, Q kl , and the Hansen coefficient, X nm . The former describes the spatial coupling between the tidal potential and a given oscillation mode of radial order k and degree l. Thus, it generally has the greatest values for small |k| and l. Under the assumption of spin-orbit alignment in the system, X nm = W lm F nm . Supplementary to Q kl , F nm describes the temporal coupling between the mode and characteristic time of periastron passage. It is given by the following expression (Fuller 2017, his Equation (5)): (1 − e cos E) l dE. (11) Here W lm is a constant depending only on the geometry of a given mode. Keeping the assumption about spin-orbit alignment, W 20 = −(π/5) 1/2 , while W 22 = (3π/10) 1/2 . Figure 6 (middle and bottom panels) shows the evolution of log(|X nm |) with increasing eccentricity for quadrupole modes with m = 0 and 2. The maximum of |X n0 | is always located at n = 1, regardless of the eccentricity, but the values of |X n0 | become greater in the entire range of n for greater eccentricities. The m = 2 modes are characterized by |X n2 |, which peaks at n > 1. Moreover, the position of the maximum moves to higher n with the increasing eccentricity. In principle, |Q kl | and |X nm | peak at different values of n; therefore, only their product and its maximum inform about the most favorable n for the occurrence of TEOs. Since Q kl does not depend on e, while X nm does, the maximum of |Q kl X nm | shifts toward greater n for increasing eccentricity (see Burkart et al. 2012, their Figures 2 and 3). Hence, we can expect that in binaries with greater eccentricity, we will observe, on average, TEOs with higher values of n. The empirical relation between the n of a TEO and the eccentricity can be seen in the top panel of Figure 6.  (Table A1  and Table A2) and their estimated eccentricities. Directly detected TEOs are denoted with filled circles. The higher the amplitude of a TEO, the bigger the circle. Open circles correspond to the TEOs deduced from the presence of 'daughter' nonharmonic TEOs. For more details, see Section 5.3. Middle: contour plot of log(|Xnm|) versus n and e for l = 2, m = 2 modes. Bottom: same as middle panel but for m = 0 modes.
First of all, theoretical predictions about the positive correlation between n and e, which we described above, are reflected in the obtained empirical distribution of TEOs. It seems plausible to claim that the observed distribution of n and e resembles the shape of |X nm | dis-tributions presented in the middle and bottom panels of Figure 6. The direct comparison between the theoretical distributions of |Q kl X nm | and the observed properties of TEOs would require the calculation of the overlap integrals, which is beyond the scope of this paper. Next, we did not find any clear relation between the amplitude of a TEO, e, and n, but yet we note that TEOs are not restricted to strongly eccentric binaries. In our sample, we can still observe high-amplitude TEOs even in the systems with e ≈ 0.1. Some of these TEOs are related to the surprisingly large harmonics of the orbital frequency, despite the low value of e. Finally, we did not find any obvious dependence between the location of the HBSs on the H-R diagram and the presence of TEOs (see Figures 9 and 10). The presented collection of evolved OGLE HBSs that exhibit TEOs is a step forward in the application of tidal asteroseismology to the RG stars and studies of the mixed-mode TEOs.

Comparison of Obtained e, i, and ω Parameters with Previous Works
To verify if the resulting parameters are consistent with true physical values, we have compared the solutions for the sample of our HBSs to the ones obtained using different methods, unrelated to the K95 model. For the comparison, we used the results of the work by Nie et al. (2017). The authors studied 81 ellipsoidal RG binaries cataloged in the OCVS. In the modeling process, based on the 2010 version of the WD code, they used both light (mainly from OGLE-II and OGLE-III) and radial velocity curves using data collected during their previous works (Nicholls et al. 2010;Nie & Wood 2014). For the details of the modeling process using the WD code, we refer the reader to Section 3 of Nie et al. (2017). The total sample of 81 ellipsoidal variables consists of 59 systems with circular orbits and 22 with eccentric ones. Among the latter group, 19 stars were present in our catalog. The three remaining systems have a low eccentricity (0.069, 0.061, and 0.054), and we have considered them as classical ellipsoidal variables.
The results of the comparison of the orbital parameters P , e, i, and ω are presented in Figure 7. The horizontal and vertical axes represent values obtained during our analysis and the one presented by Nie et al. (2017), respectively. In the diagram containing inclination, we did not include stars with a 90 • flag from Nie et al. (2017)  The results shown in Figure 7 confirmed that the K95 model is capable of finding reliable orbital parameters for HBSs containing an RG. In Section 4, we also discussed the reliability of the K95 model in the dependence on the irradiation/reflection effect.

Distribution of the K95 Model Parameters
In Figure 8, we present histograms of e, i, and ω. Each chart shows a distribution of a given orbital parameter, separately for stars located toward the GB and MCs (colored lines) and combined (black lines). An HBS must have an eccentric orbit by definition, and the higher the value of e, the easier it is to distinguish an HBS from a classical ellipsoidal variable; thus, the low number of HBSs with e 0.1 is well understood. The lack of stars with high eccentricity is the result of a natural tendency of binaries to circularize their orbit (e.g., Zahn 1975;Hut 1980). In the histogram of ω, we have not noticed any specific features. The distribution of this parameter is nearly uniform. Distribution of the observed inclination angles of the HBSs may be influenced by many factors. We discuss some of them below.
• The orbit of a binary system can be oriented in any direction; however, the distribution of the observed inclination angles would not be uniform and would favor higher values of inclination. This is because the probability of observing a system with an inclination angle between i and i+di is proportional to the area of the spherical sector between these angles. The area of the surface element on the unit sphere is proportional to the sine of the inclination angle; therefore, the area of the annular spherical sector bounded by angles i and i + di will be greater near the orbital plane than near the pole. The resulting inclination angle distribution turns out to be uniform in cos i.
• Kumar et al. (1995) mentioned that the amplitude of the light curve depends mainly on the mass and structure of the star, although if we view a system from different angles, the contribution to the observed light curve from various modes changes. For small i, the main contribution to the light curve comes from the m = 0 mode, while for higher i modes, m = ±2 starts to dominate. The change of the contribution between modes mainly impacts the shape of the light curve, but it also affects the brightness. Thus, if the observer sees only lowamplitude m = 0 modes, the brightness changes could be too small for proper classification of the star.
• Another factor that may skew the i and ω distribution is the misclassification of stars based on the shape of the light curve. The HBSs, especially with low inclinations, may mimic other types of variability, such as spotted stars (mainly Ap-type variables; e.g., Iwanek et al. 2019) or some kind of Be stars. On the other hand, for high i and for ω near 90 • , the light curve has a striking resemblance to the eclipsing binary (in the case of high eccentricity) or the classical ellipsoidal variable (for lower e).
• In the histogram of i, most prominent is the peak for i ≈ 90 • . First of all, it is observational bias, caused by the fact that the majority of OGLE HBSs were found in the catalogs of eclipsing and ellipsoidal variables, which exhibit clear minima.
Secondly, as we discussed in Section 4.1, for low eccentric orbits, due to the irradiation/reflection effect, the K95 model tends to overestimate the value of i and usually returns i ≈ 90 • , while for more eccentric orbits values of i are underestimated. This inaccuracy of the K95 model also explains a dip in the i distribution between 70 • and 90 • .

Color-Magnitude Diagrams
The majority of the previously known HBSs have been found in the Kepler mission database; therefore, those systems are mainly F-, G-, and K-class objects located on the MS. There is also a group of systems containing hotter and usually more evolved stars of classes O, B, and A (e.g., Ko laczek-Szymański et al. 2021, and references therein). About 10% of the OGLE HBS sample represents that group of stars. On the other hand, about 90% of our HBSs belong to the RGB or asymptotic giant branch (AGB) and, less frequently, to the horizontal branch. A large group of RG ellipsoidal variables was described by Soszyński et al. (2004). The majority of those systems have a circular orbit, but the authors also accentuated a group of systems with high eccentric orbits. Until our work, presented in Paper I, they were the largest collection of RG HBSs (officially cataloged as ellipsoidal variables).

Color-Magnitude Diagram for the GB
The sample of our HBSs located toward the GB is dominated by stars lying on the RGB and AGB, which are noticeable in the color-magnitude diagram (CMD) in the left panel of Figure 9. The background for the CMD has been created based on the calibrated photometric maps for several OGLE-IV fields located around the Galactic center. The I -band extinction values A I and color excess of stars E(V − I) were determined using extinction maps presented by Nataf et al. (2013). Extinction toward the GB is very heterogeneous what causes a widening of the color and magnitude distributions in the CMD, which is clearly visible on the red clump (RC; a wide group of points around [(V − I) 0 , I 0 ] = [1.05, 14.3]). If there was no extinction, most of the HBSs from our collection would be saturated because the brightness limit for the Warsaw telescope is about 13 mag in the I band. There are only two systems on the MS. Such a small number of HBSs in this region was very unexpected because most of the HBSs known before were in this place in the CMD. The absence of such HBSs in our catalog is probably the result of the data selection.

Color-Magnitude Diagram for the MCs
In the middle and right panels of Figure 9, we present the positions of the HBS samples in the CMDs from the LMC and SMC, respectively. Both CMDs have been created similarly to the CMD for the GB. We took calibrated photometric maps for several OGLE-IV fields and combined them, creating the background for the HBS sample. The E(V − I) color excess was calculated using the reddening map of the MCs (Skowron et al. 2021).
Contrary to the distribution of HBSs in the CMD for the GB, here we can easily distinguish at least two large groups of HBSs. The first one consists of hot MS and Hertzsprung-gap stars of spectral type from late O to F, which is consistent with HBSs from the Kepler sample. The second noticeable group (especially for the LMC) is located on the RGB and AGB, which is in agreement with the results for the GB. There are also less numerous groups, e.g., one HBS near RC for each location and a few stars on the horizontal branch.

Evolutionary Status of the OGLE HBSs
For a better picture of the evolutionary status of the OGLE HBSs, we describe them in the following two subsections. The first one refers to the stars located on or In the left panel, we show evolution tracks only to the He core-burning phase, while in the right panel, for the masses below or equal to 1 M , we show an evolution track only to the RGB tip (solid lines). For higher masses, we also plot the evolution after this phase (dashed lines). In both panels, with dashed black lines, we show lines of constant radius. See more details in the text.
near the MS, and in the second one, we discuss systems containing an RG star. Here we analyze stars located in the LMC only. The SMC sample is too small for such an analysis. In the GB, extinction is highly nonuniformly distributed, which hampers a precise determination of the color and, in turn, results in very high uncertainties in the effective temperature of stars. Moreover, for the LMC, we can assume that all stars are located at the same distance (the distances between stars are insignificant compared to the distance to the LMC), which allows us to determine an absolute magnitude for each star. In the GB, we do not know the distance for the majority of stars; thus, we are not able to calculate precise absolute magnitudes.

HBSs with an MS or Hertzsprung-gap Primary
The HBSs containing a hot star (T eff 8000 K) are located mainly on or near the MS (perhaps the Hertzsprung gap). The HBSs with cooler MS stars usu-ally show very small brightness changes (less than a few thousandths of magnitude; e.g., Kirk et al. 2016), which makes them very difficult to identify in the OGLE data. In panel (a) of Figure 10, we present the H-R diagram for the bluest part of our HBSs located in the LMC. To calculate the effective temperatures of stars, we used UBV photometry obtained by Massey (2002) and Zaritsky et al. (2004). We followed the prescription described in Massey et al. (1989), For two stars, the U BV photometry was unavailable. In order to estimate photometric temperatures for these stars, we used the T eff -(V − I) 0 grid collected in the work of Bessell et al. (1998) for early-type stars (their Table 1) and the surface gravity log g = 4.5.
To calculate luminosity, we used a simple Pogson's equation, where M bol = M I + BC I and M bol, = 4.74 mag (Prša et al. 2016b) are the bolometric absolute magnitudes of a star and the Sun, respectively. Assuming the distance to the LMC as d LMC = 49.59 kpc (Pietrzyński et al. 2019) and calculating the bolometric correction, BC I , for the I filter, we find M bol as Here I 0 is the mean I-band magnitude with subtracted extinction. We computed BC I values based on bolometric correction tables available on the MIST project website 6 . We assumed a metallicity [Fe/H] = −0.5 and log g = 4.5.
In Figure 10, with colored solid lines, we mark the MIST v1.2 evolutionary tracks by Choi et al. (2016) and Dotter (2016), computed with the Modules for Experiments in Stellar Astrophysics (MESA) code (e.g., Paxton et al. 2011Paxton et al. , 2013. For all tracks, we set an initial rotation of v/v crit = 0.4 and metallicity [Fe/H] = −0.4 (e.g., Harris & Zaritsky 2009). By comparing the position of HBSs in the H-R diagram with the evolutionary tracks of stars, we can see that HBSs have a wide range of possible ZAMS masses, from about 3M to as high as 40M .
There is also a noticeable trend: the cooler and less luminous the star, the farther it lies from the MS. The MS is shown as the shaded gray area in panel (a) of Figure  10. The lower and upper edges of this region represent the lines of the ZAMS and TAMS, respectively. The second notable observation is that the HBSs form a group along the line that is parallel to the line of constant radius. Moreover, we can see that the majority of HBSs lie in the channel between 5 and 20 R .
There are at least two possible explanations for the observed departure of less luminous OGLE HBSs from 6 http://waps.cfa.harvard.edu/MIST/index.html the MS. The first one is an observational bias caused by the dependence between the stellar radius and the amplitude of the heartbeat. In theory, the larger the stellar radius, the larger the tidal deformation of the stellar surface, and thus the larger the amplitude of brightness changes. Stars with M ZAMS < 10M begin their evolution on the MS with a radius not exceeding 4-5 R ; thus, the heartbeat amplitudes may be too low for detection in ground-based surveys like OGLE.
The second explanation is the increasing light contamination from the companion for less luminous primaries, which may cause a shift in the location on the H-R diagram. If we assume a system where the primary and companion stars have similar effective temperatures, the location of this system on the H-R diagram will be shifted upward relative to the position occupied by the primary itself. The dotted red line in panel (a) of Figure 10 represents the position of systems consisting of two TAMS stars with equal effective temperatures. On the other hand, if the companion has a lower effective temperature, the position of the entire system on the H-R diagram will be shifted to the right relative to the position of the primary. In Figure 10, the markers show the location of the primaries with the assumption that the light contamination from the companion is negligible, which may not be correct. This is probably the case for about half of the hot HBSs because they are also eclipsing binaries with clearly visible primary and secondary eclipses. To determine the amount of contamination from the companion, further analysis is necessary, especially with the use of spectroscopic observations, which would allow for a more accurate estimation of the effective temperature of the main component. However, such an analysis is beyond the scope of this work.

HBSs with an RG Primary
The most numerous group of OGLE HBSs is located in the RG part of the H-R diagram. Those systems generally consist of a low-mass (M ZAMS 2M ) or intermediate-mass (2M M ZAMS 8M ) primary, which is evolving through the RG region, and a less massive companion, which most likely belongs to the MS.
Unlike the H-R diagram for the hot stars described in the previous subsection, here we decided to obtain photometric temperatures based on the (V − I) 0 color derived directly from the OGLE data. We used an empirical relation from Houdashelt et al. (2000) for giant stars, which is reliable for a range 0.70 < (V − I) 0 < 1.68. Luminosities were calculated using Equations (13) and (14). We utilized the bolometric correction from the YBC database 7 , described by Chen et al. (2019).
In panel (b) of Figure 10, we present the H-R diagram for the RG part of the HBSs located in the LMC. Similarly to the blue part of the HBSs, here we also generated MIST evolutionary tracks. Solid lines represent stellar evolution until the RGB tip phase, and dashed lines show the subsequent stages (shown only for masses M ZAMS > 1M ) until the AGB tip. We used a metallicity gradient depending on the initial mass (e.g., Harris & Zaritsky 2009): (16) For each track, we set the initial rotation v/v crit = 0.0.
The obtained H-R diagram is fully consistent with the one presented in Figure 4 of Nie et al. (2017), even though we used a slightly different approach to construct it. Nie et al. (2017) used an evolutionary track from Bertelli et al. (2008), the log L/L and T eff parameters were obtained based on I-and K-band photometry using the prescription from Nie et al. (2012) and, last but not least, their H-R diagram refers mainly to the classical ellipsoidal variables (only 22 systems have an eccentric orbit), while our sample includes only HBSs. In both works, the largest number of stars occupy the region for initial masses of less than 1.85-2 M , which indicates that these systems contain the RG star that is evolving through the RGB (with a degenerate He core) or, in the case of stars more massive than the Sun, on the AGB (with a degenerate C/O core).
There is also a second group of HBSs, which is represented by stars with a larger initial mass. These stars can be in any medium stage of their life, from evolving on the RGB with a nondegenerate He core, through the He core-burning phase, to the evolution on the AGB. For stars with T eff 4250 K (log T eff 3.63), more favored is an option with AGB evolution because the tip of the RGB and the He core-burning phase for stars with M ZAMS 6M do not exceed such temperatures.

Period-Amplitude Diagrams
The orbital period and I-band amplitude distribution for the OGLE HBSs are presented in Figure 11. Our sample, contrary to HBSs from the Kepler data, mainly 7 http://stev.oapd.inaf.it/YBC/ consists of systems with a long orbital period (mostly a few hundred days) and high-amplitude brightness variations as for the HBS (a few hundredths of a magnitude and larger). The color scale reflects the effective temperatures of the stars based on the dereddened (V − I) 0 color index. Black ((V − I) 0 > 3.0) is an indication of stars located in the regions not included in the extinction maps.
The period and amplitude ranges for the cooler group of HBSs (yellow and red, types from G to M) are similar for all locations, which indicates the common nature of those systems. Hot stars (blue and cyan, types from late O to F), which correspond to HBSs located on the MS or in its vicinity, have a similar mean value of the amplitude but slightly bigger scatter than the cooler group of HBSs. However, the orbital period distribution for hotter stars is extremely different. Those systems have much shorter periods, from a few days to several tens of days. It is the result of the different sizes of the primaries. In Figure 10, we can see that the HBSs located on the blue part of the H-R diagram reach sizes from 5 to 25 R , while the HBSs from the red part of the H-R diagram have radii from 25 to even 200 R . Since stars in the system cannot come too close to each other, the smaller the stars, the smaller the minimum distance needed, and thus the shorter the orbital periods. In Figure 11, we can also see a trend in the period-color relations: the higher the value of the color, the longer the orbital period in the system.

PL Relations for RG HBSs
The RG variable stars, such as Miras, OGLE smallamplitude red giants (OSARGs), semiregular variables (SRVs), long secondary periods (LSPs), or long-period eclipsing and ellipsoidal stars exhibit classical PL relations. Within each of those classes, in the PL diagram, the stars assemble into linear sequences, marked with letters from A to E (e.g., Wood et al. 1999, Soszyński et al. 2004). The sequences differ in slope and spread, and they are also shifted relative to each other. These parameters also depend on the set of filters that were used (e.g., scatter is much smaller for near-IR filters than for optical ones). In the case of HBSs, the most crucial is sequence E, which is formed by eclipsing and ellipsoidal variables.
In Figure 12, we present PL diagrams for RG variables from the OCVS (references are listed in Table 1), including HBSs from our collection. We plotted PL diagrams separately for the LMC and GB using two types of the extinction-free Wesenheit indices W I and W JK , defined as where R I,V = A I /E(V − I) and R Ks,J = A Ks /E(J − K s ) are the total-to-selective extinction. We adopted R I,V = 1.14, R Ks,J = 0.67 for the GB sample (e.g., Dutra et al. 2002, Pietrukowicz et al. 2015 and R I,V = 1.55, R Ks,J = 0.686 for the LMC sample (e.g., Soszyński et al. 2009, Storm et al. 2011). We do not show PL diagrams for objects located in the SMC because of the low number of RG HBSs. All I-and V -band data come from the OGLE survey. In the case of MCs, values for the J and K s bands have been taken mainly from IRSF and additionally from 2MASS, and in the case of the GB, we have taken into account only data from the latter survey.
In all PL diagrams, we notice a linear trend formed by HBSs: the longer the orbital period, the more luminous the system. The second notable aspect is that HBSs seem to stick to the long-period (log P 2) group of eclipsing and ellipsoidal variables but also to the LSP stars (sequences E and D). This feature was also noticed by Soszyński et al. (2004) and highlighted by Nie et al. (2017). However, those analyses involved only a small group of HBSs located in the LMC. Here we show that those remarks are genuine for a larger sample of HBSs in the LMC, as well as for HBSs located toward the GB.
Recently, Soszyński et al. (2021) showed that LSP variables (sequence D) are systems that contain an RG and a stellar or substellar companion. In combination with eclipsing and ellipsoidal stars (sequence E), they form a group with a wide range of periods (from a few days to about 3 yr), for which the brightness changes are driven mainly by the interactions between the components of the binary system. Since the RG part of the HBSs forms a group in a similar region of the PL diagram (sequences D and E), it is likely that they are binary systems as well. Moreover, the PL diagrams show that those HBSs are a natural extension of sequence E for longer periods.
In Figure 12, the color scale of HBSs reflects the eccentricity of the system, and it is clearly seen that at a Figure 12. The PL diagrams for long-period variables (yellow, brown, pink, and red points), ellipsoidal and eclipsing binaries (green points), and HBSs (stars), located toward the LMC (panels (a) and (b)) and GB (panels (c) and (d)). On the left-hand side, the brightnesses are shown in the WI Wesenheit index, and on the other side, they are shown in WJK . The color scale of the HBSs indicates orbital eccentricity.
given brightness, the longer the period, the larger the eccentricity. In the classical ellipsoidal variables, the separation between stars in the system cannot be too small, because it may entail mass transfer via Rochelobe overflow, and it cannot be too large, because the amplitude of the brightness variations would be too slight to detect. Now, if we compare two systems with identical luminosities and periods, which indicates a similar size of the semimajor axis, but one system has a circular orbit and the other has an eccentric orbit, we can expect that the second system will be easier to detect because the minimal distance between components will be a(1−e) < a, where a is the semimajor axis; thus, the amplitude of the variations will be higher than for the first one.

SUMMARY AND CONCLUSIONS
We have presented an analysis of the HBS sample from the OCVS, described in Paper I, based on their photometric properties. The sample consists of 512 and 479 HBSs located toward the GB and MCs, respectively. An I-band light-curve variation has been modeled using a simple analytic model of the tidal deformation in the binary during periastron passage described in Kumar et al. (1995). The K95 model parameters for all HBSs were estimated using the emcee Python package, which provides the fitting with the MCMC method. We have presented the distribution of the orbital parameters eccentricity, e; inclination angle, i; and argument of the periastron, ω. The ω can be described by a flat distribution, as expected, while i can be described by the normal distribution around 65 • and a separate peak near 90 • . The majority of our HBSs have orbits with low or moderate eccentricity (0.1 e 0.35), but we also observe orbits with very high eccentricity, where e 0.5.
We have shown CMDs and H-R diagrams that indicate that HBSs, similar to the ellipsoidal and eclipsing binaries, are variables that do not belong solely to a specific evolutionary status or position on those diagrams. The most visible are two groups of HBSs. The first group of fewer than 100 systems consists of an early-type primary star lying on the MS or Hertzsprung gap, and the second group, including about 900 systems, most likely contains an RGB or AGB star. By comparing the positions of HBSs on the H-R diagram to the theoretical evolutionary tracks, we see a wide range of the initial mass of the primary, from less than 1 M to as high as 40 M . However, bear in mind that we assume a separate evolution here, which is most likely not the case for evolved RG stars.
Cool and hot HBSs differ distinctly in the period distribution and slightly in the flux variation amplitudes. The majority of hot HBSs have orbital periods in a range from a few days to 30-40 days, while the RG part of the sample is dominated by long-period binaries with a median value of about 1 yr. In the case of the I-band amplitude, the mean value is similar for both groups, but the scatter is larger for hot stars.
Like the classical ellipsoidal variables, the RG HBSs are also grouped on the PL diagrams, extending sequence E to longer periods at a given brightness, which is strong evidence confirming their binary nature. Moreover, thanks to the large number of HBSs, we proved that for a given brightness, the higher the eccentricity, the longer the period.
Using I-band photometry, we also performed timeseries analysis and found TEOs in 52 objects with a total of 76 modes. Those oscillations occur at harmonics of orbital frequencies in the range between 4 and 79.
We also provide evidence that some of them may have formed due to NLMC. Thanks to this relatively large and homogeneous sample of TEOs, we were able to construct for the first time a diagram showing the positive correlation between the TEO n and eccentricity, as predicted by theory.

ACKNOWLEDGMENT
We are grateful to the anonymous referee for many inspiring suggestions that made this manuscript more comprehensible. We thank Prof. Igor Soszyński and Prof. Andrzej Pigulski for the comments that helped to improve this manuscript. This work has been supported by Polish National Science Centre grant  Tables A1 and A2, we present the results of searching for TEOs in the OGLE HBSs.