© The Authors 2023

1. Introduction

Massive stars (M_initial > 8 M_⊙) in all evolutionary phases strongly affect their galactic neighborhood via stellar winds and ionizing fluxes. During core-H burning on the main sequence, they have spectral types O and early B. As massive stars evolve, their outer hydrogen (H) rich envelopes could be removed by stellar winds. The majority of massive stars are born in binary systems (e.g Sana et al. 2012, 2014; Moe & Di Stefano 2017), and interactions with a companion star may also lead to a removal of the outer envelope. Highly evolved massive stars that lost a large portion of their outer H-rich layers and have optically thick winds are spectroscopically classified as Wolf-Rayet (WR) stars. The WR stars are typically hotter and have stronger, optically thick winds when compared with their evolutionary predecessors. These stars come in two major subtypes, WN and WC, which are spectroscopically identified by strong emission lines of nitrogen and carbon, respectively. WR stars end their lives with a core-collapse event, likely leading to the formation of black holes (BHs; e.g., Sukhbold et al. 2016; Gal-Yam et al. 2022). Understanding the formation pathways of WR stars in nearby low-metallicity galaxies is needed for quantifying stellar feedback and compact object populations in conditions resembling the early Universe.

The Small Magellanic Cloud (SMC) galaxy has a metallicity of Z ≈ 1/7 Z_⊙ (Hunter et al. 2007; Trundle et al. 2007) and is nearby (d ≈ 61 kpc; Hilditch et al. 2005), thus providing an excellent test bed for investigating stars at low metallicity. The proximity of the SMC and its low foreground and intrinsic extinction allows for a detailed study of stellar and wind parameters of low-metallicity populations of O and WR stars. This enables us to obtain a realistic understanding of massive star evolution and feedback at low metallicity in general.

Yet what we have learned so far is perplexing. Only 12 WR stars exist in the SMC. All of them are very luminous and (except for one WO-type star) contain some hydrogen (Hainich et al. 2015; Shenar et al. 2016). According to the standard, single-star evolutionary tracks, their progenitors are O-type stars with masses ≳40 M_⊙ (Shenar et al. 2020a, and references therein). For an order of magnitude estimate, one can assume that a massive star spends 10% of its life in a WR stage. This implies that a galaxy containing approximately ten WR stars should contain approximately one-hundred massive O stars. However, there is a severe paucity of such O-stars in the SMC (Ramachandran et al. 2019; Schootemeijer et al. 2021).

Close binary evolution is expected to play a principal role in the formation of WR stars (Paczyński 1967; Kippenhahn & Weigert 1967). However, the importance of binarity is still under debate and different studies have come to different conclusions (e.g., Vanbeveren et al. 2007; Eldridge et al. 2017; Shenar et al. 2020a; Pauli et al. 2022a). Up to now, the evolution toward the WR stage in binaries has chiefly been studied for primary stars, while the evolution of secondary stars has been largely neglected. In this paper, we aim to gain insights into the evolution of both components in a binary with well-established stellar parameters.

One of the youngest and earliest O-type stars in the SMC is AzV 14, the main ionizing source of the H II region NGC 261 (see Fig. 1). Previously, single epoch optical spectra of AzV 14 have been analyzed by Massey et al. (2004) and Mokiem et al. (2006). The star was classified as O5 V with a reported spectroscopic mass of 74 M_⊙ − 90 M_⊙ being in the mass range of potential WR single-star progenitors. From a newly obtained light curve, we know now that AzV 14 is in fact a binary. In this paper, we complement recently obtained spectra of AzV 14 in the UV and optical by archival data, aiming at securely determining stellar parameters and, on this basis, to model the evolution of the AzV 14 components toward WR stages.

Fig. 1. False-color HST image composed of the images in the F475W (blue) and F657N (yellow) filters. The position of AzV 14 and the approximate size of the NGC 261 H II region are indicated by white circles.

2. Observations

Presently, six spectra of AzV 14 obtained at different epochs covering the far-UV, UV and optical wavelength ranges exist in telescope archives. One far-UV spectrum was taken with the Far Ultraviolet Spectroscopic Explorer (FUSE; Oegerle et al. 2000), two UV spectra with the Hubble Space Telescope’s (HST) Faint Object Spectrograph (FOS; Keyes et al. 1995) and Space Telescope Imaging Spectrograph (STIS; Branton & Riley 2021), and three optical spectra with the European Southern Observatory (ESO) Very Large Telescopes (VLT) Ultraviolet and Visual Echelle Spectrograph (UVES; Dekker et al. 2000) and the X-shooter (Vernet et al. 2011) spectrograph. One X-shooter spectrum was taken as part of the XShootU project¹ (from here on “X-shooter (2020)”) and one as part of our program, ID 109.22V0.001 (from here on “X-shooter (2022)”). A detailed description of the individual spectra and photometry can be found in Appendix A.

AzV 14 is close to the SMC bar, thus we adopt a distance modulus of DM = 18.9 mag (Westerlund 1997; Harries et al. 2003; Hilditch et al. 2005). The radial velocities (RV) of the different regions in the SMC are not uniform (e.g., De Propris et al. 2010). We estimate the RV shift of the NGC 261 complex by fitting Gaussians to interstellar absorption lines as v_NGC 261 = 148 ± 2 km s⁻¹. All spectra shown in this work are in the rest frame of NGC 261.

AzV 14 was observed by the NASA Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) mission (TIC 180206579) during its sectors 1, 27 and 28 in full-frame image (FFI) mode with a cadence of 30, 10 and 10 min, respectively. TESS has a low spatial resolution of 21″ px⁻¹. However, since AzV 14 is the brightest optical object in this region we can extract its light curve (see Appendix C). Furthermore, AzV 14 was detected by the Chandra X-ray telescope (Weisskopf et al. 2000).

3. Spectral analysis

3.1. Spectral line variability reveals AzV 14 as a close binary

The multi-epoch spectroscopy allows us to search for signs of binarity. In Fig. 2, we show selected He I and He II lines in the UVES and X-shooter spectra, taken at different epochs. Massey et al. (2004) and Mokiem et al. (2006) attributed the apparent emission feature in the middle of helium absorption lines of the UVES to a nebular contamination. New multi-epoch X-shooter spectra clearly show that the profiles of helium lines are complex and variable. They are well explained as composite lines originating from two components of a binary system. The UVES and X-shooter (2022) spectra must have been taken at a comparable orbital phase, as the line profile is similar in both spectra, while the X-shooter (2020) spectrum must have been obtained close to conjunction. This is confirmed by the binary ephemeris derived in Sect. 4.

Fig. 2. Comparison of selected lines in the UVES (solid blue) and the X-shooter (dashed pink and dotted yellow) optical spectra. The high resolution UVES spectrum is binned by 0.15 Å to match the X-shooter spectra. All spectra are convolved with a Gaussian having a FWHM = 0.05 Å for a better comparison. The lines in the UVES and X-shooter (2022) spectra are separated, while they are blended in the X-shooter (2020) spectrum. This provides a clear indication for binarity.

3.2. Establishing projected rotation velocity and RVs

The projected rotation velocity as well as the RV shifts of the individual binary components are necessary ingredients of spectral modeling. In the following we describe how we derived these quantities.

3.2.1. Projected rotation velocities

The helium lines in UVES and X-shooter (2022) spectra are resolved into two components with lines of similar depths and widths, indicating that the components of the binary system are similar stars. We employ the iacob broad tool (Simón-Díaz & Herrero 2014) to measure the projected rotational velocities for each binary component. For the fitting procedure the lines He IIλ4200, He Iλ4471, He IIλ4542, and He IIλ4686 were used, while cutting off the parts of the lines that are blended by the companion. Although helium lines are not only rotationally, but also pressure broadened, they were analyzed because of the lack of available metal lines. The iacob broad tool derives two different values for the projected rotational velocities, one from a Fourier transformation (FT) and one from the goodness of the fit (GoF). We quote the average value obtained by the different methods. As the estimates on the projected rotation velocities are based on blended lines, we adopt conservative error margins. We find that both binary components have similar projected rotational velocities of v_rot sin i = 90 ± 20 km s⁻¹.

3.2.2. RVs of the binary components in each spectrum

