© The Authors 2023

1. Introduction

For many reasons, massive stars (≳8 M_⊙) are of particular interest to modern astrophysics. Primarily, they are progenitors of core-collapse supernovae (e.g. Janka et al. 2007; Smartt 2009) and long γ-ray bursts (e.g. Fruchter et al. 2006; Yoon et al. 2006). For billions of years, they contributed to the chemical evolution of the entire Universe and interacted mechanically with the surrounding interstellar medium (e.g. Ouellette et al. 2007; Svirski et al. 2012), also through their intense line-driven stellar winds (Vink 2022). Furthermore, most of their remnants are compact objects, such as neutron stars (NSs; including magnetars and pulsars) and black holes (BHs), which allow the empirical study of effects of general relativity. Finally, massive stars can be observed at cosmological distances due to their enormous luminosities, and hence they dominate in the spectra of distant starburst galaxies (see Eldridge & Stanway 2022, for a recent review). These features of massive stars (and many others) demonstrate that understanding the structure and evolution of massive stars is one of the key tasks of astronomy.

As is well known, massive and intermediate-mass (≳2 M_⊙ and ≲8 M_⊙) stars reside in binary systems at a much higher rate than their lower-mass counterparts (Duchêne & Kraus 2013; Sana et al. 2012). Moreover, as shown, for example, by Sana et al. (2014) and Moe & Di Stefano (2017), O-type dwarfs in particular are often found in multiple systems. This shows that binarity is inherent in the evolution of massive stars and cannot be ignored when studying these objects or their final outcomes. Many interesting phenomena in the Universe are the result of binarity among massive and/or intermediate-mass stars. These include Be stars (Kriz & Harmanec 1975; Bodensteiner et al. 2020), so-called stripped stars (Götberg et al. 2020; Shenar et al. 2020; El-Badry & Burdge 2022), BH–BH, BH–NS, and NS–NS mergers (progenitors of gravitational wave events; Abbott et al. 2016, 2019), ‘early’ stellar mergers (Tokovinin & Moe 2020; Sen et al. 2022; Li et al. 2022), Ib and Ic supernova progenitors (Langer 2012; Yoon et al. 2012), ‘massive Algols’ (de Mink et al. 2007; Skowron et al. 2017; Sen et al. 2022), and even Wolf–Rayet stars (Shenar et al. 2019; Pauli et al. 2022).

At the same time, binary systems that contain massive main-sequence (MS) component(s) are often characterised by eccentric orbits due to their relatively young age and the presence of radiative outer layers, which are less vulnerable to tidal dissipation compared to convective envelopes (e.g. Van Eylen et al. 2016). Both observational studies of large samples of massive binaries (e.g. Moe & Di Stefano 2017) and hydrodynamical simulations of their formation (see Oliva & Kuiper 2020, and references therein) suggest significantly non-zero eccentricities at their birth.

Assuming that the periastron distance between the components is sufficiently small¹, the combined proximity effects, such as ellipsoidal distortion, the irradiation effect, and Doppler beaming and boosting, make such a system an eccentric ellipsoidal variable (EEV; e.g. Nicholls & Wood 2012). Due to the characteristic shape of the light curve of EEVs during the periastron passage (which can resemble an electrocardiogram), EEVs are sometimes referred to as ‘heartbeat stars’ (Welsh et al. 2011; Thompson et al. 2012; Beck et al. 2014; Kirk et al. 2016; Kołaczek-Szymański et al. 2021; Wrona et al. 2022a).

The orbital phase-dependent tidal potential acting on the components of EEVs can induce tidally excited oscillations (TEOs) in their interiors (Zahn 1975; Kumar et al. 1995; Fuller 2017; Guo 2021), which in turn can affect the evolution of the binary system. However, many details of TEOs in massive and intermediate-mass stars are still poorly understood, including the total number of TEOs and their frequency of occurrence. In our study, we aim to shed light on this issue based on theoretical modelling combined with machine learning (ML) techniques.

The paper is organised as follows. Section 2 provides a concise characterisation of TEOs and specifies the purpose of our paper. In Sect. 3 we present a detailed description of the adopted methodology, including the assumptions made and the software used to generate the theoretical models. We then analyse the obtained models and present our findings in Sect. 4. Finally, we summarise the entire work and draw several conclusions in Sect. 5.

2. Properties of TEOs and the purpose of the paper

Tidally excited oscillations are tidally forced eigenmodes of a star with frequencies, σ_nlm (in the co-rotating frame of the star), coinciding with integer multiples, N, of the orbital frequency, f_orb². The resonance condition can be written as follows:

$$\begin{array}{c}\hfill f_{\mathit{Nm}}\equiv N{f}_{\mathrm{orb}}-m{f}_{\mathrm{s}}\approx {\sigma }_{\mathit{nlm}},\end{array}$$(1)

where f_Nm corresponds to the frequency of the tidal forcing in the rotating frame, f_s stands for the rotational frequency of the star, while the subscripts n, l and m denote the radial order, degree, and azimuthal order of the specific eigenmode, respectively. This property of TEOs makes them relatively easy to distinguish from other types of pulsations (e.g. self-excited oscillations) in frequency spectra, provided the orbital period is known. There are numerous examples of photometrically or spectroscopically detected TEOs (e.g. Handler et al. 2002; Welsh et al. 2011; Hambleton et al. 2013; Fuller et al. 2017; Guo et al. 2019, 2020; Wrona et al. 2022b), also in massive binary systems (e.g. Willems & Aerts 2002; Pablo et al. 2017; Kołaczek-Szymański et al. 2021, 2022). Most TEOs are damped normal modes, meaning that without constant tidal forcing they would not be observed in the star. More importantly, because of their damped nature, TEOs dissipate the total orbital energy making the system tighter, more circular, and synchronised with time. On the other hand, if the TEO is naturally an over-stable mode it can transfer thermal energy from the star to the orbit via so-called inverse tides (Fuller 2021). Regardless of the type of TEOs, they undeniably contribute to the evolution of the (massive) binary system, and can therefore influence the characteristics of the phenomena and objects mentioned above. The efficiency of energy transfer between the stellar interior and the orbit due to TEOs strongly depends on the temporal behaviour of the resonance condition given by Eq. (1). It is to be expected that most TEOs are ‘chance resonances’, that is to say, resonances in which the aforementioned condition is satisfied for a relatively short time. Under such circumstances, TEOs do not have enough time to reach high amplitudes, and hence their ability to dissipate orbital energy is somewhat limited. However, if, after reaching resonance, both frequencies on the left and right sides of Eq. (1) evolve at the same rate and direction, TEOs can ‘tidally lock’ for a longer time compared to the chance resonance scenario (Fuller 2017). This unique variety of TEOs is suspected to be responsible for occasional periods of rapid evolution of the orbital parameters in binary systems (Fuller et al. 2017).

Tidally excited oscillations are not only restricted to MS stars, they can also occur in binaries with pre-MS stars (Zanazzi & Wu 2021), some compact objects (white dwarfs, Yu et al. 2021), planetary systems (Ma & Fuller 2021) and even planet-moon systems (e.g. in the Saturn-Titan system, Lainey et al. 2020).

Although the literature on theoretical studies of TEOs is indeed extensive (see e.g. Fuller 2017; Guo 2021, for recent reviews), the question of their rate of occurrence and the role they play in the evolution of massive stars is still a matter of debate. Unfortunately, only a small number of papers refer exclusively to massive EEVs. Witte & Savonije (1999a,b) studied gravity (g) and Rossby-mode TEOs in an uniformly rotating 10 M_⊙ MS star using their own two-dimensional (2D) code for different stellar rotation rates and several orbital configurations. They found that dynamical tides can effectively circularise and tighten the orbits of EEVs in just a few megayears if resonance locking occurs. However, these and many other previous works on TEOs were done under the assumption of a compact (point-like) secondary companion that is not subject to tidal perturbations during each periastron passage. This is obviously not the case in real binary systems, where both components are responsible for the tidal evolution of the orbit. As theoretically shown by Witte & Savonije (2001), for an eccentric binary system consisting of two 10 M_⊙ stars, tidal dissipation can be further enhanced due to the simultaneous excitation of tidally locked TEOs in both components. In spite of the advanced mathematical formalism, the aforementioned papers only dealt with a few assumed component masses and sets of orbital parameters. Only Willems (2003) attempted a qualitative analysis of the hyperspace of the orbital parameters favouring excitation of TEOs in massive EEVs on MS. He found that for a mass range of 2−20 M_⊙, the favourable orbital period interval lies between ∼5 and ∼12 d when both components are zero-age main-sequence (ZAMS) stars. This interval shifts towards longer orbital periods (up to ∼70 d) for components approaching the terminal-age main sequence (TAMS).

Although, as argued above, the role of TEOs in the life of massive binary systems is still not well understood, we are not aware of any published work that develops the qualitative analysis carried out by Willems (2003) based on state-of-the-art stellar evolution and oscillations codes. We would like to fill this gap by combining the Modules for Experiments in Stellar Astrophysics³ (MESA; Paxton et al. 2011, 2013, 2015; Paxton et al. 2018, 2019) stellar structure and evolution code with the GYRE⁴ (Townsend & Teitler 2013; Townsend et al. 2018) non-adiabatic stellar oscillation code. Our study aims to answer the following three questions: The first is how many resonances (given by Eq. (1)) can EEVs experience during their lifetime between the ZAMS and the TAMS, and how this picture changes with different initial parameters of the binary system. The second question regards whether we can distinguish several distinct types of EEV resonance histories that are statistically related to the initial physical and orbital parameters of binary systems. Finally, we want to determine if the resonance history correlates in any way with the properties of EEVs near the TAMS, for instance, whether systems that undergo mass transfer before reaching the TAMS are also systems with a higher total number of encountered resonances. In order to answer the last two questions, we used of ML techniques by performing a Uniform Manifold Approximation and Projection⁵ (UMAP; McInnes et al. 2018) dimension reduction analysis of the resonance histories obtained for simulated binary systems.

3. Methods

Assuming that the dynamical tide excited in the component is dominated by a single TEO close to resonance with the orbital harmonic N, one can express its amplitude of luminosity change as proportional to (after Fuller 2017, his Eq. (2))

$$\begin{array}{c}\hfill A_N\propto q{\left(\frac{R}{a}\right)}^{l+1}|Q_{\mathit{nlm}}|F_{\mathit{Nm}}{\mathcal{L}}_N,\end{array}$$(2)

