© The Authors 2024

1. Introduction

There is a general agreement that most massive stars with initial masses greater than ∼8 M_⊙ are in binary and multiple systems (Duchêne & Kraus 2013; Sana et al. 2014; Moe & Di Stefano 2017; Almeida et al. 2017). Their evolution is therefore inextricably linked to their binarity (Sana et al. 2012; de Mink et al. 2013; Eldridge & Stanway 2022). The presence of a relatively close companion not only enables mass transfer at various stages of the evolution of a massive star (e.g., Sen et al. 2022, and references therein), but it can also affect the rate of mass loss due to line-driven stellar wind (e.g., MacLeod & Loeb 2020; Vink 2022) or induce high-amplitude tidally excited oscillations (TEOs, Kumar et al. 1995; Fuller 2017; Kołaczek-Szymański et al. 2021; Wrona et al. 2022a). In particular, the intensity of stellar wind has important implications for the evolutionary track of a massive star (e.g., Smith 2014). In this context, the recent findings presented by Farmer et al. (2023) seem particularly important. The authors indicate that massive binary systems with an increased mass-loss rate are more efficient producers of potassium and other elements of similar atomic numbers than the evolutionary channels of single massive stars. Unfortunately, there are still many uncertainties regarding theoretical predictions of mass-loss rates in massive stars due to their clumped, radiation-driven winds, especially when they leave the main sequence and become blue supergiants (SGs; e.g., Vink et al. 2001; Krtička et al. 2021; Björklund et al. 2023; Gormaz-Matamala et al. 2023). Therefore, any opportunity to observe massive binary systems during the ongoing process of envelope stripping due to mass transfer or stellar wind enhanced by the presence of a companion is particularly valuable to verify theoretical studies. In this context, MACHO 80.7443.1718 is of particular interest.

MACHO 80.7443.1718 (V≈ 13.6 mag) is an eccentric (eccentricity e ≈ 0.5) single-lined spectroscopic binary system with an orbital period of 32.83 d (Jayasinghe et al. 2021, hereafter J21), located in the Large Magellanic Cloud (LMC). Discovered as variable in the MAssive Compact Halo Object survey (MACHO, Alcock et al. 1999, 2001), it was classified as an eclipsing binary. It was later reclassified by Jayasinghe et al. (2019, hereafter J19) as an eccentric ellipsoidal variable (EEV) through the analysis of the All-Sky Automated Survey for SuperNovae (ASAS-SN, Kochanek et al. 2017) data. The primary component of this system is a B0 Iae-type SG with an effective temperature of 30 000 K. Both the position of the primary component in the Hertzsprung-Russell (HR) diagram (J21) and the properties of its photometric stochastic variability (Kołaczek-Szymański et al. 2022, hereafter KS22) suggest that the star has already evolved off the main sequence (MS) and entered the Hertzsprung gap. Based on the analysis of time-series spectra, the identification of strong emission in hydrogen Balmer lines, and the presence of forbidden iron and oxygen emission lines, J21 also argued that the primary component of MACHO 80.7443.1718 could be B[e] SG (Lamers et al. 1998; Kraus 2019).

Despite attempts, J21 failed to detect the lines of the secondary component in the composite spectrum. The invisible companion is most likely an O-type dwarf. The system is unusual for at least two reasons. First, the total range of brightness variations in the system is about 0.4 mag, which is extreme for EEVs with MS or post-MS components (J19, their Fig. 2 and Wrona et al. 2022a, their Fig. 11). For this reason, J19 dubbed it an “extreme” heartbeat star – a term occasionally used to describe EEVs (e.g., Shporer et al. 2016). Neither J21 nor KS22 was able to explain such a large range of brightness variations with proximity effects using the standard models implemented in the PHysics Of Eclipsing BinariEs (PHOEBE, Prša et al. 2016; Horvat et al. 2018; Jones et al. 2020; Conroy et al. 2020) code. All these attempts yielded amplitudes at least an order of magnitude smaller than those actually observed in MACHO 80.7443.1718. Secondly, the system shows the presence of multiple TEOs with relatively high amplitudes (J19; J21; KS22). Furthermore, MACHO 80.7443.1718 is the first and so far the only EEV in which significant and relatively abrupt changes in both TEOs amplitudes and frequencies have been discovered (KS22). For simplicity, from this point on we refer to MACHO 80.7443.1718 as ExtEV (meaning “extreme eccentric variable”).

The true nature of the observed extreme brightness changes in ExtEV remains a mystery. J21 has suggested that the “heartbeat effect” seen in its light curve may be related to the presence of a circumstellar disk around the primary component. More recently, MacLeod & Loeb (2023b, herafter ML23) explored the possibility of nonlinear tidal wave breaking (MacLeod et al. 2022) on the surface of the primary component as the source of the system’s unique brightness variations. Using hydrodynamical modeling, the authors were able to obtain close-to-real amplitude only by assuming a very rapid rotation rate of the primary component.

The present study was inspired by the striking similarity in the shape of the light curve of ExtEV and the light curve of HD 38282, a massive eccentric binary system consisting of two Wolf-Rayet (WR) stars with a total mass of about 140 M_⊙ and located in the Tarantula Nebula (Sana et al. 2013; Shenar et al. 2021). Shenar et al. (2021) has shown that very dense stellar winds outflowing from both components of HD 38282 cause atmospheric eclipses when the light emitted by the more distant component is attenuated by the wind of the star that is closer to the observer (Auer & Koenigsberger 1994; Lamontagne et al. 1996; St-Louis et al. 2005). In addition, the colliding winds form a turbulent and shocked wind-wind collision (WWC) cone (Parkin & Gosset 2011; Teodoro et al. 2012; Kashi & Michaelis 2021), which may be responsible for the excess emission near the periastron passage. These combined effects can result in a light curve that mimics real EEVs, even those that undergo nonlinear tidal wave breaking on their surface. We investigate in our paper whether similar phenomena may be at work in the much less massive ExtEV system, whose total mass lies between ∼35 M_⊙ and ∼60 M_⊙, and explain the origin of this system’s unique light curve.

The paper is organized as follows. In Sect. 2, we describe the origin and preparation of the photometric and spectroscopic data used in our study. In the next section, we analyze the spectral energy distribution (SED) of the ExtEV. In Sect. 4, we discuss various properties of the ExtEV, including its potential status as a B[e] SG, the initial and current mass of the primary component, and the orbital parameters. Based on the information from the previous section, we present in Sect. 5 the modeling of the ExtEV light curve using an analytical model that accounts for phenomena typical of WWC binary systems. We discuss selected aspects of our results in Sect. 6, while a summary and final conclusions are provided in Sect. 7.

2. Observations

2.1. Time-series photometry

ExtEV is located in the continuous viewing zone of the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015) space mission. It was therefore observed during Year 1 (July 2018–July 2019), Year 3 (July 2020–July 2021) of TESS operation. Presently, the observations of the southern hemisphere scheduled for Year 5 (January–September 2023) are being completed. In the previous analysis of the TESS data of ExtEV (KS22), we performed photometry of the star on the TESS Full Frame Images obtained during Year 1 and most of the Year 3 observations. In the present paper, we use the same data as they are sufficient for our purposes. The TESS light curve phased with the orbital period is shown in Fig. 1. All details of its extraction can be found in KS22.

Fig. 1. TESS light curve of ExtEV phased with the orbital period of 32.83016 d. Orbital phase φorb = 0.0 corresponds to the time of the periastron passage derived in Sect. 4.3. The inset shows the light curve near the periastron passage. Stochastic light variations and TEOs can be seen over the entire range of orbital phases.

For the reasons outlined in Sect. 4.4, it is important to check the range of variation and the shape of the light curve of ExtEV in short-wavelength passbands. Unfortunately, we did not find any publicly available light curves in such bands, so we collected near-simultaneous time-series photometry in the Johnson U, B and V passbands at three observatories of the Las Cumbres Observatory Global Telescope (LCOGT) network, namely Siding Spring Observatory, Cerro Tololo Inter-American Observatory, and South African Astronomical Observatory. The images of the ExtEV were taken by a total of eight 1-m Ritchey-Chrétien telescopes, each equipped with the Sinistro CCD Camera¹. The field of view of the camera is 26.5′  ×  26.5′, with an angular resolution of 0.389″ per pixel. The LCOGT data cover an interval of about 68 days from 9 February to 18 April 2022, which was enough to cover just over two orbital cycles of the system. The observations were spread over the whole range of orbital phases with the desired higher density close to periastron passages.

We perform standard point spread function (PSF) photometry on the collected images using the tools included in photutils Python package² (Bradley et al. 2022). By modeling and subtracting the nonuniform background of the LMC, we managed to construct empirical models of the PSF following the procedure described in Anderson & King (2000), which we then fit using least-squares to the actual profiles of the ExtEV and several nonvariable comparison stars in its vicinity. The resulting light curves are illustrated in Fig. 2.

Fig. 2. LCOGT light curves of ExtEV in Johnson U (purple dots), B (blue dots), and V (green dots) passbands. Panel a shows the light curves. A ±0.25 mag vertical offset was applied to U and V-filter light curves to avoid overlapping. Panels b–d show the light curves phased with the orbital period. Phase zero corresponds to the time of periastron passage. In addition, we have drawn dash-dotted lines in each of the lower panels to indicate the range of variability during the heartbeat. The plot illustrates that the range is virtually the same in all three passbands.

2.2. Time-series optical spectroscopy and radial velocities

Using the High-Resolution Spectrograph (HRS, Crause et al. 2014) installed on the Southern African Large Telescope (SALT, Romero Colmenero et al. 2022), we collected a series of 18 optical spectra in low-resolution mode with an average resolving power of 14 000 between 1 October 2021 and 16 March 2022. All SALT/HRS fiber-fed echelle spectra were calibrated and reduced by the HRS MIDAS pipeline³ (Kniazev et al. 2016, 2017). The average signal-to-noise ratio (S/N) in the vicinity of the Balmer Hα and Hβ lines was around 80 and 50, respectively. Due to the poor quality of two spectra, we ultimately utilized 16 of them in our analysis and normalized them with the automatic tool SUPPNet (Różański et al. 2022).