The RVs of all spectra are determined by shifting the synthetic spectra until they fit the observations using a Markov chain Monte Carlo method combined with a least square fitting (Pauli et al. 2022b, their Sect. 3.1.2 and their Appendix A).

The UVES spectrum shows RVs of v₁ = −101.9 ± 2.6 km s⁻¹ and v₂ = 80.5 ± 3.2 km s⁻¹ for the primary and the secondary, respectively. This is comparable to those of the X-shooter (2022) spectrum, yielding RV shifts of v₁ = −112.4 ± 2.6 km s⁻¹ and v₂ = 98.6 ± 2.9 km s⁻¹. In the X-shooter (2020) spectrum the lines are not split due to an orbital phase close to conjunction. We estimated v₁ = −46.1 ± 5.8 km s⁻¹ and v₂ = 22.3 ± 5.9 km s⁻¹ for the primary and secondary, respectively. Accurate measurements of RV shifts in the optical spectra allow us to estimate the temperature and surface gravity of the individual binary components (see Sect. 3.3.1).

In the STIS spectrum, the He IIλ1640 line is split into two individual components of similar strength (see Fig. D.4), yielding RVs of v₁ = −128.6 ± 2.5 km s⁻¹ for the primary and v₂ = 113.2 ± 2.3 km s⁻¹ for the secondary. Measuring RVs in the UV spectra allows to determine the terminal wind velocity (v_∞) and thus the mass-loss rate (Ṁ) more precisely (see Sect. 3.3.4). We estimated the RVs in the low resolution FOS spectrum such that the width of the oxygen lines in the range of 1330 Å–1420 Å can be matched, arriving at v₁ = 92.9 ± 6.2 km s⁻¹ for the primary and v₂ = 113.2 ± 2.8 km s⁻¹ for the secondary.

3.3. Stellar atmosphere modeling

To analyze the spectra, we employ the Potsdam Wolf-Rayet (PoWR) model atmosphere code (Gräfener et al. 2002; Hamann & Gräfener 2003, 2004; Todt et al. 2015; Sander et al. 2015). A short characterization of the code is given in Appendix B.

3.3.1. Temperatures and surface gravities

We measure the stellar temperatures using the ratio between He I to He II lines (see Fig. D.1). To constrain the surface gravities, we used the wings of the Balmer, Hβ, Hγ, and Hδ lines in the UVES and X-shooter (2022) spectra. For the primary we obtained an effective temperature of T₁ = 43 ± 2 kK, while the secondary is slightly cooler with T₂ = 42 ± 2 kK. The surface gravity is log g = 4.0 ± 0.2 for both binary components.

3.3.2. Stellar luminosities and spectroscopic masses

Luminosity L and color excess E_B − V are determined by fitting the composite spectral energy distribution (SED) – containing the synthetic flux of both stellar components – to photometry (see Fig. D.2). The color excess is modeled as a combined effect of Galactic foreground, using the reddening law by Seaton (1979) with E_{B − V,  Gal} = 0.03 mag, and SMC background, using the reddening law by Bouchet et al. (1985) with E_{B − V,  SMC} = 0.11 mag.

The luminosities of the two binary components are chosen such that the ratio of the depths of the synthetic helium lines match all spectra. Additionally, we use the C IIIλ1175 line in the UV, which is sensitive to changes in the light ratio. Both stars have comparable luminosities in the optical and UV. The final luminosities are log(L₁/L_⊙) = 5.41 ± 0.15 and log(L₂/L_⊙) = 5.38 ± 0.15 for the primary and secondary, respectively. Accordingly, the masses of the primary and secondary are $M_{\mathrm{spec},1}={32}_{-7}^{+8}{M}_{\odot }$ and $M_{\mathrm{spec},2}={31}_{-7}^{+8}{M}_{\odot }$. The spectroscopic mass ratio is q_spec = 0.97. When employing the method of Wilson (1941), one can derive an independent measure of the mass ratio from the RVs. This yields q_Wilson = 0.95 ± 0.02 which is in agreement with the spectroscopic results.

3.3.3. CNO surface abundances

All metal lines are well-matched with standard initial abundances of the SMC (see Table B.1 and Fig. D.3). The only line which is not well reproduced is the N IVλ3479 line, which is deeper than predicted by the synthetic spectrum. This absorption line is quite narrow, suggesting that the N IVλ3479 line might be blended with an unidentified ISM line. Pristine abundances prompt us to conclude that the AzV 14 components are young unevolved stars which have not yet interacted.

3.3.4. Wind mass-loss rates

In the optical spectra of AzV 14, Hα and He IIλ4686 do not show any contribution from winds. However, in the UV one can see wind lines with P Cygni profiles. Key diagnostic lines are the N Vλλ1239, 1243 and C IVλλ1548, 1551 resonance doublets and the He IIλ1640 line. By fitting the observed N V and C IV resonance doublets in the STIS and FOS spectra consistently yields a terminal wind velocity of v_∞ = 1600 ± 200 km s⁻¹ and a mass loss rate of log(Ṁ/(M_⊙ yr⁻¹)) = −7.7 for each of the stars (see Fig. D.4).

To match the O VIλλ1032, 1038 resonance doublet in the FUSE spectrum, we include in the model a hot plasma component with T_X = 3 MK emitting X-rays with L_X, mod = 2 × 10³⁰ erg s⁻¹. According to the Chandra Point Source Catalog v.2.0 (Evans et al. 2010), the X-ray flux of AzV 14 in the 0.5 keV–7.0 keV band is F_X = 7.3 × 10⁻¹⁶ erg s⁻¹ cm⁻², corresponding to an X-ray luminosity corrected for reddening (see Fig. D.2) L_X ≈ 5 × 10³² erg s⁻¹ or log(L_X/L_bol)≈ − 6. We suggest that the majority of the observed X-rays originate from a colliding wind zone between two binary components in common with O-type binaries in the Galaxy (Rauw & Nazé 2016). A summary of the empirically determined stellar parameters is given in Table 1.

4. Orbital analysis

The TESS light curve of AzV 14 displays periodic variability (originating from ellipsoidal variability, see Fig. 3). We employed a frequency analysis and phase folding techniques to determine the dominant periodicity in the TESS light curve, yielding an orbital period of P_orb = 3.7058 ± 0.0013 d (see Appendix C).

Fig. 3. Comparison of the phased observed and synthetic light and RV curve of AzV 14. Upper panel: phased TESS light curve (dots) and best fit obtained with the PHOEBE code (red line). For the PHOEBE model we show the 1σ and 2σ deviations as pink shaded areas. Lower panel: observed (triangles) and fitted RV curves (solid lines) of the primary (green) and secondary (red) component.

Orbital parameters from the light and RV curve are determined consistently with the Physics of Eclipsing Binaries (PHOEBE) code version 2.4.5 (Prša & Zwitter 2005; Prša 2011; Prša et al. 2016; Horvat et al. 2018; Jones et al. 2020; Conroy et al. 2020). The input parameters are provided by the spectral analysis (Table 1). To reduce the free parameter space, the orbital period is fixed to the period measured from the TESS light curve. Since the light curve is sinusoidal and has minima at phases Φ = 0.0 and Φ = 0.5, the binary orbit must be close to circular (e = 0). We fit the remaining orbital parameters, including the inclination. A more detailed description can be found in Appendix E.

The best fit is achieved with an inclination of i = 35 ± 5° and an epoch of the primary eclipse of T₀ = 2459036.101 ± 0.004 HJD. The orbital masses are $M_{\mathrm{orb},1}=33.6_{-3.7}^{+5.0}$ and $M_{\mathrm{orb},1}=31.9_{-3.5}^{+4.8}$, being in agreement with the spectroscopic analysis. The remaining orbital parameters are listed in Tables E.1 and E.2. The best fitting light and RV curves are shown in Fig. 3, including the phases at which each spectrum was taken. The X-shooter (2020) spectrum was obtained in the orbital phase close to conjunction, while the UVES and X-shooter (2022) spectra were measured at comparable phases out of conjunction.

Given the inclination of the system, the rotation velocities of the two components are v_rot = 157 ± 40 km s⁻¹. Hence, the rotation period of each star is P_rot = 3.0 ± 0.8 d which is close to the orbital period, implying that the binary is tidally locked.