where q is the ratio of the masses of the two components, R stands for the radius of the component in which the TEO is excited, and a is the semi-major axis of the relative orbit. The quantity denoted Q_nlm is known as the so-called overlap integral, which describes the spatial coupling between a given oscillation mode and the actual geometry of the tidal potential (Fuller 2017, his Eq. (4)). In general, the larger the value of |n|⁶, the smaller the Q_nlm, and hence eigenmodes with a large number of nodes in the radial direction have a much lower probability of tidal excitation⁷. In addition to Q_nlm, the Hansen coefficient F_Nm (Fuller 2017, his Eq. (5)) is responsible for the temporal coupling of the forced normal mode and the Nth component in the Fourier expansion of the orbital motion. Quantitatively, it expresses the intuitive principle that for more eccentric orbits, TEOs with larger orbital harmonic numbers will be excited. This is because, as the eccentricity increases, the periastron passage takes less time for a given orbital period, so eigenmodes with higher frequencies better ‘match’ rapidly changing gravitational field, in terms of timescale. Nevertheless, for very high N, the F_Nm drops rapidly (almost exponentially). This particular property of F_Nm is responsible for the lack of excitation of the TEOs with extremely high N. It is clear here that the frequency range of TEOs in massive and intermediate-mass MS stars is limited on two sides independently by Q_nlm and F_Nm. On the low-frequency side, Q_nlm prevents the excitation of g modes with very high |n|, while on the high-frequency side F_Nm decreases sharply, strongly limiting the possible excitation of pressure (p) modes characterised by high radial orders. The last term in Eq. (2), ℒ_N, denotes the de-tuning factor given by the following formula,

$$\begin{array}{c}\hfill {\mathcal{L}}_N=\frac{f_{\mathit{Nm}}}{\sqrt{{({\sigma }_{\mathit{nlm}}-f_{\mathit{Nm}})}^2+{\gamma }_{\mathit{nlm}}^2}},\end{array}$$(3)

where γ_nlm stands for damping or growth rate of the normal mode. This Lorentzian-like factor reflects the mismatch between f_Nm and σ_nlm. Given the typical values of |γ_nlm| for g and p modes in massive and intermediate-mass MS stars (of the order of ∼10⁻⁷ − 10⁻³ d⁻¹), ℒ_N is extremely sensitive to the difference (σ_nlm − f_Nm). Hence, among many other factors, the ℒ_N undoubtedly plays a key role in the excitation of TEOs.

While the precise prediction of TEO amplitude is a difficult task⁸, we are interested in analysing the changes in resonance conditions dictated by the sum of all the contributing de-tuning factors with passing time, t. We define the following dimensionless quantity,

$$\begin{array}{c}\hfill \mathcal{L}(t)\equiv {\displaystyle \sum _{\mathit{nlm}}\sum _{N=1}^{N_{\max }}}{\mathcal{L}}_N(t).\end{array}$$(4)

In contrast to ℒ_N, associated with a single orbital harmonic, ℒ reflects the overall chance of TEOs being excited in the EEV component. However, we must stress at this point that it does not carry direct information on the amplitude of potential TEOs. The first summation in Eq. (4) applies to all the normal modes we consider in the modelling (Sect. 3.3). Obviously, the second summation in Eq. (4) should run from N = 1 to +∞, but due to time and physical constraints one has to truncate the series at some reasonably chosen N_max. A detailed description of the selection of N_max values is provided in Sect. 3.3.

In order to try to answer the questions raised in Sect. 2, we synthesised 20 000 binary evolution models and calculated ℒ(t) for both components in each of them. The whole procedure is described extensively in the next four subsections.

3.1. General assumptions

From a practical point of view, a fully consistent calculation of the evolution of binary systems taking TEOs into account is very time-consuming, as it requires time steps shorter than the times at which the resonances occur (several orders of magnitude shorter than the nuclear timescale, cf. Fig. 3). It would take an enormous amount of time to perform such consistent calculations for 20 000 binaries with hundreds of resonances occurring in each of them. Therefore, to make our project both feasible and still scientifically useful, the models were synthesised under the general assumption that each resonance encountered by the EEV components is a chance resonance. By sacrificing the ability to track resonantly locked TEOs, we are able to decouple evolutionary and seismic calculations and run them independently, greatly simplifying the whole problem. We believe that we can to do this for three reasons: (1) we are only interested in obtaining some general statistical information about the resonance conditions in a large number of simulated binary systems, (2) the phenomenon of resonance locking is rare compared to the rate of chance-resonance events, and (3) the impact of chance-resonance TEOs on the orbit is limited due to their relatively short time of existence (e.g. Witte & Savonije 1999b). In conclusion, we focus on finding candidate binaries that may or may not experience numerous TEOs, rather than precisely predicting their actual evolutionary histories, which is beyond the scope of this paper. We believe that our results will serve as a starting point for more detailed calculations performed for the most interesting cases of massive EEVs.

3.2. Synthesis of binary evolution models

Since we assumed that we could separate stellar and orbital evolution from seismic calculations, we first generated a set of binary evolutionary tracks and only then performed seismic analysis on them to find ℒ(t).

3.2.1. Initialisation of models

We used the latest open-source 1D stellar evolution code MESA (release 15140) compiled with the MESA Software Development Kit (version 21.4.1, Townsend 2021) to compute a set of 20 000 binary evolution models. The MESAbinary module (Paxton et al. 2015) allows the simultaneous evolution of binary system components and their orbits. Throughout this paper, we use the subscripts ‘1’ and ‘2’ to denote the primary (initially more massive) and secondary components, respectively.

We assumed that both components have the same chemical composition with metallicity Z = 0.02 and a solar-scaled mixture of elements taken from Grevesse & Sauval (1998). Since we were only interested in massive and intermediate-mass MS EEVs that can exhibit TEOs during their lifetime, the initial systems consisted of two stars lying on the ZAMS and were characterised by the following seven parameters, randomly drawn from the following uniform distributions, 𝒰_[α, β], on specific intervals [α, β].

The first is the mass of primary component, log(M₁/M_⊙)∼𝒰_{[log5, log30]}. A uniform distribution on a logarithmic scale was used instead of a linear scale to cover the Hertzsprung–Russell diagram (HRD) with more evenly distributed evolutionary tracks.

Second is the mass ratio, q ≡ M₂/M₁ ∼ 𝒰_{[0.2, 0.95]}, where M₂ corresponds to the mass of the secondary component. The lower limit for q was introduced due to the fact that the efficiency in driving TEOs scales with q (cf. Eq. (2)), so it is less likely to observe TEOs in a binary system at a small value of the mass ratio. Moreover, if the generated q corresponded to M₂ <  2 M_⊙, a redraw was performed.

Third is the eccentricity, e ∼ 𝒰_[0.2, 0.8], the typical range of observed EEVs. Fourth is the periastron distance, r_peri, normalised to the sum of components’ radii, ${\stackrel{\sim }{r}}_{\mathrm{peri}}\equiv r_{\mathrm{peri}}/(R_1+R_2)\sim {\mathcal{U}}_{[1,5.5]}$. However, if the generated system was initially Roche-lobe overflowing (RLOF) at the periastron, a redraw was performed. We also assumed an upper value of ${\stackrel{\sim }{r}}_{\mathrm{peri}}$ because the overall strength of tidal forces decays as $r_{\mathrm{peri}}^{-3}$ and simulating widely separated systems would contradict the aim of this paper.

Fifth is the tidally enhanced wind factor, B_wind ∼ 𝒰_[32 896]. Introduced by Tout & Eggleton (1988) for red giants residing in binary systems, it accounts for the tidal enhancement of the stellar wind mass-loss rate due to the presence of a nearby companion. The ad hoc chosen range of B_wind corresponds to a maximum amplification of the ‘nominal’ wind mass-loss rate by a factor of 1.5−10 (cf. Tout & Eggleton 1988, their Eq. (2)).

Sixth is the angular rotational velocity divided by its critical value⁹, Ω/Ω_crit ∼ 𝒰_[0.1, 0.5]. The assumed range of initial Ω/Ω_crit translates into the linear equatorial velocities between ∼50 km s⁻¹ and ∼320 km s⁻¹ in our simulations and reflects the significant non-zero rotation velocities observed in massive young MS stars (e.g. Dufton et al. 2006; Hunter et al. 2008).

The final one is the overshoot mixing parameter, f_ov ∼ 𝒰_{[0.015, 0.025]}. In our calculations, the overshooting of the material above the convective, hydrogen-burning core was treated in the exponential diffusion formalism developed by Herwig (2000). An adjustable parameter, f_ov, describes the spatial extent of the overshoot layer in terms of the local pressure-scale height, but its value for massive stars is still under debate. We adopted the range of f_ov after Ostrowski et al. (2017).

The parameters presented above were generated independently for each EEV system. Moreover, the last two parameters, Ω/Ω_crit and f_ov, were drawn independently for each component, so the final hyperspace of parameters explored in our simulations included $\{M_1,q,e,{\stackrel{\sim }{r}}_{\mathrm{peri}},\mathrm{\Omega }/{\mathrm{\Omega }}_{\mathrm{crit},1},\mathrm{\Omega }/{\mathrm{\Omega }}_{\mathrm{crit},2},f_{\mathrm{ov},1},f_{\mathrm{ov},2},B_{\mathrm{wind}}\}$. Figure 1 shows our initial sample of generated EEVs in the orbital period versus eccentricity diagram. As expected, they occupy the upper envelope of the aforementioned plane with the upper boundary dictated by the onset of periastron RLOF on the ZAMS. The rest of the necessary parameters and ‘physics switches’ were identical for each simulated binary. We briefly describe them below.

3.2.2. Integration of the evolution

Nuclear reaction rates were calculated using ‘basic.net’ option in MESA. We used a convective premixing scheme (Paxton et al. 2019, their Sect. 5) in combination with the Ledoux criterion to define the boundaries of convective instability. This specific approach of treating convection agrees with the results of modern 3D hydrodynamic simulations (Anders et al. 2022). Convective mixing was incorporated into the models via mixing length theory (MLT) in the ‘Cox’ formalism (Cox & Giuli 1968, their Chap. 14) with the value of the solar-calibrated mixing length parameter α_MLT = 1.82¹⁰ (Choi et al. 2016). As mentioned earlier, exponential overshoot mixing above the convective core was also included¹¹, but we neglected overshooting in the non-burning convection zones. For stars with masses ≥15 M_⊙, we activated the treatment of convection as ‘MLT++’ (Paxton et al. 2013, their Sect. 7.2) to reduce super-adiabaticity in convective zones dominated by radiation pressure. Since we used Ledoux criterion, semi-convection could appear in our stars with its efficiency parameter, α_sc = 0.01 (Langer et al. 1985). In our case, semi-convection sometimes occurred in chemically modified layers left by the shrinking core.