In the SALT/HRS spectrograph, in addition to the fiber centered on the target, there is a second fiber nearby which is usually used to measure the background. The spectrum measured with this fiber is later used to remove telluric lines from the spectrum of the target. Unfortunately, the ExtEV is located in a crowded field, where an extended and spatially variable emission nebula is present. This makes it difficult to find a location free of any astrophysical sources. We therefore decided to use in our analysis only the spectra from the ExtEV-centered fiber, without correcting for the signal measured by the second fiber, as it contained signal from the emission nebula and nearby fainter stars. We do not fit and remove telluric lines from any of our spectra.

We derive the radial-velocity (RV) curve for ExtEV in the following way. We first selected the spectrum with the highest S/N as the reference spectrum. We then determined the RVs of the remaining spectra relative to the reference spectrum by cross-correlating each of them with the reference spectrum. The cross-correlation function was calculated only in certain wavelength intervals centered on strong absorption lines showing no emission (Table 1). In practice, to obtain RV, we fit a parabola to the cross-correlation function around its global maximum. Then, the position of the peak of the fit parabola was taken as the measured RV. Knowing these RVs, we Doppler-shifted all spectra so that they could be averaged into a single high-S/N template. We iteratively repeated the procedure described above three times. A newly obtained template served each time as a new reference spectrum. Finally, the RV zero point was calculated by cross-correlating the final reference spectrum with a synthetic spectrum from the BSTAR grid of model stellar atmospheres (Lanz & Hubeny 2007) with effective temperature T_eff = 30 000 K, surface gravity log(g/(cm s⁻²)) = 3.25 dex, and half-solar metallicity, convolved with the rotational profile with the projected rotational velocity of 175 km s⁻¹ (J21). The resulting RVs, together with the estimated formal errors, are presented in Table 2.

3. SED fitting

Although J21 performed an SED fit to the selected photometry available for ExtEV, we decided to repeat this task because of two reasons. Firstly, as we argue below, ExtEV shows an infrared (IR) excess, and J21 did not exclude the near-IR (NIR) bands from their fit, which may have biased the result. Secondly, due to the high effective temperature of the primary component, every single measurement in the ultraviolet (UV) is particularly important. For ExtEV, we found several sources of such measurements in the literature (Table 3). We also utilize the fact that the system has been observed by the Galaxy Evolution Explorer (GALEX, Morrissey et al. 2007) satellite in the near-UV (NUV) band, although GALEX photometry for ExtEV has not been published. We therefore extract our own photometry from the “tile” in the NUV passband from the GALEX GR6/7 data release (Bianchi et al. 2014) as follows. We downloaded a tile named “LMC_39774_0422” from the MAST portal⁴, on which we performed circular aperture photometry of ExtEV using the photutils tool. We estimated the mean background signal at the ExtEV position using a nearby region without any point sources. We then converted the measured signal to magnitude using a tool provided at GALEX website⁵. A full list of ExtEV magnitudes with references can be found in Table 3.

To estimate the radius of the primary component, R₁, and consequently its luminosity, L₁, we perform the following model fitting to the observed SED. First, we converted all magnitudes presented in Table 3 to absolute spectral flux densities using the VOSA 7.5 online tool⁶ (Bayo et al. 2008). If a given photometric passband had more than one measurement in Table 3, we averaged the corresponding fluxes. Due to the IR excess of ExtEV, we considered only passbands with effective wavelengths shorter than the J filter. We then assumed that the distance to the ExtEV is equal to the distance to the LMC (49.59 kpc, Pietrzyński et al. 2019) and that the absolute spectrum of the primary component corresponds to the ATLAS9 model atmosphere (Castelli & Kurucz 2003) with the parameters T_eff = 30 000 K, log(g/(cm s⁻²)) = 3.5 dex, [Fe/H] = −0.5 dex, and solar-scaled chemical composition (Fig. 3, black line). The first two parameters were obtained by J21, while the last one is close to the mean metallicity of the LMC (which is about −0.42 dex, Choudhury et al. 2021). Based on the adopted synthetic spectrum and the extinction law of the interstellar medium (ISM) of Fitzpatrick (1999), implemented in the extinction Python package (Barbary 2016), we were able to calculate the theoretical SED with the sedpy tool (Johnson 2021). We perform a Bayesian fit of the theoretical values to the observed ones, using the standard maximum-likelihood estimation method. However, we took the so-called noise nuisance parameter, ξ, into account when calculating the total flux density errors, so they were given by ${\sigma }_{\text{total}}^2={\sigma }_{\text{phot}}^2+{\mathcal{M}}^2\exp (2\ln \xi )$, where σ_phot refers to the reported photometric error and ℳ stands for the model SED value. The purpose of this procedure is twofold. Firstly, photometric errors are usually underestimated. Secondly, in addition to the heartbeat-like feature, ExtEV also exhibits TEOs and stochastic variations in its light curve, resulting in some intrinsic scatter in the photometric measurements⁷ collected in Table 3. When fitting the theoretical SED, the following parameters were considered free: the total extinction in the V band, A_V, the ratio of total to selective extinction, R_V, R₁, and lnξ. Once the optimal solution was identified, we estimated the errors of the fit parameters using the affine invariant Markov chain Monte Carlo (MCMC) method (Goodman & Weare 2010), using the implementation provided by the emcee package (Foreman-Mackey et al. 2013).

Fig. 3. Results of the SED fitting to the multiband photometry of ExtEV. The upper panel shows the absolute fluxes calculated from magnitudes given in Table 3 (dots are color-coded according to the legend). The model spectrum of the primary component is shown with a black line, its reddened versions with a red line (with additional attenuation) and a blue line (without it). The infrared excess is clearly visible and indicated by the shaded area. The middle and lower panels show the residuals from the best fit, with and without the additional attenuation, respectively. For residuals presented on the left (UV and optical range), the ordinate axes are linearly scaled while the residuals in the IR range are normalized by the model values and logarithmically scaled.

Interestingly, fitting the SED using the extinction law of Fitzpatrick (1999) results in “U”-shaped residuals in the UV region and the U and B passbands (Fig. 3, bottom panel) regardless of the values of R_V and A_V. We also find that this situation cannot be avoided by choosing different prescriptions for the ISM extinction function. This behavior of the residuals suggests the presence of an additional source of attenuation that is not accounted for by the ISM extinction laws. Therefore, to consider the possible existence of an additional extinction component, we multiplied the interstellar extinction law by the following Gaussian term, $1-{\alpha }_{\text{G}}\exp [{(\lambda -{\lambda }_{\text{G}})}^2/2{\sigma }_{\text{G}}^2]$, where α_G, λ_G, and σ_G stand for the amplitude, mean, and standard deviation of the Gaussian profile, respectively (Fig. 4). The effects of this fitting procedure are shown in Fig. 3, while the optimized parameters can be found in Table 4.

Fig. 4. Extinction law resulting from the best fit to the observed SED of ExtEV. The black curve refers to the interstellar dust extinction function of Fitzpatrick (1999) with the parameters given in the legend. The additional extinction feature is shown with the gray line, while the red dashed curve corresponds to the total extinction function (the product of black and gray curves). With a pair of blue vertical dotted lines, we have marked the positions of the two local minima in the total extinction function. More details can be found in the main text.

The best-fit value of R₁ is 30 R_⊙, which is significantly higher than the value obtained by J21 (their R₁ ≈ 24 R_⊙). This also implies that the value of L₁ is higher than that obtained by J21, amounting to about log(L₁/L_⊙) = 5.82. The location of the ExtEV’s primary component in the HR diagram is subsequently discussed in Sect. 4.2. The rather unusual value of R_V = 1.6 and the additional attenuation component in the total extinction function can be attributed to the presence of dust in the vicinity of the ExtEV with grain sizes and chemical composition significantly different from those characterizing the Galactic ISM (e.g., Nozawa 2016). Since polycyclic aromatic hydrocarbon (PAH) molecules are responsible for the well-known absorption maximum around 2175 Å (e.g., Blasberger et al. 2017), depending on the conditions of the medium and their mixture, they can give rise to additional relatively broad absorption bands in the UV range (Malloci et al. 2004; Steglich et al. 2011). Perhaps the additional attenuation factor seen in the ExtEV SED comes from PAHs or other molecules and mineral grains not typical of the ISM.

4. Properties of the ExtEV system and its primary component

4.1. Non-B[e] status of ExtEV

B[e] supergiants are characterized by the presence of two components in their intense stellar wind: a “fast wind” located around the poles and a “slow wind” in the form of an outflowing disk emanating at a certain angle from the equatorial region (e.g., Kraus et al. 2007). In addition to the strong emission in the hydrogen lines, the primary distinguishing feature of B[e] SGs is the presence of forbidden [O I] lines at 5577, 6300, and 6363 Å (e.g., Lamers et al. 1998), which form in the immediate vicinity of the star. Other distinguishing features include emission in the Fe II and [Fe II] lines, as well as strong molecular CO emission and a significant NIR excess, which originates from relatively hot circumstellar dust with temperatures on the order of 1000 K (Kraus 2019).

Although J21 reported the presence of some of the aforementioned emission lines in their spectra and claimed that ExtEV is a B[e] SG with a disk, this emission disappeared in some of their spectra (their Fig. 16). In none of our spectra from SALT/HRS, did we observe the presence of the emission lines [O I], Fe II, and [Fe II]. This may suggest that the detections made by J21 were unlikely to involve ExtEV. Some of these could come, for example, from cosmic rays or inadequate removal of nebular and geocoronal emission lines. The absence of the three strongest [O I] lines in our spectra is presented in Fig. 5. For reasons explained in Sect. 2.2, only [O I] lines resulting from geocoronal emission are present at v_rad = 0 km s⁻¹. In the vicinity of v_rad = 313 km s⁻¹, the systemic velocity of the ExtEV (Sect. 4.3), no significant [O I] emission is observed over a period of about six months covering different orbital phases. Hence, we would like to conclude that the SALT/HRS spectra of ExtEV do not show the typical spectral features of B[e] SGs.

Fig. 5. Small sections of the SALT spectra showing [O I] emission lines at their laboratory wavelengths but not in the rest frame of ExtEV. The rest-frame wavelengths of the lines shown are indicated in each panel in the top left corner. The spectra are shifted vertically according to the orbital phase. The vertical red and blue dashed lines correspond to RVs of 0 (geocoronal emission) and 313 km s−1 (systemic velocity of ExtEV), respectively. For comparison, the bottom right panel shows the He I absorption line originating from the photosphere of the primary component of ExtEV.