Using the orbital period and the mass-ratio, we calculated the Roche radii for each star (see Table 1). By comparison to the previously determined stellar radii we conclude that both stars are underfilling their Roche lobe (R/R_RL = 0.6), further supporting our conclusion that both stars have not interacted yet.

5. Binary evolution models predict the formation of WR stars

Empirically derived stellar parameters of the binary components of AzV 14 (Table 1) are used to anchor the position of the system on the Hertzsprung–Russell diagram (HRD), and on this basis to investigate their possible future evolution. We calculated binary evolution models with the Modules for Experiments in Stellar Astrophysics (MESA) code (Paxton et al. 2011, 2013, 2015; Paxton et al. 2018, 2019). The methods and basic assumptions made in our models are described in Appendix F. In our models, a star enters WR evolutionary stages when the optical depth at its surface is τ ≥ 0.2. For further information we refer to Appendix G.

The current positions of the binary components in the HRD can be explained best by a binary evolutionary model with initial masses of M_1, ini = 35.0 M_⊙ and M_2, ini = 33.5 M_⊙, and an initial period of P_orb, ini = 3.7 d. The model predicts that AzV 14 formed 3.3 Myr ago. The corresponding evolutionary tracks are shown in Fig. 4, and the best fitting stellar parameters are listed in Table 2. According to the evolutionary model the binary components have not interacted yet, but will exchange mass in the future. During the future mass-transfer event, it is predicted that the secondary will accrete about 15 M_⊙. Both binary components evolve successively into WR stars.

Fig. 4. Evolutionary tracks of the primary (left) and secondary (right) compared to the empirical positions of AzV 14 (triangles with error ellipses). The evolutionary tracks are color-coded by the surface hydrogen abundance of the model and are overlayed by black dots corresponding to equidistant time-steps of 0.3 Myr. In the background, single star tracks are shown as dashed black lines. The instance of time in which the model is closest to the observations is marked by blue stars. The tracks are labeled by their initial stellar masses. The phases during which the model is expected to be observed as WR star, this is where τ ≥ 0.2, are highlighted by bold black frames. In addition, the positions of the WR stars AB 10 and AB 4, are marked by squares which are color-coded according to their observed surface-hydrogen abundances.

Both components will become H-poor WN type stars similar to WN stars observed in the SMC. For comparison the HRD shown in Fig. 4 includes the positions of the apparently single WR stars SMC AB 10 and AB 4. The primary will enter the WN stage at an age of 5.9 Myr and spends most of its WR lifetime (0.35 Myr) close to the helium zero-age main sequence, which is roughly at log(T_eff/K)≈5.0 (SMC AB 10 matches this position). During the WR stage, the primary will have a mass of 17 M_⊙, while being accompanied by a 50 M_⊙ main sequence star in an ≈8.5 d orbit. The mass ratio has only a second-order effect on the amount of mass removed during the mass-transfer phase (e.g., Pauli et al. 2022a, their Sect. 3.2), meaning that the primary would be located (if it goes through a stable mass-transfer event) at a comparable position in the HRD, even if it had a less massive companion that could prevent or complicate a detection. We presume that at the end of its evolution, the primary will directly collapse into a BH and remains bound.

The remaining secondary star, will continue its evolution on the main sequence. Shortly after that stage, it will initiate mass-transfer onto the BH, stripping off parts of its H-rich envelope and entering the WR phase at an age of 6.9 Myr. The resulting secondary WR star will spend most of its lifetime (also 0.35 Myr) at much lower temperatures (log(T_eff/K)≈4.75) than the primary WR star while being similarly massive with 25 M_⊙ (SMC AB 4 is at this position). This relatively cool WR star is accompanied by a BH with a mass of 16 M_⊙ in an orbit of 3.3 d.

Our binary evolutionary models of AzV 14 predict the formation of WR stars with different temperatures. Indeed, the WR population in the SMC (Hainich et al. 2015; Shenar et al. 2016, 2018) has a bimodal temperature distribution which is comparable to the predicted temperature regimes of the evolutionary models of AzV 14 (see Fig. 5). In the following, we elaborate on this idea further and discuss the implications and robustness of our findings.

Fig. 5. Observed temperature distribution of the WR stars in the SMC (gray area) compared to the time our evolutionary model of the primary (blue) and secondary (orange) are predicted to spend in the different temperature ranges during their WR phase. The observations are shown as Gaussians that have standard deviations corresponding to the observational uncertainties (see Table G.1). We excluded the binary SMC AB5 from this plot, as it has a different evolutionary origin.

6. Discussion

6.1. Exploring the parameter space by computing a small grid of binary evolutionary models

Inspired by the insights gained from the binary evolutionary models of AzV 14 (see Sect. 5), we explore if the bimodal temperature distribution of WR stars at low metallicity is also predicted by binary evolutionary models with different initial masses and initial orbital periods. Therefore, we calculated a small grid covering initial primary masses of 30, 45, and 60 M_⊙, while having a fixed mass ratio of q = 0.85. The initial primary masses are chosen to roughly represent the full luminosity distribution of the observed WR population of the SMC. The models are calculated for initial orbital periods of 5, 50, and 500 d. A HRD containing the evolutionary tracks of all the models is shown in Fig. 6.

Fig. 6. Evolutionary tracks of primary (left) and secondary (right) models with initial primary masses of 30 (green gray), 45 (blue gray), and 60 M⊙ (red gray), fixed mass ratio q = 0.85 and initial orbital periods of 5 (solid), 50 (dashed), and 500 d (dotted). We highlighted the WR phase (i.e., τ ≥ 0.2) in bold colors. Equidistant time steps of 0.3 Myr are marked by black dots during the WR stage. In the background, we marked the positions of all observed WR stars in the SMC.

Figure 7 depicts a simplified picture of the different formation channels of hotter and cooler WR stars. The flow chart is based on those models in our small grid that have a stable mass-transfer phase.

Fig. 7. Sketch of the evolutionary stages for possible formation channels of hotter and cooler WR stars under the assumption of stable mass transfer.

The evolutionary tracks of the primaries with initial masses of 30 M_⊙ are similar to those of AzV 14. All of the primary stars expand and initiate mass-transfer events, during which they lose most of the H-rich envelope, resulting in the formation of hotter WR stars (log T_eff > 4.9). The evolution of the secondaries depends on the initial parameters. The secondary in the system with initial period of 5 d accretes about ∼5 M_⊙ of material, leading to a rejuvenation of the core. After hydrogen is depleted in its core, the stellar model expands and quickly initiates mass transfer, stripping off parts of the accreted envelope. After the mass-transfer event the secondary has a temperature of log T_eff ≈ 4.7 and a high surface H-abundance of X_H = 0.5 (left side of Fig. 7). On the other hand, the secondary in the system with initial period of 50 d accretes less than 0.5 M_⊙. This is a big difference to the short-period model, which is explained by the fact that in our models accretion is only allowed when the accretor can avoid rapid rotation. In the short-period binaries, the stars are tidally locked, slowing down the rotation sufficiently for the accretor to stay below critical rotation during the accretion process. In a long-period binary, the star spins up quickly, limiting the amount of mass that can be accreted (Petrovic et al. 2005; de Mink et al. 2007; Shao & Li 2016). In such long-period systems (right side of Fig. 7) the secondary initiates mass transfer after core-H burning during its evolution toward the blue (BSG) and yellow supergiant (YSG) phase, leading to the formation of a hotter WR star with log T_eff ≈ 4.9 (i.e., the same temperature regime that is also populated by the primary models). In even longer period systems (P = 500 d) mass transfer is initiates when the star is evolving toward a red supergiant (RSG). The model has already formed a large convective envelope which expands adiabatically, making mass-transfer unstable. Potentially, this leads to a common envelope evolution which is not modeled here.

The models with initial primary masses of 45 M_⊙ and initial periods of 5 and 50 d also predict the formation of hotter WR stars (log T_eff > 4.85). However, the primary model with initial period of 500 d is on its way to becoming a RSG and has already formed a large outer convection region. Similar to the models presented above, this leads to an unstable mass transfer and possibly a common envelope evolution. In the model with an initial period of 5 d the secondary accretes a significant amount of mass (> 10 M_⊙) and initiates mass transfer after its main sequence evolution. After the mass transfer event the surface H fraction drops to X_H = 0.35 and the temperature is log T_eff > 4.7. The secondary will end its life as a cooler WR (left side of Fig. 7). For an initially wider orbit (50 d) the secondary accretes only negligible amounts of material and, similar to the primary models, evolves into a hotter WR star (right side of Fig. 7).