Upon initialisation at the ZAMS, we relaxed both components in ∼100 steps so that they rotated rigidly. Later, we allowed our stars to rotate differentially during their evolution, according to the so-called shellular approximation of rotation (Meynet & Maeder 1997). Throughout the entire evolution, we assumed that the rotation axes of the stars are perpendicular to the orbital plane. MESAstar uses the mathematical formalism of Heger et al. (2000, 2005) to apply structural corrections, perform different types of rotationally induced mixing and ‘diffusion’ of angular momentum (AM) between adjacent shells. The following rotational mixing mechanisms were taken into account in MESA: dynamic shear instability, secular shear instability, Eddington–Sweet circulation, Solberg–Høiland instability, and Goldreich–Schubert–Fricke instability (all described in detail by Heger et al. 2000). Even the combined mixing coefficients of the aforementioned rotational instabilities can be zero in some parts of the star. However, this is clearly unrealistic due to the presence of a nearby companion that induces additional mixing throughout the star. To at least approximately account for this process, we did not allow the total mixing coefficient, D_mix, to fall below 10⁵ cm² s⁻¹. This particular arbitrarily selected value is related to the mixing timescale, τ_mix ∼ (Δr)²/D_mix ≈ 15 Myr at radial distance, Δr = 0.1 R_⊙. We cannot conceal here that rotation and mixing profiles in MS stars are still poorly understood (except in the solar case). There are no definitive conclusions as to what mixing mechanisms and whether they actually occur in massive and intermediate-mass MS stars (see e.g. Pedersen et al. 2021; Pedersen 2022, or a discussion of this problem and its asteroseismic inference from B-type dwarfs).

Mass losses due to the radiation-driven stellar wind were calculated according to the prescription given by Vink et al. (2001). Their formulae take into account the presence of a so-called bi-stability jump around the effective temperature of T_eff ≈ 26 000 K, caused by ionisation and recombination of some Fe ions. Nevertheless, the presence of a bi-stability jump is still questionable and there is some evidence that the associated almost instantaneous change in the mass-loss rate may not be real (cf. Krtička et al. 2021; Björklund et al. 2023). ‘Nominal’ wind mass-loss rates in our simulations were modified in two ways: (1) the rate was amplified by the aforementioned tidal mechanism, parametrised by B_wind (Tout & Eggleton 1988) and (2) the effect of fast rotation at the surface, which can amplify the mass-loss rate, was accounted for by the simplified power-law description given by Heger et al. (2000; their Sect. 2.6). We assumed that the mass loss through the wind is completely non-conservative, meaning that there is no mechanism that could transfer some material back to the star or to a companion.

As we already mentioned above, MESAbinary allows the simultaneous integration of some stellar and orbital parameters that are coupled to each other in a binary system. We switched on the MESA controls responsible for changes in the total orbital AM caused by: (1) gravitational wave radiation, (2) wind mass loss and (3) tidal spin-orbit coupling. For the first process, the rate of orbital momentum loss was calculated assuming point masses. The mass loss through the stellar wind was completely non-conservative, so the AM lost via this channel was equal to the AM carried by the wind. The phenomenon (3), contributing to the evolution of eccentricity, orbital and spin angular momenta, was modelled using the theory of tidal interactions for radiative envelopes developed by Zahn (1977), Hut (1981) and Hurley et al. (2002), after being adapted to the shellular approximation of rotation. For stars with radiative envelopes, tidal dissipation processes are dominated by tidally excited gravity modes that propagate to the stellar surface, where they gain relatively large amplitudes and experience effective radiative damping (due to the short local thermal timescale) and non-linear damping. Consequently, they deposit their energy and AM in the outer layers of the envelope. Following earlier calculations of Zahn (1977), Hurley et al. (2002) delivered convenient formulae to describe the tidal synchronisation and circularisation timescales associated with the aforementioned phenomenon. Using these timescales combined with the formalism presented by Hut (1981), MESAbinary integrates the evolution of the eccentricity and updates spin angular frequency of each shell in the stellar model. Therefore, our calculations in MESAbinary take into account the approximate influence of the dynamical tide on the orbit, at least up to the lowest possible order. Of course, the tidal evolution formalism implemented in MESAbinary does not include the effects of resonance locking. For explicit formulae describing tidal processes in MESAbinary, we refer to Paxton et al. (2015; their Sect. 2).

We completely ignored the effects of magnetic fields, while bearing in mind that they may mainly affect the actual stellar wind mass-loss rates, the efficiency of internal mixing processes and synchronisation or circularisation timescales (e.g. via the magnetic braking mechanism). The impact of fossil magnetic fields on the evolution of massive and intermediate-mass stars was recently described by Keszthelyi et al. (2022).

All details on the parameters of our models in MESA can be found in Appendix A, where we present the contents of our MESA inlists. A concise description of the micro- and macrophysics data sources used by MESA is provided in Appendix C.

3.2.3. Termination conditions

The evolution of the binary system was carried out until at least one of the following termination conditions was met for any of the components: (1) the component reached the TAMS, in other words, the central mass abundance of hydrogen fell below X_c ≤ 10⁻⁴; (2) the eccentricity was reduced to e ≤ 0.01; (3) the rotation velocity reached Ω/Ω_crit = 0.75 at the stellar surface; or (4) episodic mass transfer between components due to the RLOF in the periastron began.

The reasons behind providing the termination conditions outlined above are as follows. Our study is exclusively dedicated to the MS phase of the evolution of EEVs, and hence the first condition has to be enforced. The second condition is self-explanatory, since we are interested in non-zero eccentricities that allow for TEO excitation¹². The third condition is related to the convergence problems that can occur in MESAbinary when one of the components nearly approaches the break-up velocity of rotation. Numerous assumptions and descriptions of rotation-related phenomena reach the limits of their applicability in MESA for Ω/Ω_crit ≈ 1. Since for Ω/Ω_crit ≳ 0.75 the deviation from spherical symmetry becomes significant, a 1D treatment of the problem is no longer adequate. For instance, the way in which such a star loses mass becomes fundamentally different from the isotropic case. We therefore decided to stop integrations under such circumstances. The last condition is related to the difficulty in correctly describing an episodic (near-periastron) RLOF, when a ‘blob’ of material could be ejected from the Roche-lobe-filling component during each periastron passage. However, this kind of orbital phase-dependent RLOF is not expected to be observed in a binary for a long time due to strong tidal forces. They should effectively suppress the eccentricity, making the system circular (and so the second condition can be quickly met).

3.3. Asteroseismic calculations

A consequence of the assumption of the aligned vectors of the orbital and spin angular momenta is a rule for selecting the geometry of modes that can be tidally excited. Under such conditions, a normal mode can be tidally excited only if mod(|l + m|,2)=0; for example, the l = 2 TEOs will be characterised only by m = −2, 0, +2. Here we restrict our study to only l = 2, m = 0, +2 modes because of two reasons. First, l = 2 modes correspond to the dominant component in the series expansion of the variable tidal potential. Modes with l >  2 undergo much weaker excitation due to the dependence on (R/a)^l + 1, which enters Eq. (2). Second, the values of F_N, −2 are very small compared to their m = 0, +2 counterparts. This can be easily seen in Fig. 2a, where we have plotted the maximum values of F_Nm for m = −2, 0, +2 and different eccentricities. F_N, −2 is approximately at least 2−3 orders of magnitude smaller than F_N, 0 or F_N, 2.

For each model of the stellar internal structure that was saved during the synthesis of binaries in MESA, we calculated the oscillation spectrum using the GYRE code in the non-adiabatic regime and the second-order Gauss–Legendre Magnus integrator. The frequencies σ_n, 2, 0 and σ_n, 2, +2 corresponding to the non-adiabatic calculations were searched by GYRE based on the preliminary adiabatic calculations. Rotational effects (Coriolis force) were taken into account using the so-called traditional approximation of rotation (e.g. Aerts et al. 2010, their Sect. 3.8). We searched for eigenvalues in the family of gravito-acoustic solutions. We assumed the necessary (differential) rotation profile inside the star from the MESA model.

As we argue in Sect. 3, it is necessary to choose a reasonable range of frequencies to scan for eigenvalues based on Q_nlm and F_Nm. Therefore, we only searched for modes with |n|≤30 and frequencies, ${\sigma }_{n,2,0}\in (f_{\mathrm{orb}},N_{\max }^{m=0}{f}_{\mathrm{orb}})$ or ${\sigma }_{n,2,+2}\in (max\{0,f_{\mathrm{orb}}-2f_{\mathrm{s},\mathrm{core}}\},N_{\max }^{m=+2}{f}_{\mathrm{orb}}-2f_{\mathrm{s},\mathrm{core}})$. In the ranges shown, $N_{\max }^{m=0}$ and $N_{\max }^{m=+2}$ refer to the limits of N due to the decrease in F_N, 0 and F_N, 2, respectively. The f_s, core is the core rotation frequency. We defined $N_{\max }^{m=0}$ and $N_{\max }^{m=+2}$ as N for which F_N, 0 or F_N, 2 is equal to 10⁻⁸ (i.e. F_Nm starts to effectively prevent the excitation of TEOs). In practice, we numerically calculated the F_Nm functions¹³ for different eccentricities and obtained the $\log N_{\mathrm{m}\mathrm{a}\mathrm{x}}^m(e)$ relations as a fit of a fourth-degree polynomial to a set of its discrete points. A summary of this process is shown in Fig. 2b. For low-e orbits, the typical range of N favourable for the excitation of TEOs reaches N ∼ 10¹, in contrast to highly eccentric orbits, which may exhibit as much as N ≈ 100 − 200 TEOs. Figure 2b also shows another feature of m = −2 modes that makes them inferior candidates for TEOs compared to axisymmetric and prograde modes – as potential TEOs they always span a narrower range of orbital harmonic numbers.