Analysis of the photometric data presented in Sect. 3 indicates that ExtEV shows a strong IR excess, particularly in the mid-IR (MIR) range (W3 and W4 WISE passbands). Therefore, we decided to investigate the infrared properties of ExtEV by reproducing the NIR and WISE color-color diagrams from Kraus (2019, their Fig. 6). The result of this analysis is shown in Fig. 6. It becomes apparent that ExtEV has markedly different IR properties compared to B[e] SGs and instead falls within the range typical of luminous blue variables (LBVs, e.g., Smith 2017). This implies that there is no hot dust in the vicinity of ExtEV, which usually accumulates relatively close to the central SG. On the contrary, the clear excess in the MIR (Fig. 6b) suggests the presence of cooler dust, probably surrounding the entire binary system at a considerable distance from both components. In the case of SGs, the (W3 − W4) – (K_s − W3) color-color diagram may serve as an indirect indicator of the circumstellar dust temperature and mass loss rate, respectively (González-Fernández et al. 2015, their Fig. 15). In the case of ExtEV, (K_s − W3)≈2.8 mag, indicating that the system is undergoing enhanced mass loss (Fig. 6c), a feature common to both LBVs and B[e] SGs.

Fig. 6. Infrared color-color diagrams for ExtEV. The position of ExtEV is marked by the cyan star in each panel. The blue dots and red triangles correspond to the B[e] SGs and LBVs located in the Magellanic Clouds, respectively. Their photometric properties were adopted from Kraus (2019). Black crosses refer to the sample of late-type SGs selected by González-Fernández et al. (2015). The shaded region in panel (c) shows the area with color index (Ks − W3)≳2 mag, where SGs with enhanced mass loss are suspected to be located.

Given all of the above, we can rule out the possibility that ExtEV is a B[e] SG. This conclusion is further strengthened by the analysis presented in Sect. 4.5. Whatever the conditions in the vicinity of the primary component of ExtEV are, with the exception of the increased rate of mass loss, they differ from those typical of B[e] SGs.

4.2. Estimating the mass of the primary component

Having estimated the effective temperature of the primary component of ExtEV (J21) and its luminosity (Sect. 3), we were able to locate it in the HR diagram (Fig. 7). In this way, we could estimate its initial mass at zero-age MS (ZAMS), M_1, ZAMS, and its current mass, M₁, by examining the evolutionary tracks consistent with the location of the primary component in the HR diagram. To do so, we used the Modules for Experiments in Stellar Astrophysics (MESA; version r22.11.1, Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023) stellar structure and evolution code, compiled with the MESA Software Development Kit (version 21.4.1, Townsend 2021). Since we were not interested in the precise modeling of the evolution of ExtEV as a binary system, but only in the range of probable M_1, ZAMS and M₁, we synthesized a series of 1000 single-star evolutionary tracks that could potentially intersect the current position of primary component in the HR diagram within 1-σ error bars. The evolution of massive stars is very sensitive not only to M_ZAMS and its metallicity, but also to the rotation velocity divided by its critical value Ω_rot/Ω_crit in the ZAMS and the mass-loss rate due to the line-driven stellar wind, Ṁ_wind. In order to realistically capture the richness of possible evolutionary tracks, we computed them with their initial parameters drawn from the distributions presented in Table 5. While the ranges of the first four parameters presented in this table require no further comment, the last parameter necessitates a deeper justification. For massive stars, the most commonly used formula for Ṁ_wind is the one given by Vink et al. (2001). Our calculations in MESA also used this formula, but we noted that it did not account for any effects related to binarity. On the one hand, there is a growing evidence (e.g., Björklund et al. 2023) that the aforementioned analytical prescription for Ṁ_wind may in some cases overestimate the mass-loss rate by a factor of 2–3. Hence, we set the lower limit of S_wind to 0.5 to account for the situation where the primary component evolved effectively as a single star. On the other hand, this study provides strong evidence indicating an enhanced mass-loss rate that is ongoing in ExtEV. Given the current state of the primary component, it is not excluded that in the past Ṁ_wind was also larger than typical (e.g., due to the presence of a close massive companion). We have therefore allowed for situations in which the actual Ṁ_wind may be up to five times more intense than predicted by the Vink et al. (2001) formula. The remaining MESA parameters and the set of “physics switches” were identical to those described by Kołaczek-Szymański & Różański (2023, their Sect. 3.2.2).

Fig. 7. HR diagrams summarizing our simulations described in Sect. 4.2. Black and blue dots show the location of the primary component of ExtEV, according to this paper and the study of J21, respectively. The gray curves show evolutionary tracks that never approach the position of the primary component determined in our study. Conversely, multicolored tracks correspond to those that are within the error box of the ExtEV (black dot) at some phase of the evolution. The color bar above the HR diagrams reflects M1, ZAMS (left panel) and M1 (right panel) and is common to both panels.

Looking at the left panel of Fig. 7 one can see that the initial mass of the primary component is most likely in the range between 27 and 55 M_⊙. A more precise value cannot be given because it would be necessary to know the past evolution of ExtEV. The current value of M₁ is a separate issue, given the strong stellar wind in this mass range. The right panel of Fig. 7 suggests that the primary component has currently a mass between 25 and 45 M_⊙. For simplicity, we assume M₁ = 35 M_⊙ in the subsequent sections of this paper and discuss the potential impact of changing this value on our findings. It is also difficult to make a definitive conclusion about the evolutionary status of the ExtEV, as our experiment indicates that the primary component may either be still at the MS, close to the terminal-age MS (TAMS), or has already ignited helium in its core. So it may be a post-MS star that crosses the Hertzsprung gap.

4.3. Spectroscopic orbit

In order to obtain spectroscopic parameters of the primary’s orbit, we used the RV data described in Sect. 2.2 together with those reported by J21⁸. The data were fit with RadVel software (Fulton et al. 2018), which uses a maximum a posteriori optimization method to find the best orbital solution and MCMC sampling to estimate errors. During the fitting, the orbital period, P_orb, was treated as a fixed parameter, while the argument of periastron, ω, the half-range of the RV changes of the primary component, K₁, and the time of the periastron passage, T_peri, were treated as free parameters. Given that we use RVs obtained with three different spectrographs (SALT/HRS, MIKE, and SOAR), during the modeling procedure we allowed the corresponding systemic velocities (γ_HRS, γ_MIKE, and γ_SOAR) to be independent. Furthermore, we noticed that the residuals from the initial fit revealed the presence of a long-term trend. We fit this trend by adding a parabolic term that takes into account both the linear ($\dot{\gamma }$) and quadratic ($\ddot{\gamma }$) changes of the systemic velocity. Due to the significant dispersion of RVs, which is not due to measurement errors, but originates from the intrinsic variability of the primary component (TEOs and stochastic variability), we also added the so-called jitter term (J_HRS, J_MIKE and J_SOAR) as an additional parameter fit separately for each dataset. This parameter accounts for the underestimated error associated with each RV measurement. As both J_MIKE and J_SOAR were consistent with zero within the errors⁹, we finally set their values to zero and fit only J_HRS. The results of the fit are given in Table 6 and shown in Fig. 8.

Fig. 8. Optimized spectroscopic solution for the RVs of the primary component of ExtEV. The points on all panels are marked with different colors depending on the spectrograph (labeled in the top panel). The blue curve represents the best-fit model. (a) RV curve including orbital solution and the long-term changes. (b) RV curve with the orbital solution subtracted. (c) RVs phased with the orbital period given in Table 6 after freeing from the long-term changes of RVs shown in panel b.

Based on the obtained solution, we calculated the mass function for this system, which amounts to f(M) = 0.74 ± 0.05 M_⊙. This value is significantly larger than $0.{51}_{-0.06}^{+0.07}$ M_⊙ given by J21. We also emphasize that the quadratic term is nonzero at the level of 2σ, implying that the ExtEV may be indeed a triple hierarchical system, with a tertiary component in a wide orbit with a relatively long orbital period. However, similar to the secondary component, the tertiary remains undetected in the ExtEV spectra.

4.4. Geometry of the system during periastron passage

Theoretically, the large amplitude of light changes observed in ExtEV could be caused by the extreme nonlinear tidal distortion of the primary component. This would mean a situation in which, when passing through the periastron, the primary component would be close to or filling its Roche lobe. This is exactly the scenario that was investigated by ML23 through hydrodynamic simulations. However, the number of assumptions made by these authors, for example, the lack of radiative losses, the adopted parameters of the system selected in the way that artificially increases tidal interaction, and the highly approximate method of transition from the hydrodynamic model to the light curve, do not allow this solution to be considered certain. Moreover, the authors claimed that the deep minimum in the light curve of ExtEV originates from a tidal bulge lifted almost directly toward the observer during the periastron passage. This conclusion was the result of a mistake that was made by using ω that was 180° different from the true value. The mistake has been corrected in erratum (MacLeod & Loeb 2023a), but the updated figures reveal a smaller amplitude of light changes than originally announced by ML23 and the huge tidal bulge is no longer pointing toward the observer near the periastron.

ML23 adopted most of the system parameters from J21, including the rotation period of the primary component of about 4.4 d. It turned out that only with an exceptionally fast rotation of the primary component, that is, significantly larger than the pseudo-synchronous rotation rate, can Roche lobe overflow (RLOF) be observed in their simulations, followed by a subsequent nonlinear collapse of a massive tidal wave generated by this close approach. However, KS22 (their Sect. 2.4) have already pointed out that there is no reasonable basis for the claim that the primary component of ExtEV rotates with a period of 4.4 d. This value is the result of misclassifying one of the “red-noise” maxima in the frequency spectrum of the TESS data of ExtEV as the rotational frequency. This is evidenced by the lack of a significant maximum around frequency corresponding to the period of 4.4 d in the frequency spectrum calculated by KS22 (their Figs. 4 and 10) using a much longer TESS light curve than that used by J21. We demonstrate in Sect. 6.2 that the rotation period of the primary component amounts to at least about 8 d. Next, the model presented by ML23 predicts that the frequency spectrum of the light curve should be dominated by the signal corresponding to the rotational frequency, which is not present in the TESS data. Finally, their hydrodynamic model is unable to reproduce one of the most characteristic features of the ExtEV light curve, that is, the TEO which is the 25th harmonic of the orbital frequency. All this makes the ML23 model questionable.