The primary models with initial masses of 60 M_⊙ populate a wide range of effective temperatures (log T_eff = 4.65 − 5.0). The broad temperature range can be explained by the changing efficiency of envelope stripping and its dependence on initial period and mass ratio. The efficiency is linked to the point in the primary’s evolution when it transitions into the core-He burning stage, which happens before or during the mass transfer (see also Klencki et al. 2022, their Fig. 3). In the primary model of the system with initial period of 5 d the mass transfer is very efficient in removing the H-rich envelope and with the help of the strong WR winds the model is able to remove a significant amount of the envelope, making the star appear hot. On the other hand, in the primary model with initial period of 500 d the mass-transfer is less efficient, making the WR star appear cooler. The secondary in the system with an initial period of 5 d accretes a significant amount of mass (> 10 M_⊙). It initiates mass transfer after the main sequence, stripping off parts of the accreted envelope, resulting in the formation of a cooler WR. The secondary in the system with an initial period of 50 d behaves similarly to its corresponding primary model and forms a hotter WR star. The secondary in the system with initial period 500 d initiates mass transfer when it is on its way to becoming a RSG and has already formed an extended outer convection zone. Mass-transfer in this model is unstable.

6.2. Bimodal temperature distribution of WR stars at low metallicity

From considering different evolutionary pathways probed by our exploratory model grid, we learned that primaries and secondaries which have not accreted a significant amount of material evolve to hotter WR stars. Only secondaries that have accreted a fair amount of mass (i.e., those with the shortest periods) evolve into cooler WR stars. There are two main factors that may be responsible for this behavior. Firstly, evolved accretors are characterized by a steep chemical gradient at the core-envelope boundary, producing cooler WR stars once the outer envelope is lost (Schootemeijer & Langer 2018, their Fig. 9). Secondly, accretors that are not fully rejuvenated tend to remain compact after the end of MS and begin the core-He burning phase as BSGs (Vanbeveren et al. 2013; Justham et al. 2014). Mass transfer from such stars leads to less efficient envelope stripping and cooler WR products (Klencki et al. 2022). This is different from the single star models (e.g., Georgy et al. 2012; Choi et al. 2016; Eldridge et al. 2017; Limongi & Chieffi 2018) which can only explain the hottest and most luminous WR stars.

Based on our models, we hypothesize that there must be a bimodal temperature distribution of faint WR stars at low metallicity. First, hotter WR stars of any luminosity can be well explained by the primary models as well as by the secondary models that did accrete negligible amounts of material (i.e., those in wide orbits). Second, cool and faint WR stars (log T_eff = 4.65 − 4.7 and log(L/L_⊙)≃5.5 − 5.9) must arise from secondaries that have accreted a significant amount of material, leading to a rejuvenation and high envelope-to-core mass ratios. Third, hotter and luminous WR stars (log T_eff = 4.65 − 4.7 and log(L/L_⊙)≃5.9 − 6.2) can be explained by primary, secondary, and single star models. We note that according to our models the cool and faint WR stars should all be accompanied by a compact companion, which can avoid detection and can help to explain observed apparently single cooler WR stars. We sketched the morphology of the WR population in the SMC in Fig. 8.

Fig. 8. HRD containing the positions of the observed WR stars in the SMC, over-plotted by a sketch stating the different evolutionary channels of single stars and binaries undergoing stable mass transfer, leading to the different marked regions.

It is less clear whether cooler WR stars at higher metallicities are formed in the same way as in low-metallcity galaxies. Two reasons should be considered. First, at high metallicity, the observed effective temperatures of WR stars can be lower, due to stronger winds (e.g., Sander & Vink 2020; Sander et al. 2023), also termed “dynamic inflation” (Grassitelli et al. 2018). Second, at higher metallicity WR stars additionally suffer from the effect of hydrostatically inflated envelopes (e.g., Gräfener et al. 2012; Sanyal et al. 2015), yielding again lower effective temperatures. Third, the stellar winds in the pre-WR stage at higher metallicity are stronger, stripping off more of the H-rich envelope, which is left after mass transfer. While the latter leads to the formation of hotter WR stars (as shown in Sect. 6.4.1), the effect on the observed temperature might partially be counter-balanced by the higher wind densities. All these effects make it difficult to distinguish between WR stars formed from stars which have accreted significant amounts of material in the past and intrinsically inflated WR stars.

Hence, in order to clearly see a bimodal temperature distribution of WR stars originating from post interaction binaries, populations at low metallcity must be considered. However, the number of WR stars decreases with decreasing metallicity (e.g., Shenar et al. 2020a), enforcing us to rely on small number statistics. In the SMC, the sample of WR stars is complete which minimizes selection biases. Therefore, despite the small number statistics, it remains the best representative sample of an WR population at low metallicity. In order to confirm or falsify our predictions, further observations of complete populations of WR stars in other low metallicity galaxies in combination with calculations of binary evolutionary models are necessary.

Our findings and other recent results (Renzo & Götberg 2021; Renzo et al. 2023) indicate that evolution of past accretors may be systematically different to those of normal (primary) stars. This impacts our understanding of stellar evolution and feedback. For instance, in the SMC hotter WR stars like AB 10 (log(Q_He II/s⁻¹) = 48.1), have strong He II ionizing fluxes, while cooler WR stars like AB 4 (log(Q_He II/s⁻¹) = 37.5) have five orders of magnitude lower ionizing flux. Broader implications for binary evolution channels, proposed in this work, are yet to be explored.

6.3. Observational fingerprints of previous evolutionary channels

Our model predicts that the observed effective temperature is one of the key diagnostics to differentiate between the evolutionary origin of WN type stars at low metallicity, as explained in Sect. 6.2. However, for WN stars originating from past accretors, a low effective temperature is not the only diagnostic.

The surface H-abundance can be considered as another clue, as all of our rejuvenated secondary model show X_H ≳ 0.3, while for their corresponding primary models the surface H-abundance are noticeably lower X_H ∼ 0.2. However, for our most luminous models (M_ini,  1 = 60 M_⊙) in the systems with the widest orbits, the primaries have low temperatures and the H-abundances are comparable to those of the rejuvenated secondaries. Hence, a high hydrogen abundance and a low temperature are not a robust criterion to identify the previous evolutionary path.

In our search for additional fingerprints to detect past accretors, we identified a difference in the predicted CNO surface composition of normal (primary) stars and the rejuvenated secondaries. It is expected that both, the primaries and nonrejuvenated secondaries, have CNO equilibrium composition (X_C = 1 − 2 × 10⁻⁵, X_N = 139 × 10⁻⁵, and X_O = 1 − 3 × 10⁻⁵). In a direct comparison, the models of the rejuvenated secondaries have a decreased surface N-abundance, while the surface O-abundance is increased by about one order of magnitude (X_C = 1 − 2 × 10⁻⁵, X_N ≈ 120 × 10⁻⁵, and X_O = 10 − 25 × 10⁻⁵), due to the formation of a shallow chemical gradient extending through the accreted envelope. Therefore, we predict that oxygen lines should be present in spectra of cooler WN stars formed from past accretors.

For the analysis of larger WN samples (e.g., Hamann et al. 2006; Sana et al. 2014; Hainich et al. 2015; Shenar et al. 2016), traditionally the N-abundance has been fixed to a value reflecting CNO equilibrium, while oxygen has not been included at all, due to its weak or absent spectral imprint in the optical regime. To predict the spectral appearance of WN stars which originate from stars that have accreted significant amounts of material in the past, synthetic spectra now need to include oxygen as well. Given the temperatures of WR stars O IV and O V are expected to be dominant ions.

In order to check this prediction, we select the WR star SMC AB 4 (WN6h), as its empirical position in the HRD can only be explained by our models of secondaries that have accreted material during mass transfer. A spectrum of AB4 taken with X-shooter at the ESO-VLT within the ESO program 106.211Z (PI: J. S. Vink) is available from the archive. We re-calculated the stellar atmosphere model of AB 4, from Hainich et al. (2015) with an enhanced O-abundances of X_O = 20 × 10⁻⁵ as well as with the CNO equilibrium value of X_O = 2 × 10⁻⁵. Indeed, the model with ten times enhanced oxygen abundance reproduces the O IV lines at wavelengths 3403 Å and 3414 Å just as observed, while the CNO equilibrium abundance fails by far (see Fig. 9).