Defining the frequency range for σ_n, 2, 0 is quite straightforward, as these are axisymmetric modes and their frequencies do not change when switching between inertial and co-rotating frames. The situation is quite different when it comes to the m = +2 modes. This time, due to the differential rotation inside the star, ${\sigma }_{\mathit{nlm}}={\sigma }_{\mathit{nlm}}(r)={\overline{\sigma }}_{\mathit{nlm}}-m{f}_{\mathrm{s}}(r)$, where r is the radial coordinate in the stellar interior and ${\overline{\sigma }}_{\mathit{nlm}}$ is oscillation frequency in the inertial frame. For some eigenmodes, σ_nlm may change its sign somewhere in the star, depending on the shape of the rotational profile. This location is known as the critical layer, where σ_nlm(r)→0, and such a mode experiences severe damping due to its very short radial wavelength (e.g. Alvan et al. 2013). To exclude these modes from our experiment, the maximum frequency of σ_n, 2, +2 was set to $(N_{\max }^{m=+2}{f}_{\mathrm{orb}}-2f_{\mathrm{s},\mathrm{core}})$ ¹⁴. This is because during evolution the core rotates almost rigidly and faster than the envelope, hence the difference $(N_{\max }^{m=+2}{f}_{\mathrm{orb}}-2f_{\mathrm{s},\mathrm{core}})\le (N_{\max }^{m=+2}{f}_{\mathrm{orb}}-2f_{\mathrm{s},\mathrm{env}})$, where f_s, env stands for the rotation frequency of the outermost part of the envelope.

More details of our calculations performed in GYRE can be found in Appendix B, where we present the explicit contents of our GYRE input file.

3.4. Derivation of ℒ(t)

The introduction of differential rotation also has consequences when it comes to interpreting the resonance condition from Eq. (1). The quantity f_s is no longer a constant value, so one has to decide which f_s to choose. Theoretical studies imply the induction of g-mode TEOs (especially of high radial order) primarily near the convective core boundary (e.g. Goldreich & Nicholson 1989) for stars with radiative envelopes. However, the resonances in our simulations are also due to p or g modes with small radial order. Therefore, we decided to apply our resonance condition to the envelope¹⁵ (not to the interface region near the core boundary), rewriting Eq. (1) more accurately as

$$\begin{array}{c}\hfill f_{\mathit{Nm}}\equiv N{f}_{\mathrm{orb}}-m{f}_{\mathrm{s},\mathrm{env}}\approx {\sigma }_{\mathit{nlm}},\end{array}$$(5)

and use it in the subsequent modelling of ℒ(t). It is essential to note at this point that the resonance condition given by Eq. (5) refers to f_Nm and σ_nlm expressed in a frame co-rotating with the outer stellar envelope. Although in principle the morphology of ℒ(t) depends on the choice of the specific resonance condition, we note that it does not affect at all resonances due to m = 0 modes and should not significantly affect resonances corresponding to p modes or low-|n|, m = +2 g modes.

Having a set of eigenfrequencies calculated by GYRE and knowing the history of the binary evolution from MESA, we performed the summation shown in Eq. (4). However, this was not a direct summation running over the models saved by MESA and GYRE, as their temporal resolution was still too coarse compared to the duration of a typical resonance. To circumvent this problem, we interpolated the temporal variations of each oscillation frequency and all necessary parameters of the binary system using Akima cubic spline functions (Akima 1970). Then, we were able to calculate the values of ℒ(t) on a uniformly spaced time grid with a constant time step of 2000 yr, which we assumed to be identical for all EEVs in our simulations. From here on, we use the term ‘resonance curve’ as a proxy for the ℒ(t) time series. Figure 3 shows a compilation of example resonance curves, although we postpone discussion of these to Sect. 4. Together with the initial parameters of binary systems, resonance curves are particularly important to us in this study.

3.5. ML analysis of the resonance curves

Although in Sects. 4.3 and 4.4 we analyse the resonance curves based on various statistics, due to their global nature we do not distinguish many details that are ‘hidden’ in the resonance curves. To characterise the morphology of all resonance curves in more detail (without having to perform a visual classification, which is almost impossible due to the number and complexity of the dataset), we applied dimensionality reduction methods. With these, we were able to explore the topology spanned by the morphological features of the resonance curves. We carried out the entire analysis described here separately for the sets of curves ℒ₁(t) and ℒ₂(t), corresponding to the primary and secondary components, respectively.

As a first step, we summarised each resonance curve with a vector Q that described its morphological features. We focused our attention on two particular features: (1) the distribution of log(ℒ) values and (2) the distribution of apparent maxima at a normalised time, t/T_max, where T_max stands for the max{t}. In practice, we calculated vectors Q_x and Q_y, which contained sets of 1000 quantiles of normalised times corresponding to local maxima of ℒ(t) and 1000 quantiles of log(ℒ), respectively. The levels of both calculated quantiles were spanned evenly from 0 to 1. Q_x describes the overall distribution of apparent maxima in time, reporting changes in the rate of resonance occurrence. We deliberately used normalisation by T_max because we want the results to be sensitive only to the relative distribution of the resonance events over the lifetime of the EEV. Otherwise, its values would be strongly correlated with the length of the resonance curve itself¹⁶, rather than with the distribution of resonances over time. On the other hand, Q_y encapsulated the combined information about the mean level of log(ℒ), any long-term trends in the resonance curve and the distribution of the heights of the maxima. In contrast to the Q_x, we did not apply any normalisation to Q_y as its absolute values carry valuable information about the strength of the resonances and the average level of the entire resonance curve. The final vector Q was constructed as the concatenation of Q_x and Q_y, which had previously been scaled using the variance in the sets of all Q_x and Q_y. The resulting Q has a total of 2000 dimensions.

In the next step, we performed a preliminary dimensionality reduction of Q by means of principal component analysis (PCA; Pearson 1901), obtaining pre-processed ‘morphology’ vectors, θ_PCA. Principal component analysis is a method that orthogonally projects the data into a coordinate system in which successive vector components explain a smaller and smaller part of the data variance. The target number of its dimensions returned by PCA for each Q was set to 10. This value was chosen experimentally by examining the percentage of the total variance of the dataset explained by successive PCA components. For ℒ₁(t), the first ten PCA components explained a total of 99.8% of variance (first component – 79% and second component – 19%). For ℒ₂(t), the corresponding value was 98.5% (in this case, the first PCA component explained 60% of the total variance, while the second explained 22%).

We then performed the final 2D embedding by applying UMAP on the collection of θ_PCA vectors. Uniform Manifold Approximation and Projection is a multipurpose non-linear dimensionality reduction technique that constructs a low-dimensional projection that preserves as accurately as possible the topological structure of the input data. For instance, in this case, a pair of embeddings of resonance curves with similar properties (in the sense of their summary statistics described above) are expected to lie in mutual vicinity on the 2D UMAP plane. We denote the UMAP results as θ_UMAP. The manifold spanned by θ_UMAP (Sect. 4.5) allowed us to effectively examine the differences in the morphology of the resonance curves and their dependence on the initial parameters of the simulated EEVs.

Unlike PCA, UMAP is a complex method, with many free parameters that need to be adjusted with care, as the resulting embedding may depend heavily on their choice. Appendix D provides all the ‘technical’ details of this process, including the values of the most important UMAP parameters adopted in this study.

4. Results

4.1. General properties of synthesised EEVs

Before going into a detailed analysis of the resonance curves, we briefly characterise the general properties of the models we synthesised using the MESAbinary and GYRE codes.

4.1.1. Evolutionary tracks in the HRD

Figure 4 shows a pair of HRDs with a compilation of all 20 000 evolutionary tracks that we obtained in our simulations for the primary (Fig. 4a) and secondary (Fig. 4b) components. Although it is impossible to clearly present thousands of evolutionary tracks on a single HRD, we have highlighted and colour-coded a small fraction of them in order to describe some of their features.

First of all, only a fraction of the primaries reached the TAMS when the central mass abundance of hydrogen dropped below 10⁻⁴ (according to the first of our termination conditions, Sect. 3.2.3). Many evolutionary tracks were interrupted at MS due to the fulfilment of one of the other termination conditions. Secondly, a number of tracks clearly change their character after crossing the line corresponding to the bi-stability jump (around T_eff = 26 000 K). This is due to the associated sharp increase in the wind mass-loss rate, as it tries to keep the stellar luminosity constant. In some circumstances, the mass-loss rate is so high that the star loses a significant part of its envelope¹⁷. This effect ‘pushes’ the star back to the high effective temperature region and is particularly pronounced for the most massive stars in our sample (cf. Fig. 4a, evolutionary tracks that ‘turn around’ and cross the bi-stability jump for a second time).

4.1.2. EEV groups in terms of the termination condition

Only four of the seven¹⁸ termination conditions actually occurred in our simulations. The majority of our EEVs (∼67.1%) ended up as MS RLOF systems in which the primary component filled its Roche lobe during the periastron passage. The next most numerous group (∼22.6%) were systems in which the primary component successfully reached the TAMS (X_c ≤ 10⁻⁴). About 10.2% of the binaries managed to circularise their orbits before any other termination condition was met. The last group contains only about 0.04% of the total sample. This is the group where the primary’s rotation velocity exceeded the maximum allowed angular velocity (Ω/Ω_crit ≥ 0.75).

Figure 5 presents these four groups of EEVs on the P_orb − e plane and allows a comparison of the initial (Fig. 5a) and final (Fig. 5b) states of the aforementioned distribution. As expected, the EEVs with the shortest orbital periods and high eccentricities tend to circularise their orbits before leaving the MS. Their trajectories in the P_orb − e diagram (Fig. 5c) follow smooth, almost vertical lines due to the strong tidal damping of eccentricity. On the other hand, the integration of the evolution of systems with large distances between components at periastron (${\stackrel{\sim }{r}}_{\mathrm{peri}}\gtrsim 3.5$) has been terminated mainly due to the exhaustion of hydrogen in the primary’s core. Although the majority of EEVs belonging to this group do not significantly change their orbital parameters during evolution, there is a subgroup of them that behaves quite differently. It can be recognised as the distinct ‘cloud’ of green dots in Fig. 5b, represented by the mainly horizontal green tracks in Fig. 5c. These are systems that were characterised by very strong stellar winds at the end of the MS phase and have lost much of their envelopes, so that their orbital period has increased significantly (Kepler’s third law).

The most numerous group of EEVs, in which the primary component has filled its Roche lobe in the MS phase, forms a kind of ‘bridge’ between the two previously mentioned groups and merges with them. The shapes of the corresponding trajectories on the P_orb − e plane may vary from system to system, depending on the interplay between tidal forces and the intensity of stellar winds, so no single ‘type’ of track can be assigned to them. However, many of them resemble the inverted Greek letter ‘Γ’ – initially, the system drifts horizontally (towards the longer orbital period) as a result of the mass loss and/or spin-orbit coupling, and then undergoes more or less rapid circularisation (moves vertically downwards) under the influence of the intense tides, which come to the fore when the primary component almost fills its Roche lobe.