From an observational point of view, a very strong tidal distortion leading to the RLOF of the primary component near the periastron should induce a dramatic gravitational darkening (von Zeipel 1924; Prša 2018) of its outermost layers. Hence, one can expect heterochromatic effects in multicolor photometry, which should reveal the dependence of the amplitude of the ExtEV light curve on wavelength. Thanks to the LCOGT photometry (Fig. 2), and the TESS photometry (Fig. 1), we can confidently conclude that the amplitude and shape of the heartbeat-like feature remains unchanged from Johnson U to NIR TESS band at approximately 0.02 mag level. This rather excludes the occurrence of strong temperature changes (on the order of a few 10³ K) on the surface of the primary component. This conclusion is further supported by the lack of depth changes in the He I lines with the orbital phase (Fig. 5), which we examined in our SALT/HRS spectra. The helium lines are very sensitive to changes in gas temperature (Gray & Corbally 2009), so any strong tidal distortions on the surface of the primary component should translate into significant changes in the depths of He I and He II lines. The depths of ionized metal lines (for example, the absorption lines of C III) do not change significantly in our spectra either. The phenomenon of light changes in ExtEV appears to be gray, inconsistent with the extreme tidal distortion and nonlinear breaking of tidal waves proposed by ML23. We note, however, that for relatively small ellipsoidal distortion of the primary component during periastron passages, we do not expect to detect any heterochromatic effects in our multicolor light curves, as both the Johnson U and TESS passbands lie on the Rayleigh-Jeans tail of its spectrum.

To check whether the primary component may undergo RLOF at periastron, we used the method presented by Sepinsky et al. (2007) to calculate its volume-equivalent Roche radius at periastron, R_L¹⁰. Unlike the Eggleton (1983) formula, the method of Sepinsky et al. (2007) takes into account the nonzero eccentricity of the orbit. We use the quantity ϑ = (R_L − R₁)/R₁ as the RLOF indicator. A negative value of this quantity indicates RLOF at the periastron, while a positive value corresponds to the situation in which the system maintains a detached geometry throughout the entire orbital cycle.

Since we do not know precisely all the parameters of the ExtEV system, it is necessary to perform a statistical simulation to understand the distribution of the quantity ϑ. The initial parameters of our simulation are listed in Table 7. Parameters such as M₁, R₁, f(M), and e and their ranges were adopted according to the results presented above. Since the inclination, i, is unknown, we assumed a realistic range of 30° to 75°. Lower inclinations would correspond to unrealistically large masses of the primary component, while with higher values, we could expect photospheric eclipses, which are not observed. The last parameter in Table 7, ℱ_rot, stands for the ratio of the rotational angular velocity of the primary component and its orbital angular velocity at periastron. The range adopted by us for ℱ_rot allows for rotation rates lower than pseudo-synchronous. The results are presented in the form of a histogram in Fig. 9. As can be seen, only ∼3% of our models are characterized by RLOF at periastron. The majority of them (∼97%) remain detached throughout the entire range of orbital phases. The most probable value of ϑ amounts to approximately 0.23, suggesting that the primary component likely fills only about 50% of its Roche lobe volume at periastron.

Fig. 9. Histogram of the distribution of parameter ϑ = (RL − R1)/R1, resulting from the simulations described in Sect. 4.4. The vertical dashed black line indicates ϑ = 0, when the system becomes semi-detached. The vertical blue lines indicate quantiles at levels of 0.16, 0.5, and 0.84. Only about 3% of samples correspond to an RLOF at periastron, while their majority (∼97%) points toward the detached geometry at periastron.

4.5. No disk around the primary component

The possibility of a periodically fading and reemerging decretion disk in the ExtEV system was strongly suggested by J21 to explain the double-peaked Balmer emission lines. These are particularly visible in the Hα and Hβ lines (Fig. 10) in the form of broad, disjoint, and asymmetric wings accompanying the narrow emission peaks that appear in our spectra due to the nebular emission of the LMC. J21 concluded that such a disk must destabilize and disappear near the periastron passage, but reappear shortly thereafter, to generate such features throughout almost the entire orbital period. However, as can be easily seen in Fig. 10 (lower panel), our spectra contradict this theory. There is no doubt that the emission in the Hβ line persists even immediately after periastron passage and is not significantly weaker than just before phase zero.

Fig. 10. Fragments of SALT/HRS spectra of ExtEV near the Hα (upper panel) and Hβ (lower panel) lines. The spectra are color-coded according to the orbital phase (the color bar is located at the top of the figure), with zero phase corresponding to the periastron passage. The narrow emission in the center of the spectrum is of nebular origin.

To further investigate the plausibility of the scenario in which the emission features of the Hα and Hβ lines are attributed exclusively to the Keplerian disk around the primary component, we carried out simplified modeling in which we tried to estimate the actual location of the circumstellar disk and its geometry necessary to reproduce observable shapes of spectral lines. The aim of this experiment is not to determine the exact parameters of the disk, but to try to answer the question of whether it is crossed by the orbit of the secondary component, assuming that the two emitting components in the Balmer lines originate from the disk. We based our model on the following assumptions: (1) the disk forms a uniform, optically thin ring around the primary component, (2) it is aligned with the orbital plane, and (3) the rotation of the disk is Keplerian, so the velocity of the gaseous component at a distance r_el from the primary can be expressed as

$$\begin{array}{c}\hfill v_{\mathrm{el}}=\sqrt{\frac{G{M}_1}{r_{\mathrm{el}}}},\end{array}$$(1)

where G is the universal gravitational constant. To obtain synthetic emission profiles, we associated each surface element of the disk with an identical, constant emitter of either Hα or Hβ photons with fixed wavelengths corrected for the radial velocity of the primary component. The resulting orbital velocity field – uniquely defined by the width of the disk, w_D, and the distance of its inner edge from the primary component, d_D – was then projected into the radial direction and used to obtain the density distribution ℋ(λ,d_D,w_D) via relation:

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{obs}}\simeq {\lambda }_{\mathrm{em}}(1+\frac{v_{\mathrm{rad}}}{c}),\end{array}$$(2)

which is a nonrelativistic approximation of the Doppler formula. Here, λ_em is the emitted wavelength, λ_obs is the observed one, v_rad denotes the radial component of the element’s velocity, and c is the speed of light in vacuum.

To develop a model that could be directly fit to the ExtEV spectra, the resulting probability distribution ℋ(λ,d_D,w_D) had to be divided into two separate blue- and red-shifted distributions, namely ℋ_B and ℋ_R, and later convolved with the intrinsic emission line profile. For this purpose, we chose the Voigt function

$$\begin{array}{c}\hfill \mathcal{V}(\lambda ,{\sigma }_G,{\gamma }_L)={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\mathcal{G}({\lambda }^{\prime },{\sigma }_G)\mathcal{L}(\lambda -{\lambda }^{\prime },{\gamma }_L)\mathrm{d}{\lambda }^{\prime },}\end{array}$$(3)

where 𝒢 and ℒ denote the Gaussian and Lorentzian profiles, respectively. Their widths, parameterized by the standard deviation σ_G and the scale factor γ_L, are used to account for thermal and pressure broadening, as well as any possible small-scale inhomogeneities or turbulent motions of the gas. The final model takes the following form

$$\begin{array}{cc}\hfill \mathcal{M}& =A_{\mathrm{B}}{\mathcal{H}}_{\mathrm{B}}(\lambda ,d_{\mathrm{D}},w_{\mathrm{D}})\ast \mathcal{V}(\lambda ,{\sigma }_G,{\gamma }_L)\hfill \\ \hfill & +A_{\mathrm{R}}{\mathcal{H}}_{\mathrm{R}}(\lambda ,d_{\mathrm{D}},w_{\mathrm{D}})\ast \mathcal{V}(\lambda ,{\sigma }_G,{\gamma }_L),\hfill \end{array}$$(4)

and is described by six parameters: d_D, w_D, σ_G, γ_L, and two additional scaling factors, A_B and A_R, which take into account the heights of both peaks, separately the blue-shifted and red-shifted. We fit this model to the SALT spectra using Monte Carlo sampling, based on a random walk in the parameter space with wide, uniform priors. We obtained the best-fitting profile by simply selecting the one that minimizes the χ² discriminant. The results are shown in Fig. 11, in which we compared three spectra at different orbital phases with fit theoretical line profiles, along with visualizations of the sizes and locations of the corresponding rings. As can be seen, our experiment results in a disk that is intersected by the orbit of the secondary component. In such a configuration, it is difficult to imagine a disk that can form around the primary component and survive periodic and permanent gravitational perturbations from the secondary component. Even if such a disk were indeed not aligned with the orbital plane, the companion star should still intersect it, effectively destabilizing its structure. We therefore conclude that the emission features in ExtEV are unlikely to be due to the presence of a Keplerian disk around the primary. The observed shapes of the Hα and Hβ lines must originate from another source.

Fig. 11. Models of the Keplerian rings around the primary component of ExtEV. Each set of figures (a, b, and c) corresponds to a different orbital phase φorb provided above each set. The right-hand panels show the observed Hα (salmon) and Hβ (aqua) lines to which the theoretical profiles (black curves) were fit. For each spectrum, two rings were fit separately, one responsible only for the Hα emission and the other for the Hβ emission. The narrow nebular emission of the LMC (gray region bounded by vertical dashed lines) was excluded from the fit. The left-hand panels show the resulting ring geometry compared to the size of the primary component (black circle) and the size of the orbit of the secondary component (dashed ellipse, secondary component size is not to scale). More details and discussion can be found in the main text.

5. Model of light variability

In the previous sections of this paper, we have argued that the nonlinear model of tidal interactions presented by ML23 has serious shortcomings that question its ability to reliably explain the variability of ExtEV (Sect. 4.4). We have also presented arguments against the presence of a disk around the primary component (Sects. 4.1 and 4.5). However, there are several indications of a significant mass loss in this system, which may be the direct cause of large-amplitude changes of the brightness observed in ExtEV. We therefore investigated an alternative model of the optical variability of ExtEV based on a high wind-driven mass-loss rate, enhanced by the presence of a nearby massive companion and its periodic periastron passages.

We decided to check whether the variability of ExtEV can be explained by a model analogous to the one that explains the photometric variability of HD 38282 system (Shenar et al. 2021), which is characterized by e = 0.506 and ω = 304.6°, very similar to those of the ExtEV’s orbit. This system and ExtEV show surprisingly similar light curves. Figure 12 presents the dimensions and geometry of the ExtEV system to scale, along with three characteristic points in the secondary’s orbit, labeled A (superior conjunction), B (periastron passage), and C (inferior conjunction). In our scenario, the decrease in brightness in the light curve of ExtEV has its origin in an atmospheric eclipse, when the light from the secondary component and the WWC cone is obscured by the dense stellar wind emitted by the blue SG (Fig. 12, phase A). Both components are sources of stellar winds with completely different properties. As a B-type SG, the primary component emits an intense but relatively slow stellar wind that collides with the fast but low-density wind of the secondary component, typical for O-type dwarfs. As a result of the collision of the two winds, a geometrically thin WWC cone is formed very close to the photosphere of the secondary component. There is also a chance that the wind from the primary component is strong enough compared to the wind from the secondary that it will cause the front of the WWC cone, particularly at periastron, to collapse directly onto the surface of the secondary component, further eroding its photosphere. The gas in the WWC cone is so hot that it can be considered completely ionized, at least in the vicinity of the secondary component. The free electrons present in it can effectively scatter the light emitted from both components and direct it toward the observer. Additionally, hot gas flowing along the WWC cone eventually cools and recombines, emitting light in specific emission lines. Phase B in Fig. 12 corresponds to the situation when the WWC cone is the densest and hottest at the time of periastron passage. This is the source of the excess emission observed in ExtEV. After passing the periastron, in phase C, the secondary component is partially obscured by the dense and cooler parts the WWC cone. In fact, we expect the WWC cone to be significantly curved near periastron passage due to the relatively fast orbital motion of the secondary compared to the velocity of stellar wind.

Fig. 12. Artistic representation of the orbit (silver ellipse) of the secondary component of ExtEV (dark blue circle) around the primary component (light blue circle near the center). The size of the orbit and both components are to scale. The orientation of the orbit corresponds to the view of the ExtEV in the sky (i.e., as seen by the observer on Earth). The primary component is the source of a slow but dense stellar wind (semi-transparent spheres, concentric with the primary component), which, after colliding with the surface of the wind of the secondary component, forms a WWC cone, a turbulent structure, with its apex located near the secondary component. Three orbital phases are labeled, corresponding to: the superior conjunction (A), periastron passage (B), and inferior conjunction (C). More details can be found in the main text.

5.1. Analytical model

The total brightness of ExtEV as a function of the true anomaly, ν, can be approximated as follows:

$$\begin{array}{c}\hfill m\approx -2.5log\{F_1+[F_2+F_{\mathrm{WWC}}(\nu )]e^{-\tau (\nu )}\}+C^{\prime },\end{array}$$(5)

where F₁, F₂, and F_WWC denote the flux from the primary component, from the secondary component, and from the WWC cone, respectively. The τ(ν) denotes the optical depth in the wind of the primary component, integrated along the line of sight from the observer to the secondary component. C′ is the zero point of the magnitude scale. Equation (5) assumes that the size of the secondary component is much smaller than that of the primary component. (According to our results, R₁ is about six times larger than R₂, which allows us to make this assumption.) Moreover, we also assumed that the brightest part of the WWC cone is in the close proximity of the secondary component. For this reason, the F_WWC in our model is also affected by the factor e^−τ because the excess flux emitted from the WWC cone can be obscured by the primary’s wind. The relative change in the total brightness of ExtEV can be obtained from Eq. (5) as

$$\begin{array}{c}\hfill \mathrm{\Delta }m\approx -2.5log\{1+[(F_2/F_1)+(F_{\mathrm{WWC}}^{\mathrm{peri}}/F_1)\psi (\nu )]e^{-\tau (\nu )}\},\end{array}$$(6)

where $F_{\mathrm{WWC}}^{\mathrm{peri}}$ is the excess emission of the WWC cone during periastron passage, while ψ(ν) stands for the function that modulates this excess emission over the orbital cycle and needs to be defined. Our model assumes isotropic and completely ionized wind of the primary component with a standard exponent of its β-velocity profile, β = 1. According to these assumptions and theoretical formulae derived by Lamontagne et al. (1996, their Sect. 3), modified later by Shenar et al. (2021, their Sect. 3.5) to account for elliptical orbits, it is possible to calculate τ(ν) analytically from the following formula:

$$\begin{array}{cc}& \tau =\frac{(1+X_{\mathrm{H}}){\sigma }_{\mathrm{e}}{\dot{M}}_{1,\mathrm{wind}}}{8\pi m_{\mathrm{p}}{v}_{1,\infty }r}\frac{2}{\sqrt{(1-{\epsilon }^2)(1-b^2)}}[arctan\sqrt{\frac{1+b}{1-b}}\hfill \\ \hfill & +arctan(\sqrt{\frac{1+b}{1-b}}tan\frac{arcsin\epsilon }{2})],\hfill \end{array}$$(7)

where X_H is the mass fraction of hydrogen in the wind of the primary component, σ_e is Thompson electron scattering cross-section, m_p is mass of the proton, v_1, ∞ denotes the terminal velocity of the primary’s wind, and r is the instantaneous distance between components. Parameters ε and b are defined as follows

$$\begin{array}{c}\hfill \epsilon =sinicos(\nu +\omega +\pi /2),\end{array}$$(8)

and

$$\begin{array}{c}\hfill b=\frac{R_1/r}{\sqrt{1-{\epsilon }^2}}.\end{array}$$(9)

Distance r between the components can be easily obtained according to the well-known formula for elliptical orbits,

$$\begin{array}{c}\hfill r=\frac{a(1-e^2)}{1+ecos\nu },\end{array}$$(10)

where the semi-major axis of the relative orbit, a, has to be calculated from Kepler’s third law. The last necessary ingredient, that is ψ(ν), can be approximated in the following way. As mentioned above, the WWC cone shines primarily due to the scattering of the primary component’s light on free electrons, whose flux varies approximately as r⁻². Furthermore, the density of WWC cone decreases approximately as r⁻¹ (e.g., Cantó et al. 1996; Ignace et al. 2022). Therefore, we can expect that these two processes (i.e., the decrease in gas density in the WWC cone and the decrease in the number of photons reaching it from the primary component with increasing r) will largely determine the changes in the brightness of the WWC cone. Hence, we can reasonably assume that (after Shenar et al. 2021)

$$\begin{array}{c}\hfill \psi (\nu )\approx {\left(\frac{r_{\mathrm{peri}}}{r}\right)}^{{\gamma }_{\mathrm{WWC}}}={\left(\frac{1+ecos\nu }{1+e}\right)}^{{\gamma }_{\mathrm{WWC}}},\end{array}$$(11)

where r_peri is the periastron distance between the components and we expect an exponent γ_WWC ≳ 3, for the reasons described above.

The reason why we cannot directly use the formulae for Δm from Lamontagne et al. (1996) is that these authors assumed optically thin wind (i.e., τ ≪ 1), which probably is not the case for the intense wind of ExtEV. Moreover, our Eq. (6) intentionally does not contain the term associated with the atmospheric eclipse of the primary component by the wind emanating from the secondary. In our preliminary analysis, we noticed that even after including this additional atmospheric eclipse, the impact of the secondary’s wind on the light curve of ExtEV is negligible.

5.2. Optimized solution

To fit the analytical variability model described in the previous section, we phased the TESS light curve (Sect. 2.1) with the orbital period and binned it in 0.002 intervals to preserve the original shape of the heartbeat-like feature as accurately as possible (Fig. 13). Due to the presence of TEOs with different amplitudes (KS22), part of their signal remains visible in the binned light curve, while the stochastic variability has been effectively averaged out.

Fig. 13. Plot summarizing the fit of our analytical variability model (Sects. 5.1 and 5.2) to the TESS light curve of ExtEV. The upper panel shows the TESS light curve binned in the orbital phase (black dots) with the best-fit model (blue curve) superimposed. The pink line corresponds to the variability model generated with PHOEBE. Phases A, B, and C from Fig. 12 are shown with vertical dashed lines and labeled. The vertical shaded region marks the range of orbital phases where we expect the secondary component to be partially obscured by the WWC cone. The lower panel shows the residuals from the best fit and their smoothed version (red curve). The zero phase corresponds to the time of periastron passage. More details can be found in the main text.

For practical reasons, we fit the Δm′ function to such binned light curve, which slightly differed from the form indicated in Eq. (6). The fit function was transformed from Δm to the following form

$$\begin{array}{c}\hfill \mathrm{\Delta }{m}^{\prime }=\mathrm{\Delta }m({\phi }_{\mathrm{orb}}+\mathrm{\Delta }{\phi }_{\mathrm{orb}})-\mathrm{median}\{\mathrm{\Delta }m({\phi }_{\mathrm{orb}}+\mathrm{\Delta }{\phi }_{\mathrm{orb}})\}+C,\end{array}$$(12)

where Δφ_orb denotes a small phase shift between the TESS data and the model due to the finite precision of our T_peri, and C is a constant that accounts for any subtle vertical shift. Since our model does not account for the “eclipse” of the secondary component by the WWC cone near the inferior conjunction¹¹, we excluded the phase range between 0.1 and 0.3 from the fitting procedure. However, we would like to emphasize that, after several tests, this did not significantly affect the fitting results.

Before starting the fitting process, the calculation of τ (Eq. (7)) required two additional steps. First, we had to select a specific value of X_H. In our case, we have set the canonical value X_H = 0.7. Secondly, Ṁ_1,wind and v_1, ∞ are highly correlated, so it is not possible to fit these two parameters independently. We can only optimize their ratio, which corresponds to the mass of the wind contained in a shell of thickness dr at a distance from the primary component where the wind reaches its terminal velocity. We denote this ratio as (dM/dr)_∞. We performed fitting and estimated its errors using a method analogous to that described in Sect. 3. This time, however, we used a different form of the nuisance parameter. We assumed that the error associated with each point in the binned light curve is equal to ξ. In this way, we have incorporated in the fitting the scatter due to TEOs, which was not originally included in our model. The free parameters for the optimization were therefore i, log(dM/dr)_∞, F₂/F₁, $F_{\mathrm{WWC}}^{\mathrm{peri}}/F_1$, γ_WWC, Δφ_orb, C, and log ξ. We emphasize that parameters such as e, ω and f(M), taken from Table 6, were fixed during fitting, strongly constraining possible solutions. We also recall that M₁ was set to 35 M_⊙, but discuss the impact of changing this value in Appendix A. The results of our fitting are shown in Fig. 13, while the optimized parameters with their errors are presented in Table 8. We include a corner plot from the corresponding MCMC analysis of the solution in Appendix B.

5.3. Interpretation of the solution

It is surprising that the simple analytical model we presented above fits the TESS light curve (Fig. 13) so well. This is evidenced by the residuals, which, due to the partially remaining TEOs, oscillate around zero, except in the vicinity of the inferior conjunction, which we write about in more detail below. Although the model has five free parameters that determine its shape¹², we want to emphasize that the fit has been obtained without modifying the orbital parameters. With the exception of the orbital inclination, the others were fixed during optimization. To determine to what extent phenomena such as ellipsoidal distortion and irradiation effect may contribute to the ExtEV’s light curve, we generated a synthetic TESS light curve using PHOEBE 2.4 software (pink curve in Fig. 13). For this purpose, we used the following set of input parameters. The P_orb, e and ω were taken from Table 6, while i was taken from Table 8. For the primary component we adopted M₁ = 35 M_⊙ and R₁ from Table 4, while for the secondary component we took M₂ = 13 M_⊙ and R₂ = 5 R_⊙. The effective temperatures of the primary and secondary components were taken from J21. The atmospheres of both components were treated as black bodies. The gravity-darkening coefficients, as well as the bolometric albedos, were set to 1 for both components. The limb-darkening coefficients were assumed, according to the tables from Claret (2017) for the TESS passband. We assumed that both components rotate pseudo-synchronously. As can easily be seen in Fig. 13, proximity effects make a secondary contribution to the high amplitude of the light variability of ExtEV.

Although we originally fit a value of log(dM/dr)_∞, one can easily estimate the resulting value of Ṁ_1,wind from this parameter, realistically assuming v_1, ∞ ≈ 1000 km s⁻¹ (e.g., Kudritzki & Puls 2000). We obtained that the primary component loses mass due to the stellar wind at a rate of ∼4.5 × 10⁻⁵ M_⊙ yr⁻¹. Such a high mass loss rate makes the primary component of ExtEV unique and explains many of its observed properties. We provide further discussion on this topic in Sect. 6.1.

Since the atmospheric eclipse is very sensitive to the inclination of the orbit, this parameter is determined with relatively high precision (Fig. B.1). Therefore, we can expect M₂ ≈ 13 M_⊙. We have verified that with this orbital configuration and the radius of the secondary component R₂ ≈ 5 R_⊙, the system is on the borderline between a pure atmospheric eclipse and an atmospheric eclipse combined with a so-called grazing photospheric eclipse. The obtained F₂/F₁ ratio of about 0.15 theoretically indicates that the spectra of the components can be disentangled. However, due to, for example, fast rotation of the secondary component, its lines may still be strongly diluted in the composite spectrum.

The residuals (Fig. 13, phase C) reveal a decrease in brightness of about 0.02 mag near inferior conjunction, which is not accounted for in our model. We interpret this decrease as attenuation of light from the secondary component by turbulent structures in the WWC cone that are present in this range of orbital phases between the observer and the secondary component.

5.4. Scatter of the residuals over the orbital cycle

With the fit variability model for ExtEV, we can subtract it from the original TESS light curve and examine the behavior of the residuals by calculating the local scatter as the standard deviation of the residuals over a specified range of orbital phases. Figure 14 shows the result of this procedure. It is clear that the scatter of the residuals (Fig. 14, bottom panel) is not random but depends on the orbital phase. This scatter is due to at least three phenomena: the TEOs (as they change their amplitude with time), the stochastic variability of the primary component (caused, for example, by internal gravity waves, Rogers et al. 2013), and the turbulent gas structures in the WWC cone. In Fig. 14 we can see that the scatter of the residuals is the largest during the periastron passage. We interpret this increase as reflecting increased turbulence of the gas in the WWC cone, causing the light curve to “flicker” with greater intensity. This seems to be a plausible explanation for several reasons. Firstly, it is at the periastron passage when the secondary component moves through the densest medium at the highest relative velocity. Secondly, the WWC cone has the highest density at this particular orbital phase, so it effectively scatters light. Finally, the WWC cone in the periastron has the smallest geometric thickness and can be subject to the so-called thin-shell instability, making the gas motion highly turbulent (e.g., Kashi & Michaelis 2021).

Fig. 14. Changes of the scatter of residuals from the fit of our model to the light curve of ExtEV. The upper panel shows the TESS light curve (black points) folded with the orbital period, while the red curve corresponds to the best-fitting model from Fig. 13. Phase φorb = 0 corresponds to the periastron passage. The middle panel presents the residuals from the fit shown in the upper panel. The lower panel shows the changes in the local scatter of the residuals, calculated as the standard deviation of the residuals in 0.007 phase intervals. The red arrows indicate a sudden decrease in the scatter of residuals around the phase of superior conjunction.

A very interesting feature of the scatter of the residuals in Fig. 14, indicated by the red arrows in the lower panel, is the clear and rapid decrease in the scatter around the phase of superior conjunction (φ_orb ≈ −0.04). This decrease can be easily interpreted based on the model we have adopted. When the secondary component with the WWC region is behind the primary component, there is a partial occultation of the WWC cone by the atmosphere, and possibly also by the photosphere, of the primary component. Since it is mainly the WWC cone that is responsible for the changes in the scatter of the residuals, the increasing contribution of the turbulent WWC emission in this phase to the total brightness of the ExtEV system suddenly decreases, leading to a significant reduction of scatter in the residuals.

6. Discussion

6.1. Wind mass-loss rate

The mass-loss rate we derived for ExtEV, 4.5 × 10⁻⁵ M_⊙ yr⁻¹, is exceptionally high compared to a typical blue SG with parameters similar to the primary component of ExtEV. To see how this compares with theoretical predictions, we calculated the values of mass-loss rate based on three selected wind mass-loss prescriptions and the physical parameters of the primary component that we determined in our study. The result of this comparison is presented in Table 9. It is clear that the wind generated by the primary component can be up to two orders of magnitude more intense than what contemporary theory predicts for a single star of this type. If our proposed scenario of light variability is correct, this implies that the primary component of ExtEV exhibits a wind mass loss rate typical of WR stars or red SGs from the tip of the asymptotic giant branch¹³ (AGB), although it is still relatively close to the MS and its spectra definitely do not exhibit spectral features of WR stars.

A significant mass-loss rate found in the ExtEV should result in noticeable changes in the orbital period. Indeed, KS22 have shown that the orbital period of the ExtEV is decreasing at a rate of approximately 11 s yr⁻¹. They also estimated that if we attribute the observed changes in orbital period solely to mass loss due to the stellar wind of the primary component¹⁴, the rate should be on the order of 10⁻⁴ M_⊙ yr⁻¹. Given that tides and TEOs also contribute to the shortening of the orbital period, the aforementioned value of Ṁ_1,wind fits very well with the proposed high mass loss in the ExtEV. We conclude that ExtEV is a WWC binary system with a deep atmospheric eclipse and strong excess emission during periastron passage.

In this context, an intriguing question arises regarding the source of such a high mass loss by the primary component. We can suspect that this is due to the presence of a close and massive companion. Because of the significant eccentricity of its orbit, along with the change in orbital phase, the primary’s Roche lobe periodically approaches the surface of the star, which can significantly enhance its stellar wind (e.g., Tout & Eggleton 1988; Lei et al. 2014; MacLeod & Loeb 2020). Moreover, the primary component shows the presence of relatively high-amplitude TEOs, which may also be responsible for the enhancement of the stellar wind. Thus, it appears that ExtEV is an extremely rare massive binary system that offers us a unique opportunity to observe a short phase of its evolution, during which eccentricity, proximity to a massive companion, rotation and the presence of TEOs significantly enhance the stellar wind, and their mutual interaction determines the evolution of the whole system.

6.2. Rotation rate of the primary component

Knowing the inclination angle of the orbit of ExtEV and assuming that the rotational axis of the primary component is perpendicular to the orbital plane, we can derive linear equatorial velocity v_eq ≈ 190 ± 38 km s⁻¹, bearing in mind that the observed v_eq sin i = 174 ± 34 km s⁻¹ (J21). This is rather an upper limit, since J21 measured v_eq sin i from a single line of He Iλ 4922 Å, neglecting its distortion and broadening caused by several TEOs. The value of v_eq implies that the primary component is likely to rotate with a period of about 8.0 ± 2.0 d, provided that it is not significantly oblate due to rotation. This is twice as long as the rotation period adopted by ML23 to obtain the large amplitude of light variations in their hydrodynamic model.

Knowing v_eq and the other parameters of ExtEV obtained in our paper, we can calculate the critical rotation period of the primary component following Gagnier et al. (2019)

$$\begin{array}{c}\hfill P_{\mathrm{crit}}=2\pi {\left[\frac{G{M}_1}{R_{\mathrm{eq}}^3}(1-{\mathrm{\Gamma }}^{3/2})\right]}^{-1/2},\end{array}$$(13)

where R_eq is the equatorial radius of the primary component if it would rotate critically. We can realistically estimate R_eq ≈ 1.5R₁ (e.g., Kopal 1959), assuming the Roche model of a rotating star and that its polar radius, R_p, does not change with rotation rate, that is R_p ≈ R₁. In turn, Γ = κL₁/(4πcGM₁) stands for the so-called Eddington factor. In the definition of Γ, κ is the flux-weighted mass absorption coefficient in the photosphere of the primary component. Since the temperature of the sub-photospheric layers increases rapidly, we can assume that κ is dominated by electron scattering in the outer layers, that is κ ≈ κ_e ≈ 0.34 cm² g⁻¹. Thus, Γ ≈ 0.5 and we are left with P_crit ≈ 9.2 ± 1.2 d¹⁵. Theoretically, this value would indicate a critical rotation of the primary component. However, the problem is that if the primary component rotates (nearly) critically, our transformation of v_eq to P_rot using the relation P_rot = 2πR_eq/v_eq, where we assumed R_eq = R₁, is no longer valid. To examine the impact of nearly critical rotation on the P_rot and P_crit presented above, we performed the following estimations. The value of R₁ obtained in Sect. 3 is a proxy for L₁, since we considered the effective temperature of the primary component to be known. The luminosity of a rotating primary component can be approximated as $L_1\approx 4\pi R_{\text{p}}^2{\sigma }_{\text{SB}}{T}_{\text{p}}^4\eta ({\mathrm{\Omega }}_{\text{rot}}/{\mathrm{\Omega }}_{\text{crit}})$, where η(1)≈0.64¹⁶ is the dimensionless factor accounting for the latitude-dependent surface temperature distribution, σ_SB is the Stefan-Boltzmann constant, and T_p stands for surface temperature at the poles. If the rotation of the primary component approaches critical rotation, the R₁ obtained by us in terms of SED fitting changes its original meaning to $R_1\sim \sqrt{\eta (1)}{R}_{\mathrm{p}}\approx \sqrt{0.64}{R}_{\mathrm{eq}}/1.5$, where we again applied the relation R_eq = 1.5R_p. Finally, it leads us to R_eq ≈ 1.87R₁. This result translates into P_rot = 14.9 ± 3.7 d and P_crit = 12.7 ± 1.7 d, which points toward subcritical rotation, contrary to the previous result. We therefore conclude that without detailed modeling, it is difficult to unambiguously answer the question about the value of Ω_rot/Ω_crit.

We also checked how the rotation period of the primary component compares to the pseudo-synchronous rotation period, P_ps, that can be expected under effective tidal action. Taking the value of e from Table 6 and using the formulae derived by Hut (1981), we obtained P_ps ≈ 8.45 d. This value is satisfactorily consistent with the rotation period of the primary component we obtained under the assumption of a subcritical rotation rate.

6.3. Remarks on X-ray signal and optical emission features

Wind-wind collision binaries are also known to emit significant X-ray flux, which can vary with the orbital phase (see Rauw & Nazé 2016, for a comprehensive review on this topic). The primary source of X-ray radiation in these systems is the front of the WWC cone, the temperature of which (depending on the properties of the system) may reach approximately 10⁷ K (e.g., Parkin & Gosset 2011). Therefore, one would expect that since ExtEV is such an extreme case of a WWC binary, it would have a high X-ray luminosity, L_X, which should be the strongest at the periastron passage. Unfortunately, J21, after carefully examining the Swift X-ray Telescope data (their Sect. 3.6), did not find any statistically significant X-ray emission at the position of ExtEV. However, the authors derived an upper limit to the X-ray brightness, L_X ≲ 2.5 × 10³⁴ erg s⁻¹. Comparing this value to other WWC binaries (Cherepashchuk 2000, and references therein), it turns out that it is difficult to detect even a strong X-ray signal with Swift due to the distance separating us from the LMC. ExtEV may therefore be characterized by significant X-ray emission, but an X-ray observatory with better capabilities than Swift is needed to detect it. Moreover, the extreme conditions in the WWC cone do not automatically result in the strong X-ray flux observed outside the system. This can happen for at least three reasons. First, dispersed gas and dust around ExtEV, but also the gas within the WWC cone, can effectively absorb X-ray photons generated by the collision of winds as is in the case of HD 38282 (Shenar et al. 2021). Secondly, if the front of the WWC cone falls near the periastron (or even throughout the orbital phase) onto the surface of the secondary component, there is a sharp drop in L_X due to the lack of the hottest and densest region of the WWC, as is the case with η Car, which is also a WWC binary with a very high orbital eccentricity (Teodoro et al. 2012). Finally, as recently shown by Mackey et al. (2023), inverse Compton scattering occurring in the WWC shocked region, often neglected during the modeling of WWC binaries, can effectively suppress the X-ray emission. This is because the mentioned physical process acts as an additional cooling mechanism that softens the spectrum emerging from the densest parts of the WWC cone.

Although the optical spectra of the ExtEV exhibit strong emissions in the Balmer Hα and Hβ lines, we do not observe their counterparts in the helium or ionized metal lines, which are typical tracers of recombination regions in the WWC cone. However, in some He I lines we do observe clear emission features located in the wings. As with X-rays, this does not necessarily imply that ExtEV is not a source of strong emissions in the aforementioned lines. As noted by Koenigsberger et al. (2022) for the HD 5980 WWC system, it is very difficult to determine whether the moderate wing emission is the result of two separate Doppler-shifted emission components or whether it is a strong and broad emission line, which is superimposed with almost equally strong absorption with a slightly narrower profile. Most likely, photons emitted in lines located in the WWC recombination zone encounter the dense, optically thick wind of the primary component, potentially leading to their effective absorption. For the reasons outlined here, the analysis of emission features in the ExtEV optical spectrum cannot be unambiguous without an advanced hydrodynamic model of the entire system, taking into account its dynamics, the structure of the WWC cone, stellar wind and radiative transfer. This issue is obviously beyond the scope of our paper and requires further examination using sophisticated, time-consuming simulations.

7. Summary and conclusions

We analyzed photometric properties, light curves from TESS and LCOGT, and SALT/HRS optical spectra for the massive eccentric binary system, MACHO 80.7443.1718 (ExtEV), located in the LMC. This system shows extremely large brightness variations (Fig. 1), which raises numerous questions about their nature (Sect. 1). Our main findings can be summarized as follows.

1. Thanks to the LCGOT light curves of ExtEV collected in the U, B, and V bands (Sect. 2.1, Fig. 2), we conclude that the range of brightness changes and the shape of the light curve do not depend on wavelength over a broad range from NUV to NIR.

2. SED of ExtEV reveals significant IR excess (Sect. 3), especially prominent in the MIR bands (Fig. 3). Although the IR properties of ExtEV indicate an enhanced mass loss rate in the system, they do not match the properties of B[e] SGs (Sect. 4.1, Fig. 6). Moreover, in our SALT/HRS spectra we did not detect any emission lines, such as [O I] lines, typical for the B[e] phenomenon (Fig. 5).

3. SED modeling resulted in the radius of the primary component R₁ ≈ 30 R_⊙, with a luminosity of about log(L₁/L_⊙)≈5.82 (Table 4). The SED indicates also for the presence of an additional attenuation around 2760 Å of unknown origin (Fig. 4).

4. Determination of the initial and current mass of the primary component is difficult, primarily due to the strong influence of rotation and mass loss on its evolution (Sect. 4.2). However, our simulations showed that the mass of the primary component at ZAMS most likely falls within the range of 27 to 55 M_⊙. The present mass of the star lies between 25 and 45 M_⊙ (Fig. 7). The evolutionary status of the primary component, as deduced from the comparison with evolutionary tracks, corresponds to either hydrogen or helium core burning.

5. We derived spectroscopic parameters of the primary’s orbit combining archival RVs of J21 and RVs obtained from spectra collected with SALT/HRS (Sect. 4.3, Fig. 8, Table 6). In particular, we derived a more precise value of the mass function f(M)≈0.74 ± 0.05 M_⊙. We also found significant changes in the γ velocity of ExtEV, an indication that the system is triple.

6. We demonstrated that the assumption that the emission features in the Hα and Hβ lines originate in a Keplerian disk around the primary component leads to the location of the disk that would strongly interact with the secondary’s orbit (Sect. 4.5, Fig. 11). We also argued that the primary component likely retains a detached geometry during the periastron passage (Sect. 4.4, Fig. 9), and its rotation period of about 8 d is consistent with pseudo-synchronous rotation (P_ps ≈ 8.45 d; Sect. 6.2).

7. The light curve of ExtEV can be successfully explained by a superposition of an atmospheric eclipse of the secondary component by the intense stellar wind of the primary and excess emission originating from scattered light on relatively dense structures of the WWC cone (Sect. 5, Fig. 13). As a result of our modeling, we precisely determined the inclination of the ExtEV’s orbit (i ≈ 66°) and the wind mass-loss rate from the primary component, log(Ṁ_1,wind/(M_⊙ yr⁻¹)) = −4.35 (Table 8).

Based on these results, we can draw the following conclusions.

1. ExtEV is not a B[e] SG because its several key properties are significantly different from those of B[e] SGs. This also means that the extreme changes in the brightness of ExtEV cannot be explained by a disk surrounding the primary component.

2. It seems unlikely that the light variability in ExtEV is caused by an extreme tidal distortion and subsequent nonlinear breaking of tidal waves on the surface of the primary component, as proposed by ML23. First, the system most likely does not undergo RLOF at periastron, and the rotation rate of the primary component is at least twice lower than that assumed in the ML23 models. Second, the shape and peak-to-peak amplitude of the light curve of ExtEV are virtually the same in Johnson UBV and TESS passbands, which rules out the presence of extreme tidal distortions of the primary component that would be accompanied by a significant gravitational darkening and hence heterochromatic effects in the light curve.

3. Ellipsoidal distortion appears to be a secondary factor in shaping the light curve of ExtEV, and therefore the star is not an extreme case of an EEV system. In other words, it is not an extreme heartbeat star. According to our study, ExtEV may be in fact a unique WWC binary system in which a considerable mass loss by the primary component causes a remarkably large range of light variability close to periastron passage. Therefore, a particular care must be taken when classifying massive binary systems as EEVs based on their light curves, as some of them may actually be WWC systems mimicking EEVs.

4. The rate of mass loss in the primary component that we obtained is two orders of magnitude greater than what the theory predicts. We can therefore suspect that the wind of the blue SG in ExtEV is tidally and rotationally enhanced, possibly with some additional interaction with TEOs. In this scenario, ExtEV is an extremely rare observational example of a massive binary system during a short but dramatic stage in its evolution, when the proximity of a massive companion leads to a sudden stripping of the primary’s envelope without the need for RLOF. If this scenario holds in ExtEV, it is likely that other massive binary systems also experience such episodes of significant mass loss, even on the MS, without having to become WR stars or red giants. Thus, ExtEV can serve as an excellent laboratory for studying the mechanisms of wind enhancement in massive stars and predicting the impact of this enhancement on their evolution.

In order to verify the variability scenario proposed in our study, we are going to perform IR interferometric observations of ExtEV and collect its UV spectra, which will help us to determine the parameters of its stellar wind, including an independent assessment of its mass-loss rate. In a final solution of the puzzling nature of these stars, it would also be beneficial to perform polarimetric measurements of this system in different orbital phases.

Acknowledgments

We would like to thank the anonymous referees for many useful comments and suggestions that helped to improve the publication. PKS would like to express his gratitude to Milena Ratajczak and Marcin Wrona for their assistance in preparing observing time proposals for ExtEV. PKS is also grateful to Morgan MacLeod for the fruitful discussion about the nature of ExtEV, and to Damien Gagnier for the discussion on some properties of rotating stars. Some of the observations reported in this paper were obtained with the Southern African Large Telescope (SALT) under program 2021-1-SCI-022 (PI: PKS). Polish participation in SALT is funded by grant No. MEiN nr 2021/WK/01. This work makes use of observations from the Las Cumbres Observatory global telescope network. PKS was supported by the Polish National Science Center grant no. 2019/35/N/ST9/03805. AP, PKS, and PŁ would like to appreciate the financial support from the Polish National Science Center grant no. 2022/45/B/ST9/03862. TR was partly founded from budgetary funds for science in 2018-2022 in a research project under the program “Diamentowy Grant”, no. DI2018 024648. This publication makes use of VOSA, developed under the Spanish Virtual Observatory (https://svo.cab.inta-csic.es) project funded by MCIN/AEI/10.13039/501100011033/ through grant PID2020-112949GB-I00. VOSA has been partially updated by using funding from the European Union’s Horizon 2020 Research and Innovation Programme, under Grant Agreement no. 776403 (EXOPLANETS-A). The authors made use of the Strasbourg Astronomical Data Center (CDS) portal and the Barbara A. Mikulski Archive for Space Telescopes (MAST) portal. This paper includes data collected by the TESS mission, which are publicly available from the MAST. This research has made use of “Aladin sky atlas” developed at CDS, Strasbourg Observatory, France. This research has made use of the VizieR catalog access tool, CDS, Strasbourg, France. This research made use of NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Matplotlib (Hunter 2007) and AstroPy (Astropy Collaboration 2022).

References

Appendix A: The impact of change in M₁

Since the present mass of the primary component, M₁, is not known with high precision (Sect. 4.2), we investigate how the change in M₁ affects the results presented in Sect. 5.2, especially the parameters presented in Table 8 obtained assuming M₁ = 35 M_⊙. For this purpose, we considered two other cases, M₁ = 30 M_⊙ and 40 M_⊙. The fitting procedure remained unchanged. The effect of this experiment is summarized in Table A.1. All parameters but inclination change only slightly with the change of M₁. The change of i between the two solutions is about 1.9° and is at the 5σ level of the uncertainty of i. Therefore, we can safely state that the conclusions drawn in our study are virtually independent of the assumption of the value of M₁ in the realistic range of this parameter.

Appendix B: Corner plot

To generate the corner plot presented in Fig. B.1 we used 32 walkers, each with 100 000 steps, from which 500 were considered as “burn-in” steps. Before generating the corner plot, each of the chains was “thinned” by half of the autocorrelation time, which amounted to about 100 steps. The solution does not reveal any significant correlations between the free parameters of the model, except for the obvious one between F₂/F₁ and i.

Fig. B.1. Corner plot resulting from our MCMC simulations described in Sect. 5.2. A series of vertical dashed lines in each histogram marks the position of 0.18, 0.5, and 0.84 quantiles.

All Tables

All Figures

	Fig. 1. TESS light curve of ExtEV phased with the orbital period of 32.83016 d. Orbital phase φorb = 0.0 corresponds to the time of the periastron passage derived in Sect. 4.3. The inset shows the light curve near the periastron passage. Stochastic light variations and TEOs can be seen over the entire range of orbital phases.
In the text

Fig. 2. LCOGT light curves of ExtEV in Johnson U (purple dots), B (blue dots), and V (green dots) passbands. Panel a shows the light curves. A ±0.25 mag vertical offset was applied to U and V-filter light curves to avoid overlapping. Panels b–d show the light curves phased with the orbital period. Phase zero corresponds to the time of periastron passage. In addition, we have drawn dash-dotted lines in each of the lower panels to indicate the range of variability during the heartbeat. The plot illustrates that the range is virtually the same in all three passbands.

In the text

Fig. 3. Results of the SED fitting to the multiband photometry of ExtEV. The upper panel shows the absolute fluxes calculated from magnitudes given in Table 3 (dots are color-coded according to the legend). The model spectrum of the primary component is shown with a black line, its reddened versions with a red line (with additional attenuation) and a blue line (without it). The infrared excess is clearly visible and indicated by the shaded area. The middle and lower panels show the residuals from the best fit, with and without the additional attenuation, respectively. For residuals presented on the left (UV and optical range), the ordinate axes are linearly scaled while the residuals in the IR range are normalized by the model values and logarithmically scaled.

In the text

Fig. 4. Extinction law resulting from the best fit to the observed SED of ExtEV. The black curve refers to the interstellar dust extinction function of Fitzpatrick (1999) with the parameters given in the legend. The additional extinction feature is shown with the gray line, while the red dashed curve corresponds to the total extinction function (the product of black and gray curves). With a pair of blue vertical dotted lines, we have marked the positions of the two local minima in the total extinction function. More details can be found in the main text.

In the text

Fig. 5. Small sections of the SALT spectra showing [O I] emission lines at their laboratory wavelengths but not in the rest frame of ExtEV. The rest-frame wavelengths of the lines shown are indicated in each panel in the top left corner. The spectra are shifted vertically according to the orbital phase. The vertical red and blue dashed lines correspond to RVs of 0 (geocoronal emission) and 313 km s−1 (systemic velocity of ExtEV), respectively. For comparison, the bottom right panel shows the He I absorption line originating from the photosphere of the primary component of ExtEV.

In the text

Fig. 6. Infrared color-color diagrams for ExtEV. The position of ExtEV is marked by the cyan star in each panel. The blue dots and red triangles correspond to the B[e] SGs and LBVs located in the Magellanic Clouds, respectively. Their photometric properties were adopted from Kraus (2019). Black crosses refer to the sample of late-type SGs selected by González-Fernández et al. (2015). The shaded region in panel (c) shows the area with color index (Ks − W3)≳2 mag, where SGs with enhanced mass loss are suspected to be located.

In the text

Fig. 7. HR diagrams summarizing our simulations described in Sect. 4.2. Black and blue dots show the location of the primary component of ExtEV, according to this paper and the study of J21, respectively. The gray curves show evolutionary tracks that never approach the position of the primary component determined in our study. Conversely, multicolored tracks correspond to those that are within the error box of the ExtEV (black dot) at some phase of the evolution. The color bar above the HR diagrams reflects M1, ZAMS (left panel) and M1 (right panel) and is common to both panels.

In the text

Fig. 8. Optimized spectroscopic solution for the RVs of the primary component of ExtEV. The points on all panels are marked with different colors depending on the spectrograph (labeled in the top panel). The blue curve represents the best-fit model. (a) RV curve including orbital solution and the long-term changes. (b) RV curve with the orbital solution subtracted. (c) RVs phased with the orbital period given in Table 6 after freeing from the long-term changes of RVs shown in panel b.

In the text

Fig. 9. Histogram of the distribution of parameter ϑ = (RL − R1)/R1, resulting from the simulations described in Sect. 4.4. The vertical dashed black line indicates ϑ = 0, when the system becomes semi-detached. The vertical blue lines indicate quantiles at levels of 0.16, 0.5, and 0.84. Only about 3% of samples correspond to an RLOF at periastron, while their majority (∼97%) points toward the detached geometry at periastron.

In the text

	Fig. 10. Fragments of SALT/HRS spectra of ExtEV near the Hα (upper panel) and Hβ (lower panel) lines. The spectra are color-coded according to the orbital phase (the color bar is located at the top of the figure), with zero phase corresponding to the periastron passage. The narrow emission in the center of the spectrum is of nebular origin.
In the text

Fig. 11. Models of the Keplerian rings around the primary component of ExtEV. Each set of figures (a, b, and c) corresponds to a different orbital phase φorb provided above each set. The right-hand panels show the observed Hα (salmon) and Hβ (aqua) lines to which the theoretical profiles (black curves) were fit. For each spectrum, two rings were fit separately, one responsible only for the Hα emission and the other for the Hβ emission. The narrow nebular emission of the LMC (gray region bounded by vertical dashed lines) was excluded from the fit. The left-hand panels show the resulting ring geometry compared to the size of the primary component (black circle) and the size of the orbit of the secondary component (dashed ellipse, secondary component size is not to scale). More details and discussion can be found in the main text.

In the text

Fig. 12. Artistic representation of the orbit (silver ellipse) of the secondary component of ExtEV (dark blue circle) around the primary component (light blue circle near the center). The size of the orbit and both components are to scale. The orientation of the orbit corresponds to the view of the ExtEV in the sky (i.e., as seen by the observer on Earth). The primary component is the source of a slow but dense stellar wind (semi-transparent spheres, concentric with the primary component), which, after colliding with the surface of the wind of the secondary component, forms a WWC cone, a turbulent structure, with its apex located near the secondary component. Three orbital phases are labeled, corresponding to: the superior conjunction (A), periastron passage (B), and inferior conjunction (C). More details can be found in the main text.

In the text

Fig. 13. Plot summarizing the fit of our analytical variability model (Sects. 5.1 and 5.2) to the TESS light curve of ExtEV. The upper panel shows the TESS light curve binned in the orbital phase (black dots) with the best-fit model (blue curve) superimposed. The pink line corresponds to the variability model generated with PHOEBE. Phases A, B, and C from Fig. 12 are shown with vertical dashed lines and labeled. The vertical shaded region marks the range of orbital phases where we expect the secondary component to be partially obscured by the WWC cone. The lower panel shows the residuals from the best fit and their smoothed version (red curve). The zero phase corresponds to the time of periastron passage. More details can be found in the main text.

In the text

Fig. 14. Changes of the scatter of residuals from the fit of our model to the light curve of ExtEV. The upper panel shows the TESS light curve (black points) folded with the orbital period, while the red curve corresponds to the best-fitting model from Fig. 13. Phase φorb = 0 corresponds to the periastron passage. The middle panel presents the residuals from the fit shown in the upper panel. The lower panel shows the changes in the local scatter of the residuals, calculated as the standard deviation of the residuals in 0.007 phase intervals. The red arrows indicate a sudden decrease in the scatter of residuals around the phase of superior conjunction.

In the text

	Fig. B.1. Corner plot resulting from our MCMC simulations described in Sect. 5.2. A series of vertical dashed lines in each histogram marks the position of 0.18, 0.5, and 0.84 quantiles.
In the text