Fig. 9. Close-up view of the O IV multiplet in the X-shooter spectrum of the WR star SMC AB 4. The observed spectrum is shown as a blue solid line, the synthetic spectra are calculated with an oxygen abundance of XO = 2 × 10−5 and XO = 20 × 10−5 and are shown as green dotted and red dotted lines, respectively.

As illustrated by this example, WN stars with enhanced oxygen abundance, such as AB 4, may be explained by the evolutionary models assuming that they accreted a significant amount of matter in the past. We suggest that the oxygen abundance of WN type stars is an important diagnostic which can help to identify WR stars originating from binary components that have accreted significant amounts of material in the past.

6.4. Sensitivity of the evolutionary model to different parameters

Estimating uncertainties in stellar evolution models is a nontrivial task, as it does not only depend on the numerical uncertainty but also on the physical assumptions made within the model itself. Stellar evolution models rely on many free parameters that can alter the evolution of a star, for instance mass-loss and mixing processes. Exploring the impact of all free parameters, quickly opens a large parameter space and needs intensive computational power. Such an investigation is far beyond the scope of the present paper. In the following, we briefly discuss two examples of how these parameters can impact our predictions on the temperature distribution of WR stars in low-metallicity environments.

6.4.1. Example 1: Enhanced mass-loss rates

In the binary model presented in this paper, the mass-transfer takes place while the primary is still core-H burning. Therefore, the evolution of the secondary is not really sensitive to the mass-loss rates of O stars. However, both stars spend a significant amount of their lifetime in the WR stage. Hence, we decided to test how the evolution would be affected if the WR mass-loss rate is increased by a factor of 3.

Figure 10 shows the resulting tracks. By comparison to Fig. 4, the predicted temperature regimes and the surface chemical composition during the WR phases changes. The primary loses a larger fraction of its H-poor envelope and spends more time at higher temperatures. At the end of core-He burning the stellar evolution model starts to contract instead of expanding. On the other hand, the stellar evolution model of the secondary shows even stronger response to the increased mass-loss rate. Instead of becoming a cooler WR star with plenty of hydrogen in the envelope, it is now able to remove large fractions of this envelope and position itself at higher temperatures in the HRD. We note that the region in the HRD where the stellar evolution model of the secondary is located is not populated by an observational counterpart.

Fig. 10. Same as Fig. 4, but now with the WR mass-loss rates enhanced by a factor 3.

The above example shows that the choice of the mass-loss recipe during different evolutionary phases can have a drastic impact on the evolution of massive stars. Only when using adequate mass-loss recipes, one can explain the cooler WR stars observed in the SMC. We note that in the case of WR winds being weaker than we assume in this work, the prediction on the morphology of the temperature distribution of WR stars would shift to lower metallicities.

6.4.2. Example 2: More efficient mixing

There are several mixing processes within stellar interiors, including convective mixing, semiconvective mixing, overshooting, rotational mixing, pulsational mixing, etc. Typically, these mixing processes are parameterized by free parameters, for example the mean free path a photon can travel, or the efficiency of a specific mixing process which cannot be predicted by theory. In the literature there are many works trying to limit the parameter space of the different free parameters (e.g., Schootemeijer et al. 2019; Higgins & Vink 2019; Gilkis et al. 2021; Michielsen et al. 2021). In this section, we want to showcase how changing one of the free parameters can impact the evolution of a star.

One of the least constrained, yet important, mixing process is the efficiency of semiconvective mixing. Schootemeijer et al. (2019) used the population of the BSGs and RGSs to constrain semiconvective mixing efficiencies, but were able only to limit it to α_sc ≳ 1. They report that higher values only barely impact the number ratio of blue and red supergiants in galaxies. In our models presented above we used α_sc = 1. For the models presented in this section, however, we increased the efficiency to α_sc = 10 in order to see how it changes our understanding of the evolution of the binary.

The evolutionary tracks of the primary and secondary with the more efficient semiconvective mixing are shown in Fig. 11. By comparing these tracks to those shown in Fig. 4, one can see that the evolution of the primary barely changes. Semiconvective mixing becomes important only in evolutionary stages after the main sequence. The primary begins mass transfer already during its main sequence evolution and spends its post-main sequence evolution as a WR star, leaving the semiconvective mixing no time to impact the evolution. On the other hand, the evolution of the secondary changes because during the accretion phase it forms extensive convective and semiconvective regions in its envelope. The secondary initiates mass-transfer after the main sequence evolution, giving semiconvection some time to efficiently mix material. In addition, due to more efficient semiconvection, the secondary expands less in the transition to the core-He burning phase and initiates mass transfer at a somewhat more advanced evolutionary stage, consequently losing less mass during mass-transfer and producing a cooler WR star with an even higher surface H abundance. This is reflected in the higher surface H-abundance and the cooler temperature. Klencki et al. (2022) recently reported something similar: Binary models at low metallicity with highly efficient semiconvective mixing will not become WR/helium stars after Case B mass transfer. Their models managed to stay at the position in the HRD at which they initiated the mass transfer and hence are core-He burning stars, disguised as O- or B-type stars. It is worth mentioning that according to our criterion based on the optical depth, the secondary of the binary evolutionary model with efficient semiconvective mixing would spend only a short time as a WR star at low temperature.

Fig. 11. Same as Fig. 4, but now with an increased efficiency of semiconvection.

These two examples illustrate that binary evolution models are quite sensitive to the input parameters. A more detailed study based on grids of detailed binary evolutionary models and population synthesis are needed to confirm or disprove our predictions on the morphology of WR populations at low metallicity.

7. Summary and conclusions

A consistent analysis of multi-epoch optical and UV spectra of one of the earlierst O-type stars in the SMC, AzV 14, reveals its binary nature. Furthermore, our analysis uncovered that the systems’ two components are very similar. The primary and secondary have temperatures of T_eff,  1 = 42.8 ± 2 kK and T_eff,  2 = 41.8 ± 2 kK, luminosities of log(L₁/L_⊙) = 5.41 ± 0.15 and log(L₂/L_⊙) = 5.38 ± 0.15 and surface gravities of log(g_1, 2/(cm s⁻²)) = 4.0 ± 0.2, respectively. The analysis of a TESS light curve confirms the binary nature of AzV 14 with an orbital period of 3.7058 ± 0.0013 d. Their spectroscopic masses of M₁ = 32 ± 8 M_⊙ and M₂ = 31 ± 8 M_⊙ are confirmed by the orbital analysis. Both binary components drive weak stellar winds with mass-loss rates of about log(Ṁ/(M_⊙ yr⁻¹)) = −7.7 ± 0.2. The relatively high observed X-ray emission of AzV 14 is attributed to colliding winds.

The new empirically derived stellar parameters and the current orbital period of AzV 14 are well explained by our evolutionary models. In particular, the spectroscopic and orbital masses of the binary components are in agreement with the evolutionary ones. According to the evolutionary models, the components of AzV 14 did not yet exchange mass. Most interestingly, the binary evolutionary model of AzV 14 predicts that the primary will evolve into a hotter WR star, while the secondary is destined to evolve into a cooler WR star.

Inspired by these results we calculated a small evolutionary model grid to investigate the conditions for the formation of hotter and cooler WR stars. Guided by these calculations, we anticipate that WR populations in low metallicity galaxies should show bimodal temperature distributions. According to our models, stars born in long-period binaries both evolve to hotter WR stars. On the other hand, for stars born in short-period binaries only the primary is destined to become a hot WR star, while the secondary evolves into a cooler WR star. In our models we assume that accretion in short-period binaries is efficient, leading to a rejuvenation of the core and a different envelope-to-core mass ratio. This eventually allows the past accretors to stay at lower temperatures during their WR stage.

These results are sensitive to the input physics, but if turned out to be correct they would have a vast impact on our understanding of binary evolution at low metallicity, stellar feedback, and, hence, on galaxy formation and evolution. To test if our models are reliable and cooler WR stars originate from past accretors we introduced an additional criterion: the surface oxygen abundance. From standard evolutionary models one expects that a WR star has a surface composition which corresponds to the CNO equilibrium value. However, in our models the cooler WR stars, formed by stars that accreted a significant amount of material in the past, all have surface oxygen abundance which are increased by a factor of about ten. We tested this prediction with archival spectra of the apparently single WN star AB4 in the SMC. We detect oxygen lines in its optical spectrum, and measure the high oxygen abundance which is in agreement with our evolutionary model predictions. This empirically supports the new evolutionary pathway to the formation of WR stars at low metallcity proposed in this paper.

Acknowledgments

The authors are appreciative to the reviewer for useful comments which helped to improve the paper, and for their suggestions. The results presented in this paper are based on observations obtained with the NASA/ESA Hubble Space Telescope, retrieved from MAST at the Space Telescope Science Institute (STScI). STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555. Support to MAST for these data are provided by the NASA Office of Space Science via grant NAG5-7584 and by other grants and contracts. The TESS data presented in this paper were obtained from MAST at the STScI. Funding for the TESS mission was provided by the NASA Explorer Program. Furthermore, its conclusions are based on observations collected at the European Southern Observatory (ESO) under the program 098.A-0049. The authors thank the managing committee of XShootU and Andrea Mehner for preparing the OBs of the XShootU project. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. D.P. and S.R.S. acknowledge financial support by the Deutsches Zentrum für Luft und Raumfahrt (DLR) grants FKZ 50OR2005 and 50OR2108. A.A.C.S. and V.R. acknowledge support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) in the form of an Emmy Noether Research Group – Project-ID 445674056 (SA4064/1-1, PI: Sander). D.M.B. gratefully acknowledges a senior postdoctoral fellowship from the Research Foundation Flanders (FWO) with grant agreement number 1286521N. R.K. acknowledges financial support via the Heisenberg Research Grant funded by the German Research Foundation (DFG) under grant no. KU 2849/9. C.K. acknowledges financial support from the Spanish Ministerio de Economía y Competitividad under grants AYA2016-79724-C4-4-P and PID2019-107408GB-C44, from Junta de Andalucía Excellence Project P18-FR-2664, and from the State Agency for Research of the Spanish MCIU through the ‘Center of Excellence Severo Ochoa’ award for the Instituto de Astrofísica de Andalucía (SEV-2017-0709). T.S. acknowledges support from the European Union’s Horizon 2020 under the Marie Skłodowska-Curie grant agreement No. 101024605. The collaboration of coauthors was facilitated by support from the International Space Science Institute (ISSI, Bern). The authors are grateful to the Lorentz Center (Leiden) for organizing a stimulating workshop. We thank Paul Crowther for useful discussions as well as advice regarding the UVES spectrum and for providing the updated IR photometry.

References

Appendix A: Spectroscopic and photometric data

The information on all observed spectra used in this work is given in Table A.1. All spectra are in a reduced state from the individual observatories pipeline. Apart from the UVES spectrum, all other spectra analyzed in this work had been corrected for barycentric motion by the automated extraction pipeline. For the UVES spectrum we had to calculate this correction ourselves. With the help of the tool from Wright & Eastman (2014) we obtained v_bary = −10.6 km s⁻¹. A summary of the estimated RVs and the calculated orbital phases of the individual observed spectra that are used in our analysis are listed in Table A.2. The photometric data employed for the SED fitting are listed in Table A.3.

Appendix B: Spectral modeling with the PoWR stellar atmosphere code

PoWR assumes that the atmosphere, including the stellar wind, is stationary, spherically symmetric, energy conserving, and in radiative equilibrium. Hot star atmospheres cannot be approximated by LTE. Instead, the equations of radiative transfer and the statistical equilibrium have to be solved consistently in the comoving frame. This is achieved by iteration with the “approximate lambda operator” technique, yielding population numbers of the individual elements.

The stellar atmosphere models which we apply here account for detailed model atoms of the elements: H, He, C, N, O, Mg, Si, P, S. The iron group elements Fe, Sc, Ti, V, Cr, Mn, Co, and Ni are treated as one generic element (“G”) using a superlevel approach (see Gräfener et al. 2002). The individual abundances for each element and the considered ionization levels adopted here are listed in Table B.1.

To calculate the detailed synthetic spectrum we assume a depth dependent micro-turbulence velocity. In our model calculations it is assumed that the micro-turbulence velocity starts with ζ_ph = 10 km s⁻¹ in the photosphere and grows linearly with the wind velocity as 0.1 v(r).

The density stratification in the photosphere is calculated from the hydrostatic equations, taking the radiation pressure fully into account (Sander et al. 2015). The velocity field in the supersonic regions, is adopted as a β-law (Castor et al. 1975) with the standard exponent of β = 0.8 (Pauldrach et al. 1986).

The stellar parameters, namely effective temperature T_eff, surface gravity g, wind mass-loss rate Ṁ, terminal wind velocity v_∞, as well as the chemical composition, are determined by a comparison of the synthetic spectrum to the observation. The effective temperature T_eff is defined as referring to the radius where the Rosseland mean continuum optical depth is τ = 2/3.

Stellar winds are inhomogeneous; this is taken into account using the approximation for optically thin clumps (microclumping). The clumping factor D describes by how much the density within the clumps is increased compared to a smooth wind with the same mass-loss rate (Hamann & Koesterke 1998). For our analysis, we adopt depth-dependent clumping that starts at the sonic radius and increases outward, until a clumping factor of D = 10 is reached at a radius of 10 R_*.

Appendix C: Extraction of the TESS light curve

AzV 14 was observed with TESS, but has no automated pipeline processed light curve yet. Hence, we followed the method of Bowman et al. (2022) and extracted pixel cutouts of size 20 × 20 using ASTROCUT (Brasseur et al. 2019) provided by MAST² to extract a custom light curve and search for photometric variability. We tried differently sized aperture masks for an optimum light curve using the LIGHTCURVE software package (Lightkurve Collaboration 2018). Since AzV 14 is a relatively faint target (V = 13.9 mag) for TESS and because the pixel size of TESS is 21″, the aperture mask that maximized the signal-to-noise of the variability in the light curve and minimized the contribution of contaminating stars was found to be only a single pixel. We note that the field around AzV 14 is quite crowded and that the employed TESS pixel contains several other stars. However, all stars identified by Gaia within this single neighboring pixels have V ≳ 16.5. This is three orders of magnitude fainter than our target. None of this stars is known to be variable and combined they contribute less than < 5% of the total measured flux. Pixel masks larger than a single pixel had larger percentage contamination and diminished signal from our target, thus justifying our extraction method. The background flux was estimated by calculating the median observed flux per frame whilst excluding pixels that predominantly contain stellar flux. The background flux was subtracted, and then we normalized the extracted light curve by dividing through the median flux.

TESS sector 1 data are known to be quite noisy and contain problematic systematics. For such a faint target this presented a significant challenge in extracting a reliable light curve. However, the light curves from sectors 27 and 28 are judged to be robust, and span a combined total of 50.8 d.

The variability in the combined sectors 27 and 28 of the TESS data for AzV 14 is well described as periodic. To determine the dominant periodicity we identified the highest-amplitude peak in the amplitude spectrum of the light curve. Next, we optimized the frequency and amplitude and determined uncertainties for this dominant periodicity by fitting a sinusoid to the light curve using nonlinear least-squares methods with the PERIOD04 software tool (Lenz & Breger 2005). We determined the dominant frequency to be 0.53969 ± 0.00019 d⁻¹, corresponding to a period of 1.85292 ± 0.00065 d. Yet, as demonstrated from our spectroscopic analysis (see Sect. 3.1), AzV 14 contains two similar massive stars. Hence, periodic sinusoidal variability in such a system implies grazing eclipses or ellipsoidal variability. In phase space, however, one expects a symmetric double wave for such variability, meaning that the true dominant frequency is 1.07938 ± 0.00038 d⁻¹, and thus the orbital period of the system is P_orb = 3.7058 ± 0.0013 d. The phase-folded light curve on this inferred orbital period is shown in Fig. 3. We also find this frequency to be significant in the TESS data, but it has a smaller amplitude, hence is not the dominant frequency. It is not uncommon for the (sub)harmonic of the true orbital frequency to have a larger amplitude in the amplitude spectra of (eclipsing) binary star or ellipsoidal variable light curves (see e.g., IJspeert et al. 2021; Southworth & Bowman 2022).

TESS has a low spatial resolution from its large pixel sizes, which when combined with our small aperture mask results in a relatively large scatter in the light curve. Nonetheless, it is possible to use the TESS data to constrain the inclination of the binary. In turn, this allows us to compare inclination dependent parameters, such as the projected rotation velocities to predictions of binary evolutionary models.

Appendix D: Additional figures

Fig. D.1. Hydrogen and helium lines in the two X-Shooter spectra and the UVES spectrum (blue) compared to the combined synthetic spectrum (red). The individual unweighted synthetic spectra of the primary and secondary components are represented as black dotted and green dotted lines, respectively.

Fig. D.2. SED comparing flux-calibrated observations (blue lines) and photometric data (open squares) to the model SED composed of both binary components (red dashed line). The individual synthetic SEDs of the primary and the secondary are shown as thin black and thick green dashed lines, respectively. Due to their similar luminosities these curves are hard to distinguish.

Fig. D.3. Close-up view of the metal lines in the X-Shooter (2020) and UVES spectra. The observed spectrum is shown as a solid blue solid line, the individual unweighted synthetic spectra of the primary and the secondary component are represented as black dotted and green dotted lines, respectively and the combined synthetic spectrum as red dashed line. In the combined synthetic spectrum ISM absorption lines, which originate from the Galactic foreground and in the SMC, are included. They are labeled as “(is)”.

Fig. D.4. Close-up view of the wind diagnostic lines in the STIS and FOS spectra, as well as Hα, which is part of the X-Shooter (2020) spectrum. The observed spectrum is shown in blue and needs to be compared to the combined synthetic spectrum, shown in red. In the combined synthetic spectrum ISM absorption lines, which originate from the Galactic foreground and in the SMC, are included and labeled as “(is)”. For the STIS and FOS spectrum we plot the weighted synthetic spectra. For the X-Shooter (2020) spectrum we show the unweighted synthetic spectra of the primary and the secondary as black dotted and green dotted lines, respectively. A close-up view of the He IIλ4686 wind diagnostic line can be found in Fig. D.1.

Fig. D.5. Close-up view of the important diagnostic lines in the far-UV FUSE spectrum. The observed spectrum is shown as a blue solid line, the weighted synthetic spectrum of the primary and secondary are shown as black dotted and green dotted lines and the combined synthetic spectrum is shown as a dashed red line. In the combined synthetic spectrum ISM absorption lines, which originate from the Galactic foreground and in the SMC, are included. They are labeled as “(is)”. We note that not all ISM lines are modeled, and that FUSE spectra are additionally contaminated by absorption lines from molecular hydrogen.

Appendix E: Light and RV curve modeling with the PHOEBE code

The PHOEBE code does not contain atmosphere tables for hot stars. Therefore, we approximated the AzV 14 spectrum as a blackbody. From comparison of a blackbody and our synthetic spectrum in the V-band, we find that this is a reasonable approximation. To model the effect of limb-darkening, we adopt a quadratic approximation for the intensity (Diaz-Cordoves & Gimenez 1992)

$$\begin{array}{c}\hfill I(\mu )=I(1)[1-a_i(1-\mu )-b_i{(1-\mu )}^2],\end{array}$$(E.1)

where μ = cos θ is the cosine of the directional angle and a_i and b_i are the linear and quadratic coefficients of the binary component i. The coefficients are determined by fitting the Eq. (E.1) to the emergent intensity distribution in our best fitting PoWR model. For the primary, we derive a₁ = 0.1474 and b₁ = 0.0985, and for the secondary, a₂ = 0.1489 and b₂ = 0.1025. The effect of gravitational darkening is approximated by a power-law with exponent β_grav = 1. For simplicity, we assume that both stars are rotating synchronously with the orbit, a reasonable approach for stars in binaries with orbital periods on the order of a few days. In our simulations, the effect of ellipsoidal variability, which is important in short-period binaries, is also included. The reflection effect (Wilson 1990) is modeled with two reflections.

As a fitting routine, we choose the built-in option of the emcee sampler (Foreman-Mackey et al. 2013), enabling us to consistently model the light and RV curve. From our PHOEBE model, we conclude that the observed light curve originates purely from ellipsoidal modulations, meaning that the flux during quadrature is the highest and lowest when the stars are aligned. The resulting parameters are listed in Table E.1 and E.2.

Appendix F: Input physics in the MESA binary evolutionary model

Our binary evolutionary models are calculated with the MESA version r15140. The chemical composition is adjusted to the SMC. Following Brott et al. (2011), we use tailored abundances for H, He, C, N, O, Mg, Si, and Fe in the evolutionary models. We note that these are similar to those used in the PoWR atmosphere models (see Table B.1). As initial parameters, we only set the initial masses of the binary components as well as the initial orbital period. We assume that the binary is initially tidally synchronized, a simplification adopted to reduce the free parameter space.

Our models include rotational mixing and treat it as a diffusive process, which encompasses the effects of dynamical and secular shear instabilities, the Goldreich-Schubert-Fricke instability, as well as Eddingtion-Sweet circulations (Heger et al. 2000). The efficiency coefficients of rotational mixing are adopted from Brott et al. (2011) and set to f_c = 1/30 and f_μ = 0.1. On top of that, we include the angular momentum transport via magnetic fields from a Taylor-Spruit dynamo (Spruit 2002).

We use the standard mixing length theory (Böhm-Vitense 1958) with a mixing length coefficient α_mlt = 1.5 and the Ledoux criterion to model convection. Our models include the effect of envelope inflation, which occurs in regions within the star that approach the Eddington limit (Sanyal et al. 2015). Semiconvective mixing is taken into account with an efficiency coefficient of α_sc = 1 (Langer et al. 1983; Schootemeijer et al. 2019) and for core-H burning we include step overshooting in such a way that the convective core is extended by 0.335H_P (Brott et al. 2011; Schootemeijer et al. 2019). Thermohaline mixing, which becomes important for late evolutionary stages and during the accretion phase of the companion star, is modeled with an efficiency coefficient of α_th = 1 (Kippenhahn et al. 1980).

The choice of mass-loss rates is crucial for stellar evolutionary models. In particular, the mass-loss rates during the WR stages are important for those binaries which avoid evolving toward the cool side of the HRD. In the literature a discussion is ongoing on the efficiency of the WR winds at low metallicity. For example, Sander & Vink (2020) argue that below some limiting metallicity the WR winds become inefficient in removing mass. This limiting metallicity is thought to be close to the metallicity of the SMC. Therefore, we adjust the WR mass-loss recipes such that they are adequate for the SMC metallicity. We include mass-loss as follows: for main sequence stars (X_H > 0.7) we use the mass-loss rates by Vink et al. (2001) (which are an order of magnitude higher than those empirically derived for AzV 14 in this work). For temperatures below the so-called bi-stability jump we use the maximum value of either Vink et al. (2001) or Nieuwenhuijzen & de Jager (1990). And for the WR phase, this is when the surface hydrogen abundance drops below X_H < 0.4, the recipe of Shenar et al. (2019, 2020b) is used. For surface abundances in the range of 0.7 > X_H > 0.4 we linearly interpolate between Vink et al. (2001) and Shenar et al. (2019, 2020b).

Mass transfer is modeled implicitly using the “contact” scheme of the MESA code. Mass transfer is treated as rotation-dependent mass accretion. This means, in our models mass transfer is conservative as long as the accretor can stay below critical rotation. As soon as it is rotating critically, we assume that all mass is lost from the system. In addition to that, we assume mass accretion from a stellar wind.

The model calculation of a binary component is stopped after core-He depletion. It is assumed that as soon as the primary model depletes helium in its core, it directly collapses into a BH (Dessart et al. 2011). This means that the orbit is not perturbed by a supernova explosion. Under the assumption that we are able to model mass-transfer from the former secondary to the BH produced by the former primary. Recently it was reported that the first supernova originating from a WR progenitor was observed (Gal-Yam et al. 2022). We estimate that in the case of a supernova event the mass of the BH in our models might be overestimated by a factor 2 to 4. However, the main parameter defining the size of the Roche lobe, and thus the start and end of mass transfer is the orbital separation and not the mass ratio (Eggleton 1983, their equation 2). We assume that the impact of this simplification on the evolution of the former secondary is negligible. During the mass-transfer phase, we assume that the BH can only accrete material with the Eddington accretion rate.

Appendix G: Identifying WR phases in stellar evolutionary models

We calculated the optical depth at the surface of a star in our evolutionary model calculation with the formula from Langer (1989)

$$\begin{array}{c}\hfill \tau =\frac{-\kappa \dot{M}}{4\pi R(v_{\infty }-v_0)}ln\left(\frac{v_{\infty }}{v_0}\right)·\end{array}$$(G.1)

In this equation, the opacity κ is approximated by the electron scattering opacity κ_es = 0.2 (1 + X_H) cm² g⁻¹, Ṁ is the wind mass-loss rate, R is the radius of the star, v_∞ is the terminal wind velocity, and v₀ is the expansion velocity at the base of the wind, which, for simplicity, is assumed to be constant v₀ = 20 km s⁻¹. Gräfener et al. (2017) found that v_∞ can be linked to the escape velocity v_esc by

$$\begin{array}{c}\hfill v_{\infty }=\alpha ·v_{\mathrm{esc}}=\alpha ·\sqrt{\frac{2GM}{R}(1-\mathrm{\Gamma })},\end{array}$$(G.2)

where Γ is the Eddington factor and α is an evolutionary dependent scaling factor. We assume α = 1.3 during the main sequence and WN phase. For factors used in the other phases we refer to Pauli et al. (2022a), their section 2.2.

Ideally the threshold to differentiate between optically thin and thick winds, He-stars or WR stars, is τ ≥ 1. However, in the calculation of the optical depth several simplifications are included. Therefore, a calibration to observations is needed. Aguilera-Dena et al. (2022) already did this for H-free WN stars and determined a threshold value of τ ≥ 1.45. On the other side, Pauli et al. (2022a) find that for H-poor WN stars a lower threshold might be needed to reproduce observations. In the SMC only H-poor WN stars (and one WO type star) are present. To find a threshold value that can reproduce the observations, we calculated the optical depths of all the H-poor WN stars in the SMC. Their stellar parameters and optical depths are listed in Table G.1. We find that all have an optical depth value of τ ≳ 0.2, which we set here as our new threshold.

All Tables

All Figures

	Fig. 1. False-color HST image composed of the images in the F475W (blue) and F657N (yellow) filters. The position of AzV 14 and the approximate size of the NGC 261 H II region are indicated by white circles.
In the text

Fig. 2. Comparison of selected lines in the UVES (solid blue) and the X-shooter (dashed pink and dotted yellow) optical spectra. The high resolution UVES spectrum is binned by 0.15 Å to match the X-shooter spectra. All spectra are convolved with a Gaussian having a FWHM = 0.05 Å for a better comparison. The lines in the UVES and X-shooter (2022) spectra are separated, while they are blended in the X-shooter (2020) spectrum. This provides a clear indication for binarity.

In the text

	Fig. 3. Comparison of the phased observed and synthetic light and RV curve of AzV 14. Upper panel: phased TESS light curve (dots) and best fit obtained with the PHOEBE code (red line). For the PHOEBE model we show the 1σ and 2σ deviations as pink shaded areas. Lower panel: observed (triangles) and fitted RV curves (solid lines) of the primary (green) and secondary (red) component.
In the text

Fig. 4. Evolutionary tracks of the primary (left) and secondary (right) compared to the empirical positions of AzV 14 (triangles with error ellipses). The evolutionary tracks are color-coded by the surface hydrogen abundance of the model and are overlayed by black dots corresponding to equidistant time-steps of 0.3 Myr. In the background, single star tracks are shown as dashed black lines. The instance of time in which the model is closest to the observations is marked by blue stars. The tracks are labeled by their initial stellar masses. The phases during which the model is expected to be observed as WR star, this is where τ ≥ 0.2, are highlighted by bold black frames. In addition, the positions of the WR stars AB 10 and AB 4, are marked by squares which are color-coded according to their observed surface-hydrogen abundances.

In the text

Fig. 5. Observed temperature distribution of the WR stars in the SMC (gray area) compared to the time our evolutionary model of the primary (blue) and secondary (orange) are predicted to spend in the different temperature ranges during their WR phase. The observations are shown as Gaussians that have standard deviations corresponding to the observational uncertainties (see Table G.1). We excluded the binary SMC AB5 from this plot, as it has a different evolutionary origin.

In the text

Fig. 6. Evolutionary tracks of primary (left) and secondary (right) models with initial primary masses of 30 (green gray), 45 (blue gray), and 60 M⊙ (red gray), fixed mass ratio q = 0.85 and initial orbital periods of 5 (solid), 50 (dashed), and 500 d (dotted). We highlighted the WR phase (i.e., τ ≥ 0.2) in bold colors. Equidistant time steps of 0.3 Myr are marked by black dots during the WR stage. In the background, we marked the positions of all observed WR stars in the SMC.

In the text

	Fig. 7. Sketch of the evolutionary stages for possible formation channels of hotter and cooler WR stars under the assumption of stable mass transfer.
In the text

	Fig. 8. HRD containing the positions of the observed WR stars in the SMC, over-plotted by a sketch stating the different evolutionary channels of single stars and binaries undergoing stable mass transfer, leading to the different marked regions.
In the text

	Fig. 9. Close-up view of the O IV multiplet in the X-shooter spectrum of the WR star SMC AB 4. The observed spectrum is shown as a blue solid line, the synthetic spectra are calculated with an oxygen abundance of XO = 2 × 10−5 and XO = 20 × 10−5 and are shown as green dotted and red dotted lines, respectively.
In the text

	Fig. 10. Same as Fig. 4, but now with the WR mass-loss rates enhanced by a factor 3.
In the text

	Fig. 11. Same as Fig. 4, but now with an increased efficiency of semiconvection.
In the text

	Fig. D.1. Hydrogen and helium lines in the two X-Shooter spectra and the UVES spectrum (blue) compared to the combined synthetic spectrum (red). The individual unweighted synthetic spectra of the primary and secondary components are represented as black dotted and green dotted lines, respectively.
In the text

	Fig. D.2. SED comparing flux-calibrated observations (blue lines) and photometric data (open squares) to the model SED composed of both binary components (red dashed line). The individual synthetic SEDs of the primary and the secondary are shown as thin black and thick green dashed lines, respectively. Due to their similar luminosities these curves are hard to distinguish.
In the text

Fig. D.3. Close-up view of the metal lines in the X-Shooter (2020) and UVES spectra. The observed spectrum is shown as a solid blue solid line, the individual unweighted synthetic spectra of the primary and the secondary component are represented as black dotted and green dotted lines, respectively and the combined synthetic spectrum as red dashed line. In the combined synthetic spectrum ISM absorption lines, which originate from the Galactic foreground and in the SMC, are included. They are labeled as “(is)”.

In the text

Fig. D.4. Close-up view of the wind diagnostic lines in the STIS and FOS spectra, as well as Hα, which is part of the X-Shooter (2020) spectrum. The observed spectrum is shown in blue and needs to be compared to the combined synthetic spectrum, shown in red. In the combined synthetic spectrum ISM absorption lines, which originate from the Galactic foreground and in the SMC, are included and labeled as “(is)”. For the STIS and FOS spectrum we plot the weighted synthetic spectra. For the X-Shooter (2020) spectrum we show the unweighted synthetic spectra of the primary and the secondary as black dotted and green dotted lines, respectively. A close-up view of the He IIλ4686 wind diagnostic line can be found in Fig. D.1.

In the text

Fig. D.5. Close-up view of the important diagnostic lines in the far-UV FUSE spectrum. The observed spectrum is shown as a blue solid line, the weighted synthetic spectrum of the primary and secondary are shown as black dotted and green dotted lines and the combined synthetic spectrum is shown as a dashed red line. In the combined synthetic spectrum ISM absorption lines, which originate from the Galactic foreground and in the SMC, are included. They are labeled as “(is)”. We note that not all ISM lines are modeled, and that FUSE spectra are additionally contaminated by absorption lines from molecular hydrogen.

In the text