Only 8 out of 20 000 EEVs underwent efficient spin up of both components due to the pseudo-synchronisation (when the rotation period of the star ‘matches’ the rate of orbital motion at periastron, so that there is no net torque over an orbital cycle; e.g. Hut 1981). These few systems are located in the upper right corner of Figs. 5a,b. Their orbits were initially highly eccentric yet relatively widely separated at periastron (${\stackrel{\sim }{r}}_{\mathrm{peri}}\approx 4.5-5.0$). Thus, in combination with the lower masses of the primary components (M₁ ≈ 5 M_⊙), there was no effective tidal dissipation. However, the envelopes of these stars tended to rotate pseudo-synchronously with the orbit (due to the relatively long nuclear timescale of the evolution of a 5 M_⊙, the primaries had enough time to do so). Consequently, this led to a very fast rotation of the primary component, exceeding the threshold of Ω/Ω_crit = 0.75.

4.1.3. Internal structure and asteroseismic properties

The shape of the resonance curve depends not only on the global properties of the components and the orbit, but also on the internal structure of the stars, which directly affects seismic properties (i.e. the spectrum of eigenmodes). Therefore, within the limited volume of this paper, we show at least one representative example of the evolution of the internal properties of the primary component for an arbitrarily chosen EEV. Figure 6 shows the evolution of a primary with mass M₁ ≈ 13.6 M_⊙ in a system with an initial eccentricity e ≈ 0.4 and an initial orbital period P_orb ≈ 4.0 d. In our simulations, this particular system finished its evolution due to the circularisation of its orbit after about 12 Myr. The HRD in Fig. 6 reveals the ‘non-standard’ evolutionary track of the primary due to the sharp change in the mass-loss rate after crossing the bi-stability jump (right panel in the top row of Fig. 6). The same panel also shows how the primary’s surface rotation rate varies over time – as the mass-loss rate increases, it loses a lot of spin AM and slows down its rotation. The eccentricity and orbital period monotonically decrease with time (middle panel in the top row of Fig. 6), except for a short episode of increase in P_orb caused by the irreversible loss of a large part of the envelope. We selected three epochs in the evolutionary history of this EEV (labelled A, B, and C on the HRD), for which we present the appearance of the rotational profiles, mode propagation diagrams and oscillation spectra of the primary component in the bottom part of Fig. 6. Epoch A corresponds to the phase of evolution just after leaving the ZAMS, epoch B is characterised by X_c, 1 ≈ 0.45, and finally, epoch C marks the situation just before the complete circularisation of the EEV. Below we briefly describe the changes occurring in each of the three types of diagram below.

The internal rotation profile of the primary is almost constant for epoch A, but by then a division between a faster-rotating core and a more slowly rotating envelope begins to emerge. The aforementioned division becomes particularly apparent in epoch B, when the core has developed a rotation rate approximately 1.25 times that of the surface layers. As can be seen, the contracting core rotates as a rigid body throughout the MS lifetime due to efficient AM transport supported by convection. The outer part of the envelope also rotates almost rigidly, but this time it is due to large-scale Eddington–Sweet meridional flows. The angular velocity gradient in the primary starts to gradually decrease as the star reaches epoch C. Various mixing processes in the chemically modified layer left by the core lead to the diffusion of AM from the core to the envelope. Moreover, the rotational profile inside the star becomes a smooth function of the radius (rather than a step-like function as for epoch B).

The majority of TEOs in our simulations belong to the g-mode family of oscillations, so it is very important to control the behaviour of the Brunt–Väisälä buoyancy frequency, N_BV, in our models. Together with the Lamb frequency for l = 2 modes, S_l = 2, they carry information about g and p mode cavities and their evanescence regions (e.g. Aerts et al. 2010, their Sect. 3.4). The evolution of N_BV and S_l = 2 is presented in the middle column of Fig. 6, which shows the position of the l = 2 p-mode and g-mode propagation cavities. The white areas that lie between the Brunt–Väisälä and Lamb frequencies correspond to the evanescence regions. During evolution, the receding convective core builds up a large g-mode trapping cavity, which is very important for their frequency spectrum. Additionally, the behaviour of the N_BV just below the stellar photosphere reveals a pair of thin subsurface convection zones, expected for this type of star (e.g. Jermyn et al. 2022). Comparing the mode propagation diagrams for epochs A and C, it can be seen that also the p modes can penetrate deeper and deeper into the primary as it gradually depletes the hydrogen in its core.

The right column in Fig. 6 contains most of the information that is directly used to obtain the resonance curve, showing the frequency range in which GYRE looked for potential TEOs (according to the criteria adopted in Sect. 3.3). We recall that that their width depends on the $N_{\mathrm{m}\mathrm{a}\mathrm{x}}^m(e)$ functions, so as the system evolves towards lower eccentricities, this range is shorter and shorter (i.e. fewer harmonics of the orbital frequency can effectively drive TEOs). We also mark the location of the tidal forcing frequencies. As can be seen, the separation between successive values of f_Nm becomes larger with passing time due to the increase in f_orb. Also shown are the eigenfrequencies found by GYRE and the linear damping rates of these modes. The presented set of three synthetic oscillation spectra reveals a typical structure for g modes with their asymptotic behaviour for large radial orders (which correspond to lower frequencies). It may appear that the dense ‘forests’ of eigenfrequencies end too early relative to the left limits of horizontal bars. However, this is not a mistake, but a direct consequence of the maximum |n| we allowed in the calculations – modes with lower frequencies would have radial orders larger than 30. During the evolution of the EEV, both the forcing frequencies and the oscillation spectrum shift, so the intersection of these two vertical line patterns is virtually inevitable in most cases. Each of these intersections is the source of a single resonance that can give rise to a noticeable TEO at the level of the photosphere.

4.2. The ‘visual inspection’ of resonance curves

The resonance curves are characterised by a striking diversity in terms of morphology, which is already partly evident in Fig. 3. The four examples of ℒ₁(t) and ℒ₂(t) shown in this figure show that the components of the EEVs can experience, firstly, a very different number of resonances and, secondly, their distribution in time can take various forms. The heights of the maxima of the resonance curves are mainly dictated by the γ_nlm of the mode to which the smallest difference corresponds, (σ_nlm − f_Nm). Statistically speaking, modes with larger |n| are more strongly non-adiabatic (have larger damping rates), and hence the maxima they induce in the resonance curves are lower (cf. Eq. (3)). Another factor determines the extent of the resonant maximum in time. It is determined by the relative ‘velocity’ with which the eigenfrequency spectrum crosses the f_Nm spectrum. By ‘velocity’ here we mean the rate of change of these two independent frequency spectra.

It should be emphasised that there are also numerous cases in which ℒ(t) drops sharply to zero at some point (cf. ℒ₂(t) curve in the bottom panel of Fig. 3) or resonances do not occur at all (see Sect. 4.3 and Fig. 8). Such a situation can occur, for example, when the oscillation spectrum lies completely outside the frequency range allowed by the F_Nm coefficients or the nuclear timescale of the secondary is much longer than the same timescale for the primary. Under such circumstances, the secondary component remains close to the ZAMS until the termination condition is met. Thus, it will not significantly change its internal structure and oscillation spectrum. This in turn means that the oscillation spectrum will not move relative to the tidal forcing frequencies, effectively reducing the number of possible resonance events.

Our sample of resonance curves includes a particular group of ℒ(t) curves that exhibit exceptionally long duration resonances compared to typical ones (we refer to them as ‘long resonances’). Figure 7 presents parts of three representative resonance curves belonging to this group. The shaded regions in the figure mark the position of the long resonances. As can be clearly seen, the typical resonance usually lasts for about 10³ − 10⁴ yr, which is approximately 100 times shorter than the duration of a long resonance (of the order of 10⁵ − 10⁶ yr). They originate from the intersection of one of the f_Nm frequencies with the σ_nlm frequency at a very small angle, in terms of their temporal evolution. As a result, they remain for a relatively long time in very close vicinity, leading to a broad resonance overlapping with narrower ones (originating from other intersections of the f_Nm and σ_nlm frequency spectra; cf. especially the middle panel in Fig. 7). The long resonances are interesting for at least two reasons. First of all, they are natural candidates for resonantly locked TEOs. However, based on our simulations, it is difficult to say whether an extended resonance would persist when the back-reaction of a TEO on the orbit is taken into account. Secondly, if the energy exchange between the eigenmode and the orbit that corresponds to a long resonance is not efficient (i.e. there is a small chance that a long resonance will be lost), such a resonance should lead to a high-amplitude TEO without the need for resonance locking. This is simply because it has enough time to reach its saturation level due to the non-linear effects. However, we do not find any significant correlations between the occurrence of a long resonance in ℒ(t) and the initial parameters of our EEVs.

4.3. Total number of resonances and the average rate of their occurrence

The first feature of the morphology of the resonance curves that we investigated is the total number of resonances that occurred in the primary and secondary components, 𝒩_res, 1 and 𝒩_res, 2. However, we did not calculate these statistics directly from ℒ₁(t) and ℒ₂(t), because some of the apparent maxima may actually be a blend of more than one resonance event. This is especially true when the involved γ_nlm differ by orders of magnitude. Then one of the resonances is characterised by a notably smaller maximum, which seems to ‘hide’ in the dominant one. Instead, we used a different approach that do not underestimate the actual number of resonances. When post-processing the generated models, we simply counted each intersection of the σ_n, 2, 0 and σ_n, 2, +2 frequency spectra with their counterparts f_N, 0 and f_N, 2, respectively. The results of such an analysis are depicted in Fig. 8.

The most important thing about Fig. 8 is that it shows the absolute number of resonances. Eccentric ellipsoidal variables can experience hundreds or even thousands of resonances during their evolution on MS. It would therefore be wrong to claim that these phenomena are rare in massive and intermediate-mass EEVs, although in general, resonances are quite short-lived compared to the nuclear timescale. The total number of resonances experienced by the primary component (Fig. 8a) shows a correlation with both its initial mass and the initial eccentricity of the system. The mildly decreasing trend of 𝒩_res, 1 towards higher M₁ originates from the fact that the mean lifetime of the star on MS shortens with increasing mass. On the other hand, the wide range of 𝒩_res is mainly due to differences in initial eccentricity. The closer the system is to a circular geometry at the beginning of evolution, the statistically lower the value of 𝒩_res, which is self-explanatory and also applies to the secondary component (Fig. 8b). Eccentric ellipsoidal variables that have managed to circularise their orbits in the MS phase (magenta dots in Fig. 8c) have on average lower initial eccentricities and thus fewer resonances. The opposite behaviour is exhibited by EEVs in which the primary component has had a chance to reach the TAMS (green dots in Fig. 8c). The secondary components experience a slightly fewer resonances compared to the primaries (Fig. 8b) and there is no clear division of the 𝒩_res, 2 distribution with respect to the termination criterion (Fig. 8d). The noticeably smaller number of resonances for secondary components with masses M₂ <  5 M_⊙ comes from the conditions of our simulations: secondaries with these masses occur in systems with decreasing mass ratios¹⁹. Hence, the large difference in nuclear timescales between the components means that the secondary component does not significantly change its eigenfrequency spectrum, resulting in a smaller number of resonances.

As we already mentioned in Sect. 4.2, some of the ℒ₁(t) and ℒ₂(t) curves do not reveal any resonances, which is why they lie in Fig. 8 on the horizontal line 𝒩_res = 0. This behaviour occurs for only 0.07% of our primaries. They are all EEVs with highly eccentric orbits that quickly filled their Roche lobes at periastron. There is much more such behaviour for secondaries, about 7%, mainly for the intermediate-mass companions of the much more massive primaries.

From an observational point of view, even more important than the total number of resonances is the rate at which they occur. Knowing this rate, for a given population of MS EEVs, we can approximately say where we have a statistically higher chance of observing TEOs. After all, our observations only correspond to one particular moment in time, not the entire evolution. Knowing the values of 𝒩_res for each component and the age of each system at termination, T_max, we can calculate the average rate of resonances as ℛ_res ≡ 𝒩_res/T_max. We show the distribution of ℛ_res, 1 and ℛ_res, 2 in Fig. 9. It is very difficult to predict what the dependence of ℛ_res on the mass of the component will look like, as it is the result of a complex interplay between many related factors. On the one hand, it can be said that massive stars should have a smaller ℛ_res because their lifetimes are shorter and they fill their Roche lobes relatively easily (in the considered range of orbital parameters). On the other hand, however, massive stars quickly change their internal structure (i.e. asteroseismic properties), so that their eigenfrequency spectra evolve rapidly, increasing the likelihood of interaction with the structure of the tidal forcing frequencies. The question is, which of these processes prevails? As can be seen in Fig. 9, it is the more massive stars that are more likely to undergo resonances. Both primary and secondary components with masses around 30 M_⊙ have on average an order of magnitude higher ℛ_res (∼10² Myr⁻¹) than components with masses around 5 M_⊙ (∼10¹ Myr⁻¹, Figs. 9a,b). Moreover, the dependence of the distributions shown in Fig. 9 on the initial eccentricity and termination condition is inherited from Fig. 8.

At this point, we can venture the conclusion that in the case of MS EEVs, TEOs should be observed mostly in the upper part of the MS (among early B- and O-type dwarfs), which still requires observational verification on a large sample of massive EEVs. Although we cannot extrapolate the obtained distributions of ℛ_res towards lower masses, these stars have an increasingly extended convective envelope, which in turn should effectively limit the photometric detection of g-mode TEOs. On the contrary, the envelopes of massive stars are radiative, which should not prevent g-mode TEOs from propagating up to the vicinity of the photosphere. Thus, they can be more easily detected by analysing the light curves, especially in the era of high-quality space-borne photometry.

4.4. Distribution of resonances over time

Since the average rate of resonances we have studied so far has effectively obliterated any differences in the corresponding temporal distribution, we can ask another important question: Are there any distinctive moments in the evolution of the simulated EEVs during which the systems experienced temporally higher resonance rates? After visually inspecting hundreds of resonance curves, we notice that the aforementioned rate changes dramatically in many cases (cf. the top panel of Fig. 3 as an example). In order to compare the temporal distribution of resonance events for the various EEVs we are dealing with, we performed this type of analysis on subgroups of systems divided according to the termination condition. We also normalised the time variable by dividing it by T_max of each resonance curve. This allowed us to present the whole evolution of components on a convenient and uniform interval, [0, 1]. Figure 10 shows the results obtained for the primary components that have managed to deplete hydrogen in their cores.

Figure 10 also demonstrates that the distribution discussed here is not uniform over time. Specific areas in this diagram are clearly distinguishable. Nevertheless, this figure still contains information on the total number of resonances, which makes it somewhat problematic to compare the shapes of these distributions for different masses of the components. We, therefore, prepared histograms of the times of resonances for five intervals of the primary’s initial mass (every 5 M_⊙). Separate sets of histograms were generated for the primary and secondary components and the three main termination conditions²⁰. All histograms are shown in Fig. 11.

The most diverse structure of the temporal distribution of resonances is shown by systems in which the primary component has completed its evolution in our simulations at the TAMS (Figs. 11a,b). In fact, for all mass ranges, the distribution has two distinct maxima. The smaller of the two is located near the ZAMS, while the other is just before reaching the TAMS. Their presence can be explained by the rate of change in the stellar eigenspectrum, which is the highest (after averaging over all modes) at the aforementioned moments of evolution. In particular, the rapid changes in the radius of the star when it is close to complete depletion of hydrogen in its core cause a very high ‘concentration’ of resonances in the final MS phase. The height of this dominant maximum decreases towards higher masses, but at the same time it becomes wider and wider. The distribution for the secondary components also reveals this kind of maximum near the TAMS, which is particularly well pronounced for companions of primaries with masses ≳25 M_⊙. Given these facts, an interesting conclusion can be drawn. Massive and intermediate-mass EEVs with at least one component leaving the MS should experience an increased rate of encountered resonances. One might therefore suspect that there is a statistically higher chance of observing TEOs in more evolved EEVs. The properties of the histograms for EEVs in which RLOF eventually occurred at the periastron (Figs. 11c,d) are very similar to the case described above. The only evident difference between the two is the reduction in maximum of the distribution near the TAMS. This is due to the fact that the primary is likely to start the RLOF earlier than it reaches the TAMS, preventing the occurrence of a large number of resonances in a relatively short time, as mentioned earlier.

Eccentric systems that are subject to effective circularisation (Figs. 11e,f) behave quite differently from the two previous cases. They experience the vast majority of their resonance phenomena at the beginning of evolution, and then reduce the number of resonances almost monotonically, as the orbital eccentricity becomes smaller and smaller with time. Hence, the chance of observing TEOs in initially relatively tight EEVs (cf. Fig. 5a) is largest in the vicinity of the ZAMS, which stays in contrast to the systems described above.

4.5. Investigation of the morphology of resonance curves using UMAP

All the analysis described above is based solely on the distribution of resonances in time, meaning we neglected the actual morphology of the resonance curves, for example differences in the height and width of resonance maxima, mean level of ℒ(t), and long-term trends in ℒ(t). Using the dimensionality reduction techniques presented in Sect. 3.5, we constructed 2D UMAP embeddings of the space of resonance curves in terms of their morphological features. Figures 12 and 14 show the results obtained for the resonance curves of the primary and secondary component, respectively. We recall that the idea of the low-dimensional embedding performed here is to preserve the distances between two points in the original space as accurately as possible, so that the distances in the 2D plane reflect the distances in the full (original) space of morphological features (the vector of 2000 quantiles, Q). In other words, a pair of distant points in Figs. 12 and 14 should correspond to resonance curves with notably different morphologies and, vice versa, a pair of resonance curves with similar properties is expected to lie in mutual vicinity on the 2D UMAP plane. Thanks to this key property of UMAP and many other dimensionality reduction methods we can effectively explore the entire space of resonance curve morphologies.

4.5.1. UMAP plane for primary components

We begin with a discussion of Fig. 12. Firstly, the presented 2D embedding does not indicate the presence of any well-separated groups among the resonance curves for the primary components. This is an observation that is true over the entire range of different values of the UMAP free parameters (Appendix D) as well as for the different summary statistics of the resonance curves that were considered during the preliminary experiments. The morphology of the resonance curves changes smoothly depending on to the initial parameters of the simulated EEVs.

Secondly, as can be seen in Figs. 12b,c, the initial eccentricity and normalised periastron distance are parameters strongly correlated with the overall morphology of the resonance curves of the primary components. Moreover, their gradients in the UMAP plane are approximately orthogonal. Therefore, the pair of these parameters is the primary factor that determines the shape of ℒ₁(t). The termination condition (Fig. 12e) generally follows the behaviour of ${\stackrel{\sim }{r}}_{\mathrm{peri}}$ except at small periastron distances, when the morphology remains similar but the simulations were terminated due to hydrogen depletion in the primary’s core or near-complete circularisation of the orbit. The initial mass of the primary component (Fig. 12a) and its initial angular velocity of rotation (Fig. 12d) are second-order factors shaping the morphology of the ℒ₁(t) resonance curves. In the inner part of the plane, M₁ is distributed almost randomly. The clear exception is the boundary of the plane, which can be roughly divided into two parts of mostly high or low initial mass of the primary. A similar conclusion can be drawn for the initial Ω₁/Ω_crit, 1. This time, however, the upper right part of Fig. 12d reveals a well-defined group of high initial Ω₁/Ω_crit, 1 and ${\stackrel{\sim }{r}}_{\mathrm{peri}}$.

It is difficult to include here a complete presentation of the changes in the morphology of ℒ₁(t) as a function of their position on the UMAP plane. Therefore, we only focus on some extreme points to present some boundary cases. Figure 13 shows examples of ℒ₁(t) from different areas of the morphological plane. Primary components with lower masses and high initial eccentricities are generally characterised by resonance curves with high mean levels and a rich set of resonances, as shown in Fig. 13a. The resonance curve depicted in Fig. 13b represents intermediate-mass fast-rotating primary with large initial ${\stackrel{\sim }{r}}_{\mathrm{peri}}$ and low initial eccentricity. Here, the base level of ℒ₁(t) increases by an order of magnitude and then the system experiences a large number of resonances, during the evolution near the TAMS. The increase in the mean value of ℒ₁(t) is characteristic of stars with a high initial rotation rates. Primary components lying between points a and b generally do not manifest this characteristic. Moving along a straight line on the plane from a to b, the increase in the number of resonances during the evolution near the TAMS becomes more and more apparent. Figure 13c shows an example of a system with a small initial eccentricity and a short periastron distance at the ZAMS that is rapidly circularising. As expected, the ℒ₁(t) resonance curves for such objects have a small number of resonance maxima and a low base level. The case corresponding to the larger initial eccentricity is shown in Fig. 13e, where the total number of resonances is much greater. In this case, the evolution is mainly distinguished by a decrease in the frequency of resonances with time, related to the efficient circularisation of the orbit, and therefore a decrease in the mean level of ℒ₁(t). Finally, the resonance curve in Fig. 13d is representative of the most EEVs between points c and d. They are characterised by an approximately uniform distribution of resonances over time and an almost constant base level of the resonance curve.

4.5.2. UMAP plane for the secondary components

The situation for the secondary component (Fig. 14) is quite different from the previous case. The UMAP manifold obtained for the set of ℒ₂(t) reveals slightly more complex structure than the shape of embedding in Fig. 12. Since the time span of ℒ₂(t) is largely determined by the mass of the primary component, the resonance curves for the secondary components were terminated at times not necessarily related to their actual evolutionary status and are statistically shorter than they could be for primaries of the same mass. For this reason, secondary components experience, on average, fewer resonances, but, at the same time, their resonance curves can take more diverse forms compared to ℒ₁(t). Undoubtedly, the main factor shaping the morphology of the ℒ₂(t) is the initial eccentricity (Fig. 14c), which with the exception for two small areas, varies smoothly across the UMAP plane. The other parameters (Figs. 14a–d) play a secondary role, showing the complex and fine structures on the plane. As can easily be seen in Fig. 14a, the extreme cases of ℒ₂(t) in terms of their morphology belong almost exclusively to the EEVs with intermediate-mass secondary components that gather at the periphery of the plane. Some of these objects even form slightly better separated groups, isolating from the central part of the area.

Five limiting examples of resonance curves for secondary components are shown in Fig. 15. Panels a, c and e together with the area they approximately enclose contain resonance curves morphology very similar to that described for primary components. The resonance curves belonging to the ‘clouds’ of points labelled as b and d in Fig. 15 are completely different. They are distinguished by the complete absence of resonances during a certain period of the evolution of the system. The differentiating feature of these cases is the disappearance of resonances from some time to the end of evolution (Fig. 15b) or the presence of resonances only around the middle of the considered evolution time (Fig. 15d).

5. Summary and conclusions

In our paper we aimed to investigate the temporal variation of conditions that favour the excitation of TEOs in EEVs with massive and intermediate-mass MS components (Sect. 2) and see how their picture changes with different initial system parameters. In order to achieve this goal, we simulated the evolution of 20 000 EEVs using the MESA software in combination with the GYRE stellar oscillations code (Sect. 3). Our calculations started at the ZAMS and were terminated if one of the conditions presented in Sect. 3.2.3 was met. We considered only modes with l = 2, m = 0, +2 because they are expected to be dominant TEOs. We also assumed that all TEOs are due to chance resonances, that is to say, we neglected the effect of TEOs on the orbit. Knowing the evolution of the orbital parameters of simulated EEVs and the temporal changes in the eigenmode spectra of the components, we were able to derive resonance curves ℒ₁(t) and ℒ₂(t), defined by Eqs. (3) and (4). The equations reflect the overall resonance conditions and, thus, indirectly also the chance of TEOs occurring, separately for the primary and secondary components of our simulated EEVs.

After visually inspecting the obtained resonance curves, calculating basic statistics for them, and applying ML-based methods to the entire dataset, our main results can be summarised as follows.

1. Resonance curves are characterised by striking diversity in terms of their morphology (Sect. 4.2). The components of EEVs can experience a very different number of resonances, and their distribution over time can take various forms, including a lack of resonances over long periods of time. We also distinguished a group of resonance curves that exhibit prolonged resonances, about two orders of magnitude longer than typical (Sect. 4.2, Fig. 7). These long resonances are the potential sources of high-amplitude and resonantly locked TEOs.

2. Resonances between tidal forcing frequencies and the spectrum of stellar normal modes are not rare events among massive and intermediate-mass MS EEVs (Sect. 4.3). Although the total number of resonances depends mostly on the initial orbital parameters, it is typically of the order of 10² − 10³ for a given system during the entire MS phase (Fig. 8). We emphasise at this point that these numbers are lower limits for the actual 𝒩_res in EEVs because we considered only l = 2 TEOs. Taking higher-degree modes into account will certainly increase the reported values of 𝒩_res.

3. On average, the more massive a star is, the higher the rate of resonances it experiences (Sect. 4.3). For the most massive stars in our sample (≈30 M_⊙), the average rate of resonances can reach ∼10² Myr⁻¹, which is approximately an order of magnitude higher than for intermediate-mass stars (Fig. 9).

4. The distribution of resonances over time is not homogeneous and depends primarily on whether the system circularises before the primary reaches the TAMS or RLOF occurs at the periastron (Sect. 4.4, Fig. 11). We noticed a particular moment in the evolution of our EEVs near the TAMS, when the components undergo an increased number of resonances in a relatively short time (Figs. 11a,b).

5. The low-dimensional representation of the morphology of the resonance curves, summarised by quantile-based statistics and subsequently processed by UMAP, shows that its manifold forms a rather smooth distribution without well-defined (separated) groups (Sect. 4.5, Figs. 12 and 14). Less differentiated, at least in terms of the adopted method, are the resonance curves of the primary components, for which the initial eccentricity and the normalised periastron distance largely determine their morphological features. Although secondary components experience far fewer resonances, their shapes are generally more complex due to the predominant influence of the primary component on evolution time.

In light of the results obtained in our study, we can draw several interesting conclusions. Firstly, statistically speaking, TEOs are more likely to be discovered in more massive EEVs, as their components have a higher average rate of resonances. This does not necessarily mean that a higher absolute number of more massive EEVs exhibiting TEOs than less massive ones will be observed²¹. However, when comparing two particular systems, one with intermediate-mass components and the other with much higher-mass components, it is for the latter that we have a statistically higher chance that some resonance is currently underway there. Secondly, it seems that TEOs should be especially visible in EEVs that contain a component approaching the TAMS. Given these facts, the ‘extreme-amplitude’ massive EEV MACHO 80.7443.1718 (Jayasinghe et al. 2021; Kołaczek-Szymański et al. 2022) fits this picture almost perfectly. Its primary component is a B0.5 Ib-II supergiant leaving the MS, and, more importantly, many high-amplitude TEOs have now been detected in this extreme system. It is possible that what the primary component of MACHO 80.7443.1718 is currently undergoing corresponds to the resonance curve shown in Fig. 13b, that is to say, it is in the phase of a high resonance rate caused by relatively fast changes in its radius and orbital parameters. Moreover, the amplitudes of TEOs observed in MACHO 80.7443.1718 vary notably over time, suggesting that we may witness rapid changes in the resonance conditions for the primary component of this particular EEV. It would therefore be very valuable to carry out an observational study for a large sample of EEVs to verify whether TEOs are common in massive and intermediate-mass EEVs whose components have already depleted most of the hydrogen in their cores. With high-quality space-borne photometric observations, both operational, such as the Transiting Exoplanet Survey Satellite (Ricker et al. 2015) or BRITE-Constellation (Weiss et al. 2014), and planned missions, such as Planetary Transits and Oscillations of Stars (Rauer et al. 2014), it is definitely a feasible task.

The excitation of g-mode TEOs, which propagate deep inside the star, may be an underestimated mechanism for AM transport inside the components of EEVs. It has long been suspected that self-excited oscillations and internal gravity waves²² efficiently redistribute AM in the radial direction of the star (e.g. Rogers et al. 2013; Rogers & McElwaine 2017). Although the majority of our resonances have relatively short durations (of the order of 10³ − 10⁴ yr), they can be quite frequent (especially near the TAMS), and hence the question of their contribution to AM transport and mixing processes becomes urgent for the components of EEVs. Performing calculations that would treat the evolution of the orbit, components, and TEOs in a fully self-consistent way seems particularly interesting for massive eccentric systems leaving the MS.

We have already entered the era of observational studies of the distant star-bursting galaxies and stellar populations in low-metallicity environments that shaped the Universe in its early epochs. Recalling that the metal-poor stars were much more massive than their current metal-rich counterparts (e.g. Hosokawa et al. 2013; Susa et al. 2014), we can ask what effect metallicity has on the occurrence of TEOs in massive EEVs and how the results we have presented depend on metal contents. Therefore, future studies of the importance of TEOs in massive, metal-poor EEVs seems worthy of further investigation, especially in light of the ongoing James Webb Space Telescope mission²³ (JWST; Gardner et al. 2006), which is certain to bring many discoveries regarding the stellar astrophysics of early stellar populations, including stars in eccentric binary systems.

---

Acknowledgments

P.K.S. is indebted to his brother, Adam Karol Kołaczek-Szymański, who oversaw the purchase and assembly of a dedicated PC workstation to enable the efficient calculation of models in MESA and GYRE. Without his generous help, this project would literally never have been completed. The authors are thankful to Prof. Andrzej Pigulski for many important suggestions and fruitful discussions that made this manuscript more comprehensible, and to the anonymous referee for many inspiring comments that helped to improve the manuscript. P.K.S. was supported by the Polish National Science Center grant no. 2019/35/N/ST9/03805. T.R. was partly founded from budgetary funds for science in 2018–2022 in a research project under the program “Diamentowy Grant”, no. DI2018 024648. Much of this work was developed and written during IAU Symposium 361, “Massive Stars Near & Far”, held in Ballyconnell, Ireland, 8–13 May, 2022. P.K.S. is very grateful to the organisers for the opportunity to participate in this event.

References

Appendix A: MESA input files

MESAbinary needs at least three input files (hereafter ‘inlists’) to start the calculations. Two of them provide all the necessary parameters to perform the evolution of each component separately²⁴, and the last one describes the evolution of the orbit as well as other processes that depend on binarity²⁵ (e.g. spin-orbit coupling). In the following appendix, we present example MESAstar and MESAbinary inlists, which we used to generate a set of the evolutionary tracks of the components of the binary system. However, we have intentionally omitted any controls related to the names of the files or directories where the results should be stored. All values that changed in the inlists depending on the simulated binary system were enclosed in square brackets – [ ⋅ ]. The parameters adopted below resulted in a typical number of about 2,400 zones in the radial direction of the star. The number of models calculated per single EEV was typically of the order of several hundred, mainly depending on the termination condition.

A.1. MESAstar inlist

&star_job

!OUTPUT

history_columns_file = "my_history_columns.list"

profile_columns_file = "my_profile_columns.list"

show_log_description_at_start = .false.

save_photo_when_terminate=.false.

!MODIFICATIONS TO MODEL

new_rotation_flag=.true.

change_rotation_flag=.true.

new_omega_div_omega_crit=[Ω/Ω_crit]

num_steps_to_relax_rotation=100

relax_omega_max_yrs_dt = 1d4

relax_omega_div_omega_crit=.true.

set_initial_cumulative_energy_error = .true.

new_cumulative_energy_error = 0d0

/ ! end of star_job namelist

MYAMPeos

/ ! end of eos namelist

MYAMPkap

use_Type2_opacities = .true.

Zbase = 0.02

/ ! end of kap namelist

MYAMPcontrols

!SPECIFICATIONS FOR STARTING MODEL

initial_z=0.02d0

!CONTROLS FOR OUTPUT

terminal_interval=100

write_header_frequency=1

photo_interval=100000

history_interval=5

star_history_dbl_format = "(1pes40.6e3, 1x)"

profile_interval=10

max_num_profile_models=5000

write_pulse_data_with_profile=.true.

pulse_data_format="GYRE"

add_double_points_to_pulse_data=.true.

!WHEN TO STOP

max_model_number = 5000

xa_central_lower_limit_species(1)="h1"

xa_central_lower_limit(1)=1d-4

omega_div_omega_crit_limit=0.75

!MIXING PARAMETERS

mixing_length_alpha=1.82d0

use_Ledoux_criterion=.true.

num_cells_for_smooth_gradL_composition_term = 0

alpha_semiconvection=0.01d0

okay_to_reduce_gradT_excess=.true.

mlt_make_surface_no_mixing = .true.

overshoot_scheme(1)="exponential"

overshoot_zone_type(1) = "burn_H"

overshoot_zone_loc(1) = "core"

overshoot_bdy_loc(1) = "top"

overshoot_f(1) = [f_ov]

overshoot_f0(1) = 0.005

do_conv_premix=.true.

set_min_D_mix=.true.

min_D_mix=1d5

!ROTATION CONTROLS

am_D_mix_factor=0.0333333d0

D_DSI_factor = 1

D_SH_factor = 1

D_SSI_factor = 1

D_ES_factor = 1

D_GSF_factor = 1

!ATMOSPHERE BOUNDARY CONDITION

atm_option="table"

atm_table="photosphere"

!MASS GAIN OR LOSS

hot_wind_scheme="Vink"

hot_wind_full_on_T=1.2d4

cool_wind_full_on_T=0.9d3

Vink_scaling_factor=1d0

no_wind_if_no_rotation=.true.

mdot_omega_power=0.43d0

max_mdot_jump_for_rotation=5d0

rotational_mdot_kh_fac = 1.0d3

!MESH ADJUSTMENT

max_delta_x_for_merge = 0.01d0

max_dq=1d-3

min_dq=1d-16

min_dq_for_split=1d-16

!ASTEROSEISMOLOGY CONTROLS

num_cells_for_smooth_brunt_B = 0

!STRUCTURE EQUATIONS

use_dedt_form_of_energy_eqn = .true.

!TIMESTEP CONTROLS

min_timestep_factor=0.5d0

max_timestep_factor=2.0d0

dH_div_H_limit=0.5d0

delta_lgL_phot_limit = 0.05d0

/ ! end of controls namelist

A.2. MESAbinary inlist

&binary_job

!OUTPUT/INPUT FILES

show_binary_log_description_at_start = .false.

binary_history_columns_file = "my_binary_history_columns.list"

!STARTING MODEL

evolve_both_stars=.true.

change_ignore_rlof_flag = .true.

new_ignore_rlof_flag = .true.

/ ! end of binary_job namelist

MYAMPbinary_controls

!SPECIFICATIONS FOR STARTING MODEL

m1=[M₁]

m2=[M₂]

initial_eccentricity=[e]

initial_period_in_days=-1

initial_separation_in_Rsuns=[a]

!CONTROLS FOR OUTPUT

history_interval=5

photo_interval=100000

terminal_interval=100

write_header_frequency=1

!TIMESTEP CONTROLS

fa=0.02d0

fa_hard=0.03d0

fr=0.10d0

fj=0.001d0

fj_hard=0.005d0

fe=0.02d0

fr_dt_limit = 1.0d2

fdm = 1d-3

fdm_hard = 5d-3

dt_softening_factor = 0.3d0

varcontrol_ms=5d-4

varcontrol_post_ms=5d-4

dt_reduction_factor_for_j=5d-2

!MASS TRANSFER CONTROLS

do_enhance_wind_1=.true.

do_enhance_wind_2=.true.

tout_B_wind_1 =[B_wind]

tout_B_wind_2 =[B_wind]

!ORBITAL JDOT CONTROLS

do_jdot_gr=.true.

do_jdot_ls=.true.

do_jdot_ml=.true.

do_jdot_mb=.false.

!ROTATION AND SYNC CONTROLS

do_tidal_sync=.true.

sync_type_1="Hut_rad"

sync_type_2="Hut_rad"

!ECCENTRICITY CONTROLS

do_tidal_circ=.true.

circ_type_1="Hut_rad"

circ_type_2="Hut_rad"

anomaly_steps=300

/ ! end of binary_controls namelist

Appendix B: GYRE input file

The GYRE stellar oscillations code requires a single input file that collects all the user-specified parameters of the asteroseismic calculations being performed²⁶. Below is our example file gyre.in. As in Appendix A, we have omitted any keywords related to specific file names and have highlighted variables by enclosing them in square brackets.

&constants

/

MYAMPmodel

model_type = "EVOL"

file_format = "MESA"

/

MYAMPmode

l = 2

m = 0

n_pg_min = -30

n_pg_max = 30

tag = "m0"

/

MYAMPmode

l = 2

m = 2

n_pg_min = -30

n_pg_max = 30

tag = "m2"

/

MYAMPosc

inner_bound = "REGULAR"

outer_bound = "VACUUM"

adiabatic = .true.’

nonadiabatic = .true.’

/

MYAMProt

coriolis_method = "TAR"

Omega_rot_source = "MODEL"

/

MYAMPnum

ad_search = "BRACKET"

nad_search = "AD"

diff_scheme = "MAGNUS_GL2"

/

MYAMPscan

grid_type = "LINEAR"

freq_min = [$f_{\min }^{m=0}$]

freq_max = [$f_{\max }^{m=0}$]

n_freq = [$N_{\mathrm{freq}}^{m=0}$]

freq_units = "CYC_PER_DAY"

grid_frame = "INERTIAL"

freq_frame = "INERTIAL"

tag_list = "m0"

/

MYAMPscan

grid_type = "LINEAR"

freq_min = [$f_{\min }^{m=+2}$]

freq_max = [$f_{\max }^{m=+2}$]

n_freq = [$N_{\mathrm{freq}}^{m=+2}$]

freq_units = "CYC_PER_DAY"

grid_frame = "COROT_I"

freq_frame = "COROT_I"

tag_list = "m2"

/

MYAMPgrid

/

MYAMPad_output

/

MYAMPnad_output

summary_file_format = "TXT"

summary_item_list = "freq,l,m,n_p,n_g,n_pg"

freq_units = "CYC_PER_DAY"

freq_frame = "INERTIAL"

The frequency scan limits, $f_{\min }^{m=0}$, $f_{\max }^{m=0}$, $f_{\min }^{m=+2}$, and $f_{\max }^{m=+2}$, were calculated as described in Sect. 3.3. The total numbers of discrete frequency points, $N_{\mathrm{freq}}^{m=0}$, and $N_{\mathrm{freq}}^{m=+2}$, were obtained as follows,

$$\begin{array}{c}\hfill N_{\mathrm{freq}}^{m=0,+2}=\lceil (f_{\max }^{m=0,+2}-f_{\min }^{m=0,+2})/(0.005{\mathrm{d}}^{-1})\rceil ,\end{array}$$(B.1)

where ⌈  ⋅  ⌉ denotes the ceiling function.

Appendix C: Data used by MESA

Our work uses the MESA stellar evolution code, which incorporates a vast compilation of knowledge, mainly from micro- and macrophysics, collected by many authors. The MESAeos module is a mixture of OPAL (Rogers & Nayfonov 2002), SCVH (Saumon et al. 1995), FreeEOS (Irwin 2004), HELM (Timmes & Swesty 2000), PC (Potekhin & Chabrier 2010), and Skye (Jermyn et al. 2021) equations of state. Radiative opacities are taken primarily from OPAL (Iglesias & Rogers 1993, 1996), with low-temperature data from Ferguson et al. (2005) and the high-temperature, Compton-scattering dominated regime by Poutanen (2017). Electron conduction opacities are from Cassisi et al. (2007) and Blouin et al. (2020). Nuclear reaction rates are from JINA REACLIB (Cyburt et al. 2010) and NACRE (Angulo et al. 1999), and additional tabulated weak reaction rates are from Fuller et al. (1985), Oda et al. (1994), and Langanke & Martínez-Pinedo (2000). Screening is included via the prescription of Chugunov et al. (2007). Thermal neutrino loss rates are taken from Itoh et al. (1996). Roche lobe radii in binary systems are computed using the fit of Eggleton (1983).

Appendix D: Adjustable parameters of the UMAP

UMAP, as a highly flexible method, is prone to returning misleading results in the case of inappropriately set free parameters. On the one hand, they can lead to the appearance of spurious groups and, on the other hand, to the loss of finer topological structure. The vital UMAP parameters that need adjustment are n_neighbors, min_dist, n_components and metric. The n_neighbors parameter is the most important, as it controls the balance between the local and global structure present in the data that will be mapped to the embedding. We experimented with different values of this parameter ranging from 5 to 1,000 (the default is 15) and concluded that the resonance curves (summarised by the proposed statistics) always form a single group, almost independently of the choice of n_neighbors. Later, min_dist sets the minimum distance between two different points on the embedding. We tested its values from 0.0 to 0.5 and do not observe any significant effect on manifold. We took the last two parameters, namely n_components that specifies the number of dimensions of the embedding, and metric specifying the metric used for similarity calculation, as default values. Finally, we used the following set of free parameters: n_neighbors = 500, min_dist = 0.1, n_components = 2 and metric = ’euclidean’.

All Figures