© The Authors 2026

1. Introduction

Due to tidal distortion of the components, eccentric binaries exhibit periodic brightenings near the periastron passage. As a matter of curiosity, we note that probably the first case, when mild light brightening was observed near the periastron passage in a binary with an eccentric orbit, is the study of 96 Her (HD 164852) based on radial velocities (RVs) from photographic spectra and ground-based UBV photometry from Hvar (Koubsky et al. 1985). An attempt was made to model the effect by Hadrava (1986), and in relation to binary pulsars by Kumar et al. (1995). This class of binaries was recognised only much later, using the Kepler space telescope data by Welsh et al. (2011) and dubbed ‘heartbeat stars’ by Thompson et al. (2012). Initially, mostly low-mass stars (M ≤ 2 ℳ_⊙) were found (Shporer et al. 2016), while later a sample of massive heartbeat stars from Transiting Exoplanet Survey Satellite (TESS) was found and analysed by Kołaczek-Szymański et al. (2021).

Pejcha et al. (2013) predicted that quadruple systems in a 2+2 configuration may be very common among heartbeat stars and that their formation is more prevalent than that of triples or binaries. Although some massive heartbeat multiple systems have been observed, we are unaware of any study of a quadruple system with the 2+2 configuration with a heartbeat binary plus an eclipsing binary, as is the case in this study.

Only a handful of studies are devoted to massive heartbeat stars, for example Pablo et al. (2017) studied O9 III + B1 III/IV binary ι Ori and found multiple tidally excited oscillations (TEOs) in an O type star. MacLeod & Loeb (2023) studied the blue supergiant MACHO 80.7443.1718 with the most extreme heartbeat; Kołaczek-Szymański et al. (2022, 2024) noted that its light variation is unlikely to be caused by an extreme tidal distortion and that the ellipsoidal distortion may be of secondary importance in shaping the light curve.

HD 135160 is a bright ($V=5\stackrel{\text{m}}{.}8$) B0.5V-B1V star (also known as HR 5661, CPD−60°5698, HIP 74750, TIC 455463415 or MWC 234; α_2000.0 = 15^h16^m36$\stackrel{\text{s}}{.}$69, δ_2000.0 = −60°54′14$\stackrel{″}{.}$4) and is the brightest member of a visual binary system WDS J15166-6054A, with 8$\stackrel{\text{m}}{.}$55 component B at a separation of 1$\stackrel{″}{.}$2, and 11$\stackrel{\text{m}}{.}$62 component C at a slowly varying separation of 20″ to 27″ over more than a century. It has also been reported to be a Be star by Payne (1927). Thackeray & Emerson (1969) published two radial velocities (RVs) and argued that HD 135160 shares a common space motion with δ Cir. In their survey of the multiplicity of high-mass stars, Chini et al. (2012) classified HD 135160 as a spectroscopic binary with two visible spectra, based on the two spectrograms they had.

Our interest in this star started with an email that one of us, Mark Blackford, sent to Petr Harmanec on 15 August 2021. He was looking for a suitable comparison star for his planned photometry of δ Cir and looked at the TESS photometry of HD 135160. He found that HD 135160 is an eclipsing binary with a period of 5${}^{\mathrm{d}}.$853 and gravitational distortion of one star, and noted the presence of humps that did not correlate with the eclipse period and that small oscillations are also present throughout the orbit.

The eclipsing binary was officially reported in the paper by IJspeert et al. (2021), which was devoted to an overview of the OBA-type binaries observed by TESS. These authors identified HD 135160 as an eclipsing binary with a period 5${}^{\mathrm{d}}.$84360, having a superior conjunction at BJD 2458633.01530, with the depths of the primary and secondary minima of 0$\stackrel{\text{m}}{.}$0217 and 0$\stackrel{\text{m}}{.}$0062, respectively. Prša et al. (2022) published a catalogue of 4584 eclipsing binaries discovered in TESS data from sectors 1 to 26. For HD 135160 they quote a somewhat different value of the period, 5${}^{\mathrm{d}}.$85205(62), and the time of the first primary eclipse of BJD 2458627.15570(62). They estimated the depths of the primary and secondary minima as 0$\stackrel{\text{m}}{.}$022 and 0$\stackrel{\text{m}}{.}$004, respectively.

In the hope of obtaining an orbital solution for the eclipsing binary and determining the system properties, we started systematic spectroscopic observations of the object with the CHIRON echelle spectrograph at the Cerro Tololo Inter-American Observatory (CTIO). In this paper, we report the surprising results of our effort.

2. Observations and data reductions

Throughout this paper, we give all times of observations using the reduced heliocentric Julian dates defined by RJD = HJD − 2400000.0. We denote T_{0,  perpass} for the time of periastron passage and T_{0,  supconj} for the time of superior conjunction.

2.1. Spectroscopy

We obtained 70 high-resolution echelle spectra of HD 135160 from the CHIRON spectrograph at the CTIO over a time interval RJD 59472.4654 – 60916.4865, with a spectral resolution of R ∼ 25 000 and exposure times of 30 s, 100 s, 300 s, and 500 s in the spectral range 450–890 nm. Their overview is given in Table 1. Although a small number of echelle spectra are available in online databases, they could not be easily combined with our data, and we have decided not to use them. All spectra were converted to 1D images and wavelength-calibrated with the pipeline of the CTIO Observatory. All consecutive reductions, i.e. rectification, removal of cosmic-ray spikes and occasional flaws, and RV measurements based on the tracking paper method, were carried out with the latest version 2 of the program reSPEFO¹.

During the initial inspection of the spectra, up to three different components were found in the hydrogen and helium lines. In the two strongest He II lines at 4686 and 5411 Å, a spectral line of only one component is visible, while two or even three lines were seen in the He I lines at 4923, 5016, and 6678 Å, Hα, and Hβ. For all these lines, RVs were independently measured by Michal Zummer and Petr Harmanec, and mean values were adopted with errors being the standard deviation. They are tabulated in detail in Table 2, only available at the Strasbourg Astronomical Data Centre (CDS).

2.2. TESS photometry

Although primarily designed for the monitoring of transiting exoplanets, TESS has been invaluable in obtaining a continuous series of very accurate observations over about a month for a given area in the sky. This provided densely covered and precise light curves of eclipsing and ellipsoidal binaries, including the discoveries of new ones. We have extracted data using the lightcurve python library² (Lightkurve Collaboration 2018) from TESS sectors 12 and 65, using the Simple Aperture Photometry flux (SAP) to avoid possible removal of real astrophysical variability by detrending. To facilitate further analysis, we have normalised the light curves (LCs) of both available sectors separately to the median values.

3. Towards orbital solutions

The plots of the RVs measured versus the phase of the 5${}^{\mathrm{d}}.$853 period indicate that they do not vary with that period. This prompted us to carry out a period search using a program based on the Deeming (1975) method for amplitude periodograms, for the RVs of the He II 4686 Å line, which appeared single in all our spectrograms. The corresponding periodogram, shown in Fig. 1, clearly indicates a new period of 8${}^{\mathrm{d}}.$234 (0.1214 c d⁻¹). We also plot the RV curves based on the measurements in reSPEFO in Fig. 2. Segments of the spectra containing the Hα and He I 6678 Å lines, ordered according to the phase of the 8${}^{\mathrm{d}}.$234 period, are shown in Fig. 3 with up to 3 distinguishable stellar components.

Fig. 1. Amplitude periodogram of RVs measured on He II 4686 Å line. The highest peak at frequency of 0.1214 c d−1 corresponds to period of 8${}^{\mathrm{d}}.$234.

Fig. 2. Top: RV curve of the 8${}^{\mathrm{d}}.$234 primary measured in reSPEFO on He II 4686 Å line. Bottom: RV curves of the 8${}^{\mathrm{d}}.$234 primary (blue) and secondary (red) based on the reSPEFO RVs of the clean He I 5016 Å line. A preliminary ephemeris TRVmin. = RJD 60003.740 + 8${}^{\mathrm{d}}.$2347 × E was used in both plots. Errors are taken as a mean of two independent measurements and are too small to be seen.

Fig. 3. Parts of the CHIRON spectra with Hα and He I 6678 Å ordered according to their phase for ephemeris T0,  perpass = RJD59995.6424 + 8.23459 × E. Only some of the spectra are plotted, corresponding to around a third of the orbital period, so that change of the position of lines of 3 components (Aa, Ab, and Ba) is visible. Distances between spectra do not correspond to the distance in phase. Difference in S/N can be seen between 30 s, 300 s, and 500 s exposures.

3.1. Quadruple system

While a triple star system has only one stable configuration 2+1, quadruples are stable in 3+1 and 2+2. As two close periods of 5${}^{\mathrm{d}}.$853 and 8${}^{\mathrm{d}}.$234 are present, we have concluded from the argument of dynamical stability that HD 135160 consists of 2 binaries. A schematic illustration of the multiple system HD 135160, based on current knowledge, can be seen in Fig. 5. Only three of four stellar components are visible in the spectra, corresponding to early-type stars Aa, Ab, and Ba.

3.2. Preliminary orbital elements based on RV measurements

To obtain the input values of the orbital elements of the 8${}^{\mathrm{d}}.$234 subsystem for a more sophisticated treatment, we derived trial orbital solutions using program FOTEL (Hadrava 1990, 2004a) for the RVs, measured in reSPEFO on the He II 4686 line and the He I 5016 Å line, which is clean of neighbouring blends. These preliminary orbital elements are shown in Table 3. For the He II 4686 line, we allowed calculation of individual systematic γ velocities for each series of observation (see journal of spectra in Table 1), which resulted in a significant decrease of the rms errors. This solution is in the last column of Table 3. A plot of the local γ’s from the observing intervals versus time is shown in Fig. 4 and clearly documents the orbital motion of the Aa-Ab binary in a long orbit. Our data still do not cover one whole mutual orbit of subsystems A and B, but based on the period analysis and subsequent modelling of the long orbit, it seems probable that the value of this long period is between 1600 and 2200 days.

Fig. 4. Locally derived systemic (γ) velocities of binary A, based on the orbital solution with FOTEL from RVs measured on He II 4686 Å line.

In several spectra, we also observed faint lines that belong to the primary Ba of the eclipsing 5${}^{\mathrm{d}}.$853 subsystem, but reliable RV measurements could not be obtained in reSPEFO.

3.3. Spectra disentangling

To obtain more reliable RVs for all three spectral components, we used the program KOREL for spectra disentangling (Hadrava 1995, 1997, 2004b), and derived solutions for both subsystems A and B on a long orbit in the neighbourhood of He I 4923, 5016, 6678 Å, Hα, and Hβ lines. We then calculated the mean RV, the error being the standard deviation for components Aa and Ab from the He I lines and Hβ, and with all of the above for component Ba. Different spectral lines were chosen for components Aa+Ab and Ba after multiple different trials, and the best suitable combination was chosen. The disentangled Hα line caused greater errors in the RVs of components Aa and Ab, due to its large width, and was therefore omitted. The results for RVs can be seen in Fig. 6. We re-ran calculations multiple times, allowing the program to fit the long orbit, with results that suggest a period between 1600 and 2200 days, as predicted from the period search of γ velocities. Examples of disentangled components Aa, Ab, and Ba of the He I lines are visible in Fig. 7.

Fig. 5. Scheme of the HD 135160 system based on this study.

Fig. 6. Left: PHOEBE 2 model of primary component Aa (blue) and secondary component Ab (red) of binary A, using RVs calculated from KOREL as a mean of He I 4923, 5016, 6678 Å, and Hβ. Ephemeris of T0,  perpass = RJD59995.64 + 8.23447 × E is used. Right: PHOEBE 2 model of primary component Ba of an eclipsing binary B, using radial velocities calculated from KOREL as mean of He I lines, Hα, and Hβ. Ephemeris of T0,  supconj = RJD58627.15 + 5.852864 × E is used. Error bars are shown as a standard deviation of RVs computed from selected lines. RVs are corrected for systemic gamma velocities calculated for each segment of observations.

Fig. 7. Examples of normalised disentangled components Aa, Ab, and Ba on lines He I 4923, 5016, 6678 Å.

3.4. Estimates of the radiative properties of system components

We used the program PYTERPOL, developed by J. Nemravová (Nemravová et al. 2016; Nasseri et al. 2014), which allows the determination of radiative parameters and individual RVs of a multiple star system using the comparison of real and synthetic spectra. PYTERPOL was applied to all available blue segments of CTIO spectra, over the wavelength range 4630–5428 Å, to estimate the radiative properties and relative luminosities of the three stellar components Aa, Ab, and Ba present in the spectra. The results are summarised in Table 4, and an example of the fit for one spectrum close to the elongation of the Aa-Ab binary is shown in Fig. 8. The program also returns the RVs of the components. However, we do not use them in the following analyses, since in this particular case, three relatively similar spectra and a weak tertiary are not ideal, although they essentially confirm the preliminary results from direct measurements in reSPEFO for the Aa primary and Ab secondary. The disentangled spectra and RVs from KOREL give more accurate results for all 3 components Aa, Ab, and Ba, and the line-by-line analysis also provides some estimates of their errors for the following analyses in Sects. 4.1 and 4.2.

Fig. 8. An example of PYTERPOL fit of the spectrum from RJD 60793.672, observed near one elongation of the Aa-Ab binary.

4. Analyses of TESS light curves

An inspection of the TESS light curves of HD 135160 from sectors 12 and 65 with 120 s exposure times, with both the ellipsoidal (heartbeat) pair A in eccentric orbit and the eclipsing pair B visible (see Fig. 10), shows immediately that there are brightenings in the light curve that repeat periodically with the 8${}^{\mathrm{d}}.$234 period, and their positions correspond to instants near the periastron passages of subsystem A. This effect is now known as the ‘heartbeat effect’. Phase plots of TESS photometry from both sectors with the ephemeris of the final orbital solution for both subsystems A and B, are shown in Fig. 11.

Small periodic flux variations are observed throughout the orbit. We have carried out a period search on parts of the TESS data without the eclipses and brightenings near the periastron passages in the 8${}^{\mathrm{d}}.$234 orbit. The corresponding periodogram can be seen in Fig. 9 and shows two distinct peaks at 14.14 c d⁻¹, and 16.70 c d⁻¹, corresponding to periods of 0${}^{\mathrm{d}}.$071, and 0${}^{\mathrm{d}}.$060. However, inspection of the data plots versus time shows that the amplitude of these changes varies with time. No simple clear periodicity is obvious on the level of accuracy of the TESS data.

Fig. 9. Amplitude periodogram of TESS data from RJD 60070.22-60072.61. The highest peak at frequency of 14.14 c d−1 corresponds to period 0${}^{\mathrm{d}}.$071. The second highest peak at frequency of 16.70 c d−1 corresponds to period 0${}^{\mathrm{d}}.$060.

4.1. Eclipsing binary B

To obtain LC of eclipsing binary B, we manually removed data points related to 8${}^{\mathrm{d}}.$234 brightenings of star A, shown in Fig. 10 for both available normalised TESS sectors 12 and 65. The ‘clean’ normalised flux of binary B, phased with a period of 5${}^{\mathrm{d}}.$853, can be seen in Fig. 11 right. We have used PHOEBE 2³ (Physics of Eclipsing Binaries; Prša et al. 2016) for further modelling of the radial and light curves of the system B.

Fig. 10. Left: TESS light curve of HD 135160 from sector 12 (top) and 65 (bottom) with the 120 s exposure times using SAP flux. Primary (green) and secondary eclipses (red) of binary B with a 5${}^{\mathrm{d}}.$853 period and heartbeat effect of binary A (blue) with a 8${}^{\mathrm{d}}.$234 period are shown. Primary eclipse at RJD 58644.71 and secondary eclipse at RJD 60069.69 of binary B is reduced by heartbeat of binary A. Right: Zoomed plots of visible oscillations with a period of 0${}^{\mathrm{d}}.$071. This small variability is not always clearly visible due to changing signal noise. LCs are normalized per sector to their median values.

Fig. 11. Left: Phase folded flux of binary A after subtraction of eclipsing model for ephemeris T0,  perpass = RJD59995.6424 + 8.23459 × E. Maximum of flux occurs shortly after the periastron passage. Right: Light curve of binary B from TESS without heartbeat from binary A. Ephemeris of T0,  supconj = RJD58627.15594 + 5.8528657 × E is used. Both available sectors 12 and 65 were used in all plots and are normalized to the median value per sector.

After preliminary trials with different parameters, we used the Markov Chain Monte Carlo method (MCMC) with 10 000 steps and 95 walkers to identify possible degeneracies and probe the probability space. MCMC was re-run multiple times with slightly different initial values to be sure that we properly covered the global minimum of the solution. As for the preliminary FOTEL solutions, we derived local systemic velocities for the KOREL RVs of component Ba and pre-whitened these RVs for the long orbit for use in PHOEBE 2. Using radiative properties for the component Ba from PYTERPOL and TESS LC, we have derived the elements listed in Table 5, assuming that the orbital period P_B remains constant. We note that some of the tabulated errors might be underestimated. We used a third light fraction of 0.85, as the value computed by PYTERPOL does not converge to a reasonable solution. For star Ba, we had to use the logarithmic limb darkening law and fit its coefficients, since the interpolated limb-darkening law is not yet available for hot stars in PHOEBE 2. We note in passing that if the component Ba is a main-sequence star, the radius of Bb, which follows from the solution, is too large for a main-sequence star with T_eff ≈ 6000 K and implies that the component Bb is an evolved, giant star.

4.2. Ellipsoidal binary A with an eccentric orbit

As obtaining separate LC and fitting radii of binary A is much more difficult, we decided to first probe the orbital parameter space using the MCMC method with 10 000 steps and 95 walkers, and derived the final orbital solution for binary A based only on the radial velocities. The results of the final orbital solution are in Table 5.

A possible way to obtain the LC of the ellipsoidal binary A is to subtract the synthetic LC of binary B from the original TESS data. As data are normalised to 1, we subtracted the model from the real data and added 1. This can be seen in Fig. 11 left. It would be theoretically possible to use both RV and TESS flux to model the system for binaries A and B separately. We have assumed blackbody atmospheres and temperatures of T_Aa = 31 700 K and T_Ab = 22 000 K from PYTERPOL, and assumed that both Aa and Ab are main-sequence stars. Using estimates of radii R_Aa = 6.8 ℛ_⊙, R_Ab = 4.3 ℛ_⊙, and masses M_Aa = 17 ℳ_⊙, M_Ab = 13.6 ℳ_⊙ based on Harmanec (1988), which leads to an inclination of i = 57°, we have derived a solution in Fig. 11 left.

However, in the case of HD 135160, we point out that the shape of the brightened light curve is slightly different from the usual heartbeat. This is shown by the disparity between the data and the model in PHOEBE 2 (Fig. 11 left), which cannot explain its shape. As no reliable estimators for modelling heartbeat light curves are available, we tried multiple combinations of different parameters, changing radii, inclination, and rotational periods of stars, but none of these models came close to the observed shape of the heartbeat. Even changes in temperature did not change the shape of the light curve.

We also tried to obtain orbital parameters, including the inclination of the LC of binary A using the model by Kumar et al. (1995), with no success as the model required i ≈ 30°, which would lead to unphysical masses of stars Aa and Ab. We can only speculate that, with high-mass stars, other secondary effects are present during the periastron. Contrary to the PHOEBE 2 model, real stars have finite response times to tidal forces, which might cause a delay in maximum flux. If these estimates of masses based on main sequence stars are reliable, binary A will lead to a total mass of the system on the high end of currently discovered heartbeat stars.

5. Discussion and conclusions

We have presented a detailed study of a newly discovered massive 2+2 quadruple system HD 135160, with a unique architecture consisting of an eccentric ellipsoidal (heartbeat) binary and an eclipsing binary pair, based on our four-year systematic echelle spectroscopy and photometry from TESS. We demonstrate orbital motion of pair A on a common orbit, with a period between 1600 and 2200 days (4.4–6 years). Three components, Aa, Ab, and Ba, are massive B-type stars with temperatures of 31 700 ± 300, 22 000 ± 1000, and 16 070 ± 500 K, respectively, while the fourth star is probably of a spectral type F with a temperature of 6000 ± 500 K. Short and shallow eclipses from subsystem B are slightly more than 2% of total flux, which is explained by the fact that binary B corresponds to less than 15% of total flux from the system. The huge difference in the temperatures of Ba and Bb results in a strong reflection effect.

We also found small cyclic light variations on a time scale of 0${}^{\mathrm{d}}.$071 (14.14 c d⁻¹), with a variable amplitude. This could indicate a beat of two close periods, but we are unable to study these changes in more detail on the level of accuracy of the TESS data. We note that component Aa lies within the β Cephei instability strip where stellar pulsations are expected. Although the observed period of 0${}^{\mathrm{d}}.$071 is shorter than typical β Cep periods (0.1–0.3 d), such a short period is not entirely unprecedented (Stankov & Handler 2005). Additionally, tidal distortions in this extreme heartbeat system can, in principle, modify the stellar mode cavity and affect pulsation frequencies (Fuller et al. 2020, 2025). The nature of the observed 0${}^{\mathrm{d}}.$071 variation therefore remains unclear and warrants further investigation.

There is no evidence of Hα emission in any of our spectra; we suspect the earlier classification of the object as a Be star resulted from misinterpreting double absorption lines, which are occasionally present, as a central emission line.

As shown by the precise radial velocities obtained over four years, subsystem A is on a highly eccentric orbit, which causes periodic brightenings slightly after the periastron passage. Contrary to numerous other heartbeat stars, we are unable to obtain the inclination (and masses) of the subsystem, as the shape of the brightened peak is not easily modelled. We suggest that the brightenings may not be caused by extreme tidal distortions near the periastron, but that other secondary effects may cause a distinct shape of the light curve. We modelled both binary pairs separately, neglecting any mutual interactions such as light-travel effects or eclipse timing variations.

We hope that future spectral coverage of the entire orbit and dedicated interferometric observations might lead to the determination of precise radial velocities of components Aa, Ab, and Ba, and the orbital parameters of the long orbit. It will then be possible to investigate the evolutionary stage that leads to this remarkable configuration.

Data availability

Table 2 is available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/705/A60

Acknowledgments

This study is based on the new spectroscopic observations with the CHIRON echelle spectrograph of the Cerro Tololo Inter-American Observatory (CTIO) and on space photometry from the Transiting Exoplanet Survey Satellite (TESS) which are publicly available from the Mikulski Archive for Space Telescopes (MAST). Funding for the TESS mission is provided by NASA’s Science Mission Directorate. We gratefully acknowledge the use of the latest publicly available version of the programs FOTEL and KOREL, written by P. Hadrava, the reduction program for spectroscopic reductions reSPEFO, written by A. Harmanec, and the program for spectra modelling written by J. Nemravová. M. Zummer thanks M. Wrona for advice with this paper, A. Prša and K. Hambleton Prša for advice with PHOEBE 2 modelling. Finally, we acknowledge the use of the electronic database from the CDS, Strasbourg, and the electronic bibliography maintained by the NASA/ADS system.

References

All Tables

All Figures

	Fig. 1. Amplitude periodogram of RVs measured on He II 4686 Å line. The highest peak at frequency of 0.1214 c d−1 corresponds to period of 8${}^{\mathrm{d}}.$234.
In the text

Fig. 2. Top: RV curve of the 8${}^{\mathrm{d}}.$234 primary measured in reSPEFO on He II 4686 Å line. Bottom: RV curves of the 8${}^{\mathrm{d}}.$234 primary (blue) and secondary (red) based on the reSPEFO RVs of the clean He I 5016 Å line. A preliminary ephemeris TRVmin. = RJD 60003.740 + 8${}^{\mathrm{d}}.$2347 × E was used in both plots. Errors are taken as a mean of two independent measurements and are too small to be seen.

In the text

Fig. 3. Parts of the CHIRON spectra with Hα and He I 6678 Å ordered according to their phase for ephemeris T0,  perpass = RJD59995.6424 + 8.23459 × E. Only some of the spectra are plotted, corresponding to around a third of the orbital period, so that change of the position of lines of 3 components (Aa, Ab, and Ba) is visible. Distances between spectra do not correspond to the distance in phase. Difference in S/N can be seen between 30 s, 300 s, and 500 s exposures.

In the text

	Fig. 4. Locally derived systemic (γ) velocities of binary A, based on the orbital solution with FOTEL from RVs measured on He II 4686 Å line.
In the text

	Fig. 5. Scheme of the HD 135160 system based on this study.
In the text

Fig. 6. Left: PHOEBE 2 model of primary component Aa (blue) and secondary component Ab (red) of binary A, using RVs calculated from KOREL as a mean of He I 4923, 5016, 6678 Å, and Hβ. Ephemeris of T0,  perpass = RJD59995.64 + 8.23447 × E is used. Right: PHOEBE 2 model of primary component Ba of an eclipsing binary B, using radial velocities calculated from KOREL as mean of He I lines, Hα, and Hβ. Ephemeris of T0,  supconj = RJD58627.15 + 5.852864 × E is used. Error bars are shown as a standard deviation of RVs computed from selected lines. RVs are corrected for systemic gamma velocities calculated for each segment of observations.

In the text

	Fig. 7. Examples of normalised disentangled components Aa, Ab, and Ba on lines He I 4923, 5016, 6678 Å.
In the text

	Fig. 8. An example of PYTERPOL fit of the spectrum from RJD 60793.672, observed near one elongation of the Aa-Ab binary.
In the text

	Fig. 9. Amplitude periodogram of TESS data from RJD 60070.22-60072.61. The highest peak at frequency of 14.14 c d−1 corresponds to period 0${}^{\mathrm{d}}.$071. The second highest peak at frequency of 16.70 c d−1 corresponds to period 0${}^{\mathrm{d}}.$060.
In the text

Fig. 10. Left: TESS light curve of HD 135160 from sector 12 (top) and 65 (bottom) with the 120 s exposure times using SAP flux. Primary (green) and secondary eclipses (red) of binary B with a 5${}^{\mathrm{d}}.$853 period and heartbeat effect of binary A (blue) with a 8${}^{\mathrm{d}}.$234 period are shown. Primary eclipse at RJD 58644.71 and secondary eclipse at RJD 60069.69 of binary B is reduced by heartbeat of binary A. Right: Zoomed plots of visible oscillations with a period of 0${}^{\mathrm{d}}.$071. This small variability is not always clearly visible due to changing signal noise. LCs are normalized per sector to their median values.

In the text

Fig. 11. Left: Phase folded flux of binary A after subtraction of eclipsing model for ephemeris T0,  perpass = RJD59995.6424 + 8.23459 × E. Maximum of flux occurs shortly after the periastron passage. Right: Light curve of binary B from TESS without heartbeat from binary A. Ephemeris of T0,  supconj = RJD58627.15594 + 5.8528657 × E is used. Both available sectors 12 and 65 were used in all plots and are normalized to the median value per sector.

In the text