© The Authors 2023

1. Introduction

Massive stars (M ≳ 8 M_⊙) have a pivotal role in the evolution of Galaxies. They influence their environments throughout their entire evolution via multiple feedback mechanisms. Core hydrogen-burning massive stars are hot, leading to ionisation feedback (Hollenbach & Tielens 1999; Matzner 2002) and forming spectacular H II regions. Their high ultraviolet (UV) luminosity drives powerful stellar winds, depositing mechanical feedback and heavy elements into the interstellar medium (Rogers & Pittard 2013). When a massive star evolves, it expands and becomes a red supergiant (RSG) and cools down. These stars still have strong winds, powered by radiation pressure on dust grains. Most massive stars are believed to end their lives with a supernova explosion that rapidly deposits substantial amounts of energy in the surrounding medium (Rogers & Pittard 2013). During these events, neutron stars (NSs; Baym et al. 2018; Vidaña 2018) or black holes (BHs; Orosz et al. 2011; Miller-Jones et al. 2021) are formed. All these feedback mechanisms make massive stars, which are the drivers of the chemical evolution of the interstellar medium (Burbidge et al. 1957; Pignatari et al. 2010; Thielemann et al. 2011; Kasen et al. 2017; Kajino et al. 2019) and the arbitrators of stellar formation; massive stars can initiate or abort star formation (Mac Low et al. 2005).

Massive stars preferentially form in multiple systems (Sana et al. 2012, 2013, 2014; Moe & Di Stefano 2017). The evolution of stars in binary systems differs from single star evolution as the two components can interact and exchange mass, for instance when one star in the binary (i.e. the more massive one) expands and fills its Roche lobe (Pols 1994; Vanbeveren et al. 1998; Nelson & Eggleton 2001; Smith 2014). This mass can then, depending on other factors such as rotation, accrete onto the companion star. Due to the fate of massive stars as either NSs or BHs, close massive star binaries are the progenitors of gravitational wave observations (Abbott et al. 2016, 2019). However, the final mass of a star dictates its fate (Heger et al. 2003), and in a binary system the final mass depends on the previous evolution and interaction (Shenar et al. 2020a; Massey et al. 2021; Pauli et al. 2022a).

The duration of mass transfer strongly depends on the evolutionary stage of the star (Paczyński 1967; Vanbeveren & Conti 1980; Wellstein & Langer 1999). If the expanding star is still core H-burning, it first drops out of nuclear equilibrium, making the first phase of mass transfer rapid as it takes place on the thermal timescale (fast Case A). As soon as the star regains nuclear equilibrium, it continues to expand on the nuclear timescale, overflowing its Roche lobe and transferring mass to the companion star (slow Case A). If the mass donor during the slow Case A mass transfer is less massive than the mass gainer, the systems are called Algol-like system (named for β Per, Paczyński 1971; Batten 1989; Li et al. 2022). Alternatively, stars that have finished core H-burning expand on the thermal timescale, filling their Roche lobe quickly, and lose a large fraction of their envelope (Case B). As most of the mass transfer events are on short timescales, it is practically impossible to catch a system during the fast Case A or Case B mass transfer phase. As slow Case A mass transfer occurs on the nuclear timescale, it is possible to find observational counterparts.

Metallicity (Z) has many effects on stellar evolution. For instance, stars are more compact and hotter, and line-driven winds get weaker for lower metallicity stars (Mokiem et al. 2007). Low-Z environments are of particular interest when looking to understand the origin of gravitational waves as the observed mergers take place in high-redshift galaxies. However, it is only possible to study low-Z stars in nearby irregular dwarf galaxies, such as the Small Magellanic Cloud (SMC, Z_SMC = 1/7 Z_⊙) as these stars can be resolved. To date, there are few detailed studies on massive stars in pre- and post-interaction binaries in the SMC (de Mink et al. 2007; Pauli et al. 2022b).

In this paper we remedy this by studying in detail one of the most massive binaries in the low-Z SMC galaxy, namely SSN 7 (Sabbi et al. 2007, also known as MPG 435; Massey et al. 1989, or NMC 26; Niemela et al. 1986). Our target is located within the core of the young cluster NGC 346 (Dufton et al. 2019; Fig. 1), which contains a large population of O stars (Niemela et al. 1986; Massey et al. 1989; Rickard et al. 2022). The cluster is nearby (d = 61 kpc; Hilditch et al. 2005) and has been studied several times with different instruments. Previous studies have categorised SSN 7 as a spectroscopic binary with two moving components (SB2; Dufton et al. 2019), but it has never been analysed consistently with a stellar atmosphere code.

Fig. 1. Optical (left) and infrared (right) image of NGC346 in the SMC. Each image shows an area on the sky of approximately 3.5′×3.5′ (≈62 pc × 62 pc at d = 61 kpc). The positions of the main ionising sources of NGC 346 (see Sect. 6.4), namely SSN 7 (ON3 If*+O5.5 V((f))) and SSN 9 (O2 III(f*)), are highlighted by red and green squares, respectively. The false-colour optical image was taken with the HST ACS. Credit: A. Nota (ESA/STScI). The false-colour infrared image was taken with the James Webb Space Telescopes (JWST) Near Infrared Camera (NIRCam). Credits: NASA, ESA, CSA, O. Jones (UK ATC), G. De Marchi (ESTEC), and M. Meixner (USRA). Image processing: A. Pagan (STScI), N. Habel (USRA), L. Lenkic (USRA), and L. Chu (NASA/Ames).

The archival data, including spectra that cover the full range of the electromagnetic spectrum as well as photometric observations, are detailed in Sect. 2. In Sect. 3 we explain how the spectroscopic data is used to obtain radial velocity (RV) shifts, used to conduct an orbital analysis. Furthermore, this data is used in Sect. 4 to make a consistent spectroscopic analysis, yielding stellar and wind parameters. The empirically derived stellar and orbital parameters are put into an evolutionary context in Sect. 5 and its implications are discussed in Sect. 6.

2. Observations

NGC 346 SSN 7 has been observed numerous times across almost the entirety of the electromagnetic spectrum, ranging from the far-ultraviolet (FUV) to the infrared (IR). In this work we complement UV observations from the Hubble Space Telescopes (HST) UV Space Telescopes Imaging Spectrograph (STIS) and multi-epoch datacube observations obtained by the Multi Unit Spectroscopic Explorer (MUSE), which is mounted on the Very Large Telescope (VLT), with archival spectra from X-shooter, GIRAFFE, the UV-visual Echelle Spectrograph (UVES), and far-ultraviolet (FUV) spectra from the Far Ultraviolet Spectroscopic Explorer (FUSE). A summary of all spectra employed in this work is given in Table 1, where each spectrum is assigned an ID number.

The UV (including the FUV) spectra are rectified by division through the synthetic model continuum, while the optical and near-IR spectra are rectified by hand. The photometric data used for the fitting of the spectral energy distribution (SED) is given in Table 2.

The core of NGC 346, including our target and other bright O stars (namely MPG 396 (O7V), MPG 417 (O7.5V), MPG 470 (O8.5III), and MPC 476 (O6V)), was observed with the Chandra X-ray telescope (Weisskopf et al. 2000). The core region has a combined X-ray flux, already corrected for reddening, of F_X,  obs = 1.5 × 10³⁴ erg s⁻¹ cm⁻² in the 0.3 keV–10.0 keV (Nazé et al. 2002). We use this as an upper constraint on the X-ray flux that we can include in our atmospheric models. However, we are cautious in our use of this value, as close O-type binaries in the Milky Way suggest that the majority of the X-rays might originate from colliding wind zones of close binaries (Rauw & Nazé 2016).

3. Orbital analysis

NGC 346 SSN 7 is a known SB2 binary; both components show contributions in the He I and He II lines. The ephemeris of this system was obtained by Niemela & Gamen (2004) who obtained RVs by fitting helium lines, yielding an orbital period of 24.2 d and an eccentricity of 0.42. The high eccentricity is surprising for an O-star system with an orbital period below 50 d, as theory predicts that the orbit in such systems is quickly circularised. Unfortunately, we do not have access to their data. Still, we have 17 optical spectra, of which 11 (the MUSE spectra) cover half of the reported period, allowing us to make further restrictions to the estimate of the orbital period and the eccentricity, and to constrain a mass ratio.

3.1. Line associations to individual binary components

To prepare for the extraction of the RVs of each component, we identified which lines originate from which star. The most prominent line, which is detected in almost all of the used optical spectra of SSN 7, is the He IIλ 4686 line. In Fig. 2 we illustrate how the morphology of the line changes within the individual observations. For clarity, only a handful of observations are shown. The emission feature of He IIλ 4686 line shifts opposite to the absorption feature across multiple epochs. All lines that shift with the emission feature of the He IIλ 4686 line are associated with the spectroscopic primary, and all lines that shift with the He IIλ 4686 absorption line are associated with the spectroscopic secondary.

Fig. 2. Multi-epoch spectral observations of SSN 7. The X-shooter spectrum (ID 09) is shown as a solid line, the GIRAFFE spectrum (ID 17) as a dashed line, the GIRAFFE spectrum (ID 18) as a dot-dashed line, and the UVES spectrum (ID 19) as a dotted line. The shifting emission and absorption features are associated with the primary and secondary, respectively.

A full table of the considered lines is presented in Table A.3, including which line can be assigned to the primary, secondary, or to both. Spectral lines that only show a contribution of either the primary or the secondary are scarce.

3.2. Radial velocities

Measuring RVs in a binary, where both stars contribute to most of the diagnostic lines, is a non-trivial task. We employed a Markov chain Monte Carlo (MCMC) method combined with a least-squares fitting method to derive the RVs of each component. In each MCMC step, the individual synthetic spectra, which were obtained from the spectroscopic analysis (see Sect. 4), are shifted by a different RV value until the combined synthetic spectrum matches the observations. Thus, we can obtain RVs from complex absorption profiles. For a more detailed explanation we refer to Pauli et al. (2022b, their Sect. 3.1.2 and Appendix A).

The MCMC method returns a final probability distribution around the true solution. The errors quoted are the larger error margin of the 68% confidence interval of the final distribution: commonly used codes for inferring binary parameters from observables are not able to include asymmetric uncertainties. To minimise uncertainties introduced by the intrinsic variability of the lines, we fitted several lines that show distinguishable contributions of the primary and secondary components and calculated their mean RV.

The mean RVs obtained for each spectrum, and binary component are listed in Table 3. The fitted RVs of all lines used in this process are listed in Tables A.1 and A.2. The change in sign of the RV of the primary (from ∼100 km s⁻¹ to ∼−80 km s⁻¹) within ∼2 d calls in to question the reported period of 24.2 d and limits the orbital period to a few days.

3.3. Radial velocity curve modelling

When fitting RVs obtained from spectroscopic data with observations separated by several years, it is possible to find multiple orbital solutions consistent with the observations. To obtain a constraint on alternative ephemeris, a Monte Carlo sampler suitable for two-body systems with sparse and/or noisy radial velocity measurements, namely The Joker (Price-Whelan et al. 2017), is applied to the measured RVs. As The Joker is only able to model the RV curve of one component at a time, its output is used in a second step as input parameters for the Physics of Eclipsing Binaries (PHEOBE) code (Prša & Zwitter 2005; Prša et al. 2016; Horvat et al. 2018; Jones et al. 2020; Conroy et al. 2020), which can model the orbits of both binary components simultaneously.

3.3.1. Obtaining ephemerides with The Joker

The Joker uses a rejection sampling analysis on a densely sampled prior probability density function (pdf) and produces a sample of multimodal pdfs which contain the most important information about the ephemerides. The default pdfs of The Joker code (Price-Whelan et al. 2017, their Sect. 2) are used here as the priors. As outlined in Sect. 3.2, the RVs obtained from the MUSE spectra suggest that the orbital period should be of the order of a few days (see Table 3). Hence, a log-uniform prior for the orbital period in the range of 0.1 d–40 d is adopted for the fitting procedure. The upper limit was chosen as roughly two times the orbital period obtained by Niemela & Gamen (2004) to enable their solution to be considered, while the lower limit was set to a reasonably small value. For the eccentricity we used the β distribution from Kipping (2013) with a = 0.867 and b = 3.03. For the argument of the pericentre and the orbital phase a uniform prior without any restrictions is used. The RVs of the primary, which are more reliable since they have lower error margins, are employed as input data.

From the 1 × 10⁷ samples created, 340 pass the rejection step, equivalent to an acceptance rate of 0.0034%. Afterwards, a standard MCMC method is used to sample around the most likely solution found from the rejection step in order to fully explore the posterior pdf. We find that the ephemerides are described best with a circular orbit (e = 0) and an orbital period of P_Joker = 3.07438 ± 0.00002 d. The complete set of emphemerides obtained with The Joker are listed in Table 4.

3.3.2. Obtaining orbital parameters with PHOEBE

The ephermerides obtained in Sect. 3.3.1 were determined by only using the RVs of the primary. This was done to get a first constraint on the orbital parameters; however, valuable information such as the mass ratio (q), the projected orbital separation (a sin i), and the projected Keplerian masses (M sin i) can only be obtained when fitting the RVs of the two components simultaneously. Therefore, the PHOEBE light and RV curve modelling software (version 2.3) was employed to model the observed RVs, with the built-in option of the MCMC sampler +emcee+ (Foreman-Mackey et al. 2013).

The initial probability distributions of the period P, the argument and phase of the pericentre ω and M₀, and the barycentre velocity v₀ are assumed to be of Gaussian shape. The Gaussians are centred at the best-fitting value as obtained by The Joker (see Table 4) and have a standard deviation similar to their respective largest error margin. The ephemerides obtained with The Joker favour a circular orbit, and hence the eccentricity is fixed at e = 0. To constrain the current mass ratio we employed the method of Wilson (1941) yielding q_wilson = 1.34 ± 0.07. Because of a lack of information about other parameters, we used a uniform distribution in the range of q_orb,  i = 0.33 − 3.0, where a starting value q_orb,  i = q_wilson is assumed. The initial projected orbital separation is approximated to be of the order of a sin i_i = 10 R_⊙ and that it is uniformly distributed over the range of a sin i_i = 0.1 R_⊙ − −100 R_⊙.

The best fit of the RV curve obtained with the PHOEBE code is shown in Fig. 3. The final orbital parameters are listed in Tables 5 and 6. The RV fit yields masses without information about the inclination. By using the results from the spectroscopic analysis (see Sect. 4) we estimate the inclination to be i_orb = 16 ± 1°.

Fig. 3. Observed (triangles) and synthetic (solid lines) RV curves of the primary (green) and secondary (red) component, as obtained by the PHOEBE code. The dashed black line indicates the barycentric velocity offset v0.

To date, there is no observed light curve of NGC 346 SSN 7 reported in the literature. Since the PHOEBE code is capable of modelling a light curve, it can be used to make a prediction on the shape of a possible light curve (see Appendix B for more details). From the predicted synthetic light curve, we conclude that a detector with a precision below < 18 mmag in the V band is needed to see any kind of variability. This precision is hardly achievable for stars in the SMC with current ground-based telescopes. Hence, it is not surprising that no light curve has been detected so far.

4. Spectral analysis

4.1. Stellar atmosphere modelling

Expanding stellar atmospheres of hot stars (T > 15 kK) are in non-local thermal equilibrium (non-LTE) and can only be calculated using state of the art stellar atmosphere codes such as the Potsdam Wolf-Rayet (PoWR; Gräfener et al. 2002; Hamann & Gräfener 2004; Oskinova et al. 2011; Hainich et al. 2014, 2015; Shenar et al. 2015; Sander et al. 2015; Hainich et al. 2015) code. Within PoWR, the radiative transfer equations and the rate equations are solved iteratively in the co-moving frame. The code is able to produce a synthetic spectrum that accounts for mass-loss, line blanketing, and wind clumping.

For our stellar atmosphere models, we used tailored surface abundances that are representative of the chemical composition of the SMC. For the modelling of SSN 7 within NGC 346, the initial H abundance is adopted from Asplund et al. (2005), while the initial abundances of C, N, O, Mg, and Si are based on the works of Hunter et al. (2007) and Trundle et al. (2007). The abundances of P and S are taken from Scott et al. (2015) and are scaled down to the metallicity of the SMC (Z_SMC = 1/7 Z_⊙). The iron group elements (Fe, Sc, Ti, V, Cr, Mn, Co, and Ni) are either taken from Trundle et al. (2007) or, when not available, taken from solar values (Scott et al. 2015) and multiplied by a factor 1/7 to account for the lower metallicity of the SMC (Table 7). The numerous complex iron lines are handled with a ‘superlevel’ approach (Gräfener et al. 2002).

Microturbulence within a stellar atmosphere leads to a Doppler-broadening of spectral lines. The PoWR models used assume a depth-dependent microturbulence starting in the photosphere at ξ_ph = 10 km s⁻¹ and increases outwards linearly with the wind velocity ξ(r) = 0.1 ⋅ v(r).

Our PoWR models account for wind inhomogeneities and introduced optically thin clumps (microclumping) within the stellar atmosphere. The density contrast of the clumps within the wind is described by the clumping factor D. For NGC 346 SSN 7, a depth-dependent clumping factor is used, starting at the sonic radius and reaching D = 10 at ten times the stellar radius.

For the velocity stratification within the wind, a double β-law (Hillier & Miller 1999) of the form

$$\begin{array}{c}\hfill v(r)=v_{\infty }[f{(1-\frac{r_0}{r})}^{{\beta }_1}+(1-f){(1-\frac{r_1}{r})}^{{\beta }_2}]\end{array}$$(1)

is adopted. The inner part of the velocity stratification is described by an exponent β₁ = 0.8, which is typical for O stars (Pauldrach et al. 1986). The outer part is described with an exponent β₂ = 4 and only contributes with a fraction of f = 0.4. The parameters r₀ and r₁ are set within the code, such that the wind is smoothly connected.

The effective temperatures quoted in this work are related via the Stefan-Boltzmann law to the radius at which the Rosseland mean continuum optical depth is τ = 2/3. Stellar and wind parameters can be obtained by adjusting the synthetic spectrum until it matches the observed spectra simultaneously. This also means that the abundance of a specific element (e.g. H, C, N, O) can be adjusted if needed to match the strength of observed lines.

4.2. Resulting spectroscopic parameters

Deriving stellar and wind parameters from stellar spectra in which lines are blended by two stars of a binary is a difficult task. The following subsections describe the strategy employed and the final parameters that are listed in Table 8.

4.2.1. Rotation rates

Before starting the spectral fitting process, accurate values for the projected rotational velocities (v sin i) of each component are needed. These impact the depth and shape of all synthetic lines and hence conclusions drawn from the spectroscopic fit. We employ the IACOB-BROAD tool (Simón-Díaz & Herrero 2014) which uses a combination of a Fourier transform (FT) and a goodness of fit (GOF) method to obtain accurate measurements of v sin i.

When obtaining projected rotation rates, the preference is to use metal lines over helium lines as the latter are, in addition to the rotational broadening, also pressure broadened. For the primary, we use the N IVλ 4057 as well as the N IVλ 6381 line. As there are no metal lines that solely can be attributed to the secondary, we use the prominent He Iλ 6875 and the He Iλ 7065 line. The projected rotational velocities of the primary and secondary component are best described with v₁ sin i = 135 ± 10 km s⁻¹ and v₂ sin i = 185 ± 10 km s⁻¹, respectively.

Most of the observed lines are formed in the rotating photosphere of a star, while fewer but still important lines in the UV are formed within the static (non-rotating) wind. To account for this, the synthetic PoWR models used in this work are calculated with a rigidly rotating photosphere and a non-rotating wind, as described in Shenar et al. (2014).

4.2.2. Luminosity, temperature, and surface gravity

To constrain the luminosity of each binary component, we carefully looked at the ratios of isolated and blended helium and metal lines. One of the key diagnostic lines was the C IIIλ 1175 in the UV, which originates purely from the secondary, giving an extra constraint on the luminosity ratio. A selected set of the lines used is illustrated in Fig. 4. Determining the luminosity ratio is an iterative process and is revisited after each modelling change, as changes in temperature and surface gravity also affect the shape of key diagnostic lines. The luminosity of each component is log(L₁ / L_⊙) = 5.75 and log(L₁ / L_⊙) = 5.78 for the primary and secondary, respectively. The synthetic SED of each component and the combined synthetic SED compared to the observed spectra and photometry is shown in Fig. 5.

Fig. 4. Close-ups on the carbon and helium lines used to constrain the luminosity ratio (solid blue) compared to the combined synthetic model (dotted red). Left: observed FUSE spectra (ID 01). The weighted synthetic spectra of the primary and secondary components are shown as a dotted black and green line. Right: observed X-shooter spectra (ID 09 and 31). The unweighted synthetic spectra of the primary and secondary components are shown as a dotted black and green line.

Fig. 5. Observed (blue) and combined synthetic (red) SED. The individual synthetic SEDs of the primary and secondary are shown as the dotted black and green line, respectively. The photometry, listed in Table 2, is shown as open squares.

The surface gravity of each of the component stars is determined by fitting the wings of the Balmer lines in the spectra with the biggest RV shifts. We find that the primary and the secondary have similar surface gravities of log g = 3.7 ± 0.1 (see Fig. 6).

Fig. 6. Close-ups on the Balmer lines in the X-shooter spectrum (ID 09; solid blue) compared to the combined synthetic model (dotted red). The unweighted synthetic spectra of the primary and secondary components are shown as dotted black and green line.

The temperatures of the stars are derived using the ratio between the He I and He II lines associated with each component (see Table A.3). To get a more accurate constraint on the primary’s temperature, the ratio between the different observed nitrogen lines are used, namely N IVλλλ 3479, 3483, 3485, N IVλ 4058, N Vλλ 4604, 4620, N IIIλλ 4634, 4641, and N IVλλλ 6212, 6216, 6220. We find for the primary T_eff,  1 = 43.6 ± 2 kK, while the secondary is best fit with T_eff,  2 = 38.7 ± 2 kK.

4.2.3. Surface abundances

During the determination of the temperature of the primary, we find that it is possible to match the ratio of the He I and He II lines, but these lines are systematically too weak. Hence, we lowered the primary’s surface hydrogen abundance to X_H,  1 = 0.60 ± 0.1. Furthermore, to fit the observed nitrogen lines, which are also crucial for the temperature determination, we find that the nitrogen abundance of the primary is enhanced to $X_{\mathrm{N},1}={138}_{-30}^{+10}\times {10}^{-5}$ (i.e. 2 N_⊙).

In addition, the spectrum contains carbon and oxygen lines that show clear contributions from the primary. To constrain its surface oxygen abundance, we used the O IVλλλ 1337-1343-1344 triplet in the UV, the O IV multiplet around ∼3000 Å, and the O III lines in the range from 1409 Å–1412 Å. We find that oxygen in the primary’s atmosphere is depleted with $X_{\mathrm{O},1}=4_{-3}^{+6}\times {10}^{-5}$. Regarding the carbon abundance, we used C IVλ 1169 in the UV and the C IVλλ 5801, 5812 doublet in the optical, yielding a surface carbon abundance of $X_{\mathrm{C},1}=1_{-0}^{+9}\times {10}^{-5}$ in the primary.

Determining the surface abundance of the secondary is more complicated as it shows little contribution to the observed metal lines. We can see in some spectra that the secondary contributes to the N IIIλλ 4634, 4641 doublet. Hence, we tried to use this as a constraint on its surface nitrogen abundance and find the best fit with $X_{\mathrm{N},2}={60}_{-10}^{+20}\times {10}^{-5}$. We can also see a contribution to the C IVλλ 5801, 5812 doublet in the optical, and the C IIIλ 1175 in the UV. These lines can be well reproduced with the initial carbon abundance of $X_{\mathrm{C},2}={21}_{-5}^{+0}\times {10}^{-5}$. Unfortunately, we cannot see a clear contribution of the secondary to the oxygen lines. As the nitrogen abundance is increased and theory predicts that the amount of CNO material should be X_CNO ≈ 137 × 10⁻⁵ for stars in the SMC (Kurt & Dufour 1998), we scaled down the oxygen abundance of the secondary to X_O,  2 = 80 × 10⁻⁵. We note that this does not have any impact on the quality of our fit.

Figure 7 shows the all-important metal lines used to fit the abundances of the primary and secondary. The changed surface abundances of the primary, especially the lower surface hydrogen abundance and the increased nitrogen, hint at previous binary interaction that stripped off parts of the envelope (see Sect. 6). These metal lines are also impacted by microturbulence (see Sect. 4.1), Fig. 8 shows the adopted microturbulence values for each component result in a composite spectrum consistent with the UV observations.

Fig. 7. Close-ups on the metal lines in the X-shooter spectrum (ID 09 and 31; solid blue) compared to the combined synthetic spectrum (dashed red). The unweighted synthetic spectra of the primary and the secondary component are shown as a dotted black and green line, respectively.

Fig. 8. Close-up on the iron forest of the HST FOS (ID 03) spectrum (solid blue) compared to the combined synthetic model flux (dashed red). The weighted synthetic model spectra of the primary and secondary are shown as a dotted black and green line, respectively. The iron forest is used to check the iron abundance and to constrain the microturbulence within the stellar atmosphere.

4.2.4. Spectroscopic stellar masses

Given the fundamental stellar parameters (T_eff, log g, L) of each binary component, we can calculate the spectroscopic masses. We find that the mass of the primary component is $M_{\mathrm{spec},1}={32}_{-13}^{+19}{M}_{\odot }$, surprisingly low for its luminosity. The mass of the secondary is $M_{\mathrm{spec},2}={55}_{-21}^{+32}{M}_{\odot }$. The low mass of the primary strengthens our argument that the binary components must have interacted in the past.

To further strengthen this hypothesis, we calculated the Roche radius, using the results from the orbital analysis in combination with the spectroscopic masses and find that both stars are currently slightly overfilling their Roche lobe with R₁/R_RL = 1.01 and R₂/R_RL = 1.03. This strongly suggests that the system interacted in the past, and that both stars are still currently in contact.

4.2.5. Wind velocities and mass-loss rates

So far only photospheric lines were considered in the spectroscopic analysis. However, in the spectra there are several lines that are sensitive to wind parameters, like the N Vλλ 1239, 1243 and C IVλλ 1548, 1551 doublets in the UV, and Hα and He IIλ 4686 in the optical. Both stars are very bright, massive, and hot, meaning that both of them are expected to contribute to the wind diagnostic lines. Within the FUV FUSE observation, no measurable P V or S IV wind resonance lines are present. This is consistent with our composite model.

The primary motion is reflected by the emission parts of Hα and He IIλ 4686, while the secondary’s motion is evident from the absorption components of the same lines. We conclude that the primary must have a stronger mass-loss rate compared to the secondary. We used the Hα and He IIλ 4686, in combination with N Vλλ 1239, 1243 and C IVλλ 1548, 1551, to calibrate the primary’s mass-loss rate to log(Ṁ₁/(M_⊙ yr⁻¹)) = −5.4 ± 0.1. With this mass-loss rate, the resonance doublets in the UV were fully saturated in the synthetic spectrum of the primary, yet the observed profiles were not fully saturated. The secondary contributes roughly 60% of the UV flux, so we fit the combined spectra with a partially saturated secondary wind. The secondary is fit with a mass-loss rate of log(Ṁ₂/(M_⊙ yr⁻¹)) = −7.3 ± 0.1, which results in the combined synthetic spectrum matching the observed spectrum (Fig. 9).

Fig. 9. Close-ups on the key wind diagnostic lines (solid blue) in the FUV, UV and optical spectra overplotted by the combined synthetic model (dashed red). In the upper panels the weighted and in the lower panels the unweighted synthetic spectra of the primary and secondary components are shown as dotted black and green line, respectively. Upper left: the O VI line of the FUSE (ID 1) spectrum. This line is sensitive to X-rays. Upper middle: the N V resonance doublet in the HST FOS (ID 3) spectrum. Upper right: the C IV resonance doublet in the HST FOS (ID 3) spectrum. Lower left: the He IIλ 4686 line in the X-shooter spectrum (ID 9), showing strong wind emission features from the primary. Lower right: Hα of the X-shooter spectrum (ID 31), again, the primary dominates the wind emission of this line.

The terminal wind velocity of the primary is adjusted to match the blue edge of the N Vλλ 1239, 1243 absorption trough, giving v_∞,1 = 2500 ± 100 km s⁻¹. For the secondary, we experimented with different values of v_∞ to see how the slope of the C IVλ 1548,1551 resonance doublet changes. The best fit that matches the profile in the STIS and FOS spectrum was obtained when using the same terminal wind velocity for the secondary as for the primary (v_∞,2 = 2500 ± 300 km s⁻¹).

4.2.6. Required X-ray flux

The O VIλλ 1032, 1038 resonance doublet in the FUV shows weak indications of a P Cygni profile, which cannot be modelled by default in a stellar atmosphere model. Cassinelli & Olson (1979) were the first to suggest that these high ionisation stages must be ‘super-ionised’ by a hot X-ray plasma that originates from within the wind.

We include, in our primary and secondary, a hot plasma that is smoothly distributed throughout the wind, starting at 1.1 R_* and having a temperature of T_X = 3 MK. We find that the O VIλλ 1032, 1038 resonance doublet (Fig. 9) can be explained with an X-ray luminosity of log(L_X,  1/(erg s⁻¹)) = 32.4 and log(L_X,  2/(erg s⁻¹))  =  32.1 for the primary and secondary, respectively. This is much lower than the X-Ray luminosity observed from the core of NGC 346 and comparable to X-rays needed for other O stars (Nebot Gómez-Morán & Oskinova 2018).

5. Stellar evolution modelling

In the previous sections we determine the stellar and wind parameters of the binary components, and refine the binary period, and we establish that the system is most likely in contact. With this information at hand, we proceed to model the past, present, and future evolution of this remarkable system.

The evolution of SSN 7 in this work is modelled using the Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2011, 2013, 2015, 2018, 2019) code (v.15140). Our aim is to understand, from the predictions of the stellar evolutionary models, the current evolutionary stage of our target rather than finding a perfect fine-tuned fit for our empirically derived parameters. The orbital and spectroscopic parameters (Sects. 3.1 and 4.1) limit the possible previous evolutionary paths. To find a fitting model, we explore a parameter space in the following ranges: donor masses in the range M_1,  i = 45M_⊙ − 65M_⊙, initial accretor masses with M_2,  i = 30M_⊙ − 55M_⊙, and initial orbital periods in the range P_i = 2 d − 7 d.

5.1. Stellar and binary input physics

The employed binary evolutionary models are constructed using the input files provided by Marchant et al. (2016)¹. A brief summary of the assumed input physics follows.

In the binary models, following the work of Brott et al. (2011), tailored initial abundances are used. The initial abundances are X_H = 0.7460, X_He = 0.2518, and Z = 0.0022, with Z being the total metal fraction. The individual initial metal abundances are listed in Table 9. The abundances are comparable to fit in our atmospheric model (see Table 7).

Mixing is included as follows. Convection is modelled as standard mixing length theory (MLT; Böhm-Vitense 1958) using the Ledoux criterion and α_mlt = 1.5. Semiconvective mixing is assumed to be efficient with α_sc = 1 (Langer et al. 1983; Schootemeijer et al. 2019). For core hydrogen-burning stars, the step overshooting is assumed to be up to 0.335H_P (Brott et al. 2011; Schootemeijer et al. 2019). Furthermore, rotational mixing is included as a diffuse process, including dynamical shear, secular shear, and Goldreich-Schubert-Fricke instabilities, as well as Eddingtion-Sweet circulations (Heger et al. 2000). The efficiency parameters are set to f_c = 1/30 and f_μ = 0.1. Thermohaline mixing is modelled with α_th = 1 (Kippenhahn et al. 1980). Lastly, the models take into account momentum transport originating from a magnetic field as expected from a Taylor-Spruit dynamo (Spruit 2002).

To avoid numerical complications in the late evolutionary stages of the star, the evolution of both components is only modelled until core helium depletion. Furthermore, for core helium burning stars with core masses above > 14 M_⊙ a more efficient mixing theory, MLT++ (Sect. 7.2 of Paxton et al. 2013), is employed.

Mass-loss via a stellar wind is included in the models in the following way. For hydrogen-rich (X_H ≥ 0.7) OB stars, the recipe of Vink et al. (2001) is used. As soon as the temperature of any model drops below the bi-stability jump (Vink et al. 2001, their Eqs. (14) and (15)), the maximum of the mass-loss rate of Vink et al. (2001) or Nieuwenhuijzen & de Jager (1990) is taken. During the WR stage, when the surface hydrogen abundance drops below X_H ≤ 0.4, the mass loss rates of Shenar et al. (2019, 2020b) are used. For surface abundances between 0.7 > X_H > 0.4, the mass-loss rate is linearly interpolated between the mass-loss rates of Vink et al. (2001) and Shenar et al. (2019, 2020b).

To reduce the free parameter space, it is presumed that the components of SSN 7 are initially tidally locked and that the orbit is circularised. This is a reasonable approach for a binary with a period as short as 3.1 d.

Mass transfer by Roche-lobe overflow is modelled implicitly using the contact scheme of the MESA code. The mass transfer is modelled conservatively as long as the mass-gainer is below critical rotation. Whenever the mass-gainer spins up to critical rotation, it will stop the accretion, and we assume that the material is directly lost from the system. In addition to the mass transfer by Roche-lobe overflow, mass accretion from the stellar wind is taken into account.

At the instance of time when the donor star depletes helium in its core, its model is stopped. It is assumed that the primary directly collapses into a BH without a supernova explosion and without a kick. Further mass transfer via Roche-lobe overflow and wind from the companion onto the BH accretion onto the BH is limited by the Eddington accretion rate.

5.2. Resulting binary models

The best fit of the empirically derived fundamental stellar parameters of both binary components is achieved with an initial primary mass of 52 M_⊙, an initial secondary mass of 38 M_⊙, and an initial orbital period of 2.7 d. The corresponding tracks of the mass donor and mass gainer in the Hertzsprung–Russell diagram (HRD) are shown in Fig. 10. The binary evolutionary model that matches the spectral fit parameters best predicts that the observed stars are currently in contact in a slow Case A mass-transfer phase. This prediction would make SSN 7 the most massive Algol-like system known so far, making it an ideal laboratory to study the ongoing effects and efficiencies of mass transfer in binaries hosting massive stars.

Fig. 10. Evolutionary tracks of the donor (left) and the accretor (right) colour-coded by their surface H-abundance. The track of the mass gainer is only colour-coded until the primary dies and then continues as a solid white line, assuming that the binary is not disrupted by a supernova explosion. The tracks are overlayed by black dots that show equidistant time steps of 30 000 yr to highlight phases in which the stars spend most of their time. The spectroscopic results are shown as green and red triangles surrounded by error ellipses. Start and end phases of fast and slow Case A, as well as Case AB mass-transfer are indicated by arrows. The iso-contours of equal radii are indicated in the background as dashed grey lines.

The parameters of the closest fitting binary evolution model are listed in Table 8. Errors quoted are sophisticated guesses on the uncertainties of the binary evolutionary models from the experiences we gained when tailoring the evolutionary models to match the empirically derived stellar parameters. For more accurate error margins, a fine grid of detailed stellar evolution models is needed.

The empirically derived fundamental stellar parameters as well as the observed period of 3.1 d can be explained by our binary evolutionary model. Our favourite model predicts that the binary has an age of 4.2 Myr. The only discrepancy found is between the predicted and observed rotation rates. We discuss this in detail in Sect. 6.2.2.

It is predicted that the primary, after it depletes the hydrogen in its core, expands and initiates another more rapid mass-transfer phase, also called Case AB mass-transfer. During this phase it strips off large fractions of its H-rich envelope and evolves bluewards in the HRD towards the regime where WR stars can be found. During its WR phase, our model is predicted to have some H remaining in the envelope (X_H ∼ 0.2). During this time, the secondary, which will be rejuvenated from mass accretion, is still core hydrogen burning.

After the primary, which is the current mass donor, depletes the helium in its core, we assume that it directly collapses into a BH. The formed BH has a mass of 25.5 M_⊙, which should be considered as an upper limit since we assumed a direct collapse. However, this assumption allows us to make predictions about the future evolution of the secondary. Our model predicts that the secondary will evolve off the main sequence and will initiate mass transfer on the BH, stripping parts of its H-rich envelope (X_H ∼ 0.4). This star does not evolve bluewards to the WR regime and the Of/WNh population. This is linked to the increased envelope-to-core mass ratio, leading to the formation of a steep chemical gradient (Schootemeijer & Langer 2018; Pauli et al. 2023). After a short time (of the order of a few hundred thousand years), the remaining star will also collapse to a BH with a mass of 25 M_⊙. The two black holes orbit each other every 4.85 d. Following Peters (1964) we calculated the merger timescale of the binary to be τ_merger = 18.5 Gyr. This is longer than the age of the Universe. We note that our estimated value should be considered a lower limit as the BH masses of the individual components could be lower, which would result in even longer merger timescales.

6. Discussion

6.1. Orbital, spectroscopic, and stellar evolutionary masses

The most reliable mass estimate is the projected orbital mass of the system (see Sect. 3). Since we do not have a light curve available, it is only possible to compare it with the spectroscopic result by choosing an inclination. However, this introduces additional uncertainties and makes the orbital masses inaccurate. A better way to compare the results is by using the derived mass ratios.

In Sect. 3 we obtained q_wilson = 1.34 ± 0.07 using the fitting method from Wilson (1941) and $q_{\mathrm{PHOEBE}}=1.{47}_{-0.09}^{+0.13}$ when modelling the RV curve in detail with the PHOEBE code. From our spectroscopic analysis, we derived a somewhat higher mass ratio of $q_{\mathrm{spec}}=1.7_{-1.0}^{+2.9}$. The large error margins arise from the uncertainties from fitting the surface gravity and from the uncertainties introduced in the multiple solutions that can be found for the luminosity ratio (i.e. M_spec ∝ R² ∝ L). Hence, it could be that the spectroscopic mass ratio is overestimated. However, the conclusion that the primary component is more massive than the secondary component remains.

During the tailoring process of our evolutionary model towards the empirically derived stellar parameters, we exploited various combinations of the mass ratio of the two components. In these models, the mass and the position of the primary were always only explainable when it was in a slow Case A mass-transfer phase. The secondary in these models had different initial masses, and thus the mass during the Case A mass transfer was also different (i.e. lower) compared to our best-fitting model.

6.2. Effects of the inclination on the final parameters

6.2.1. Fundamental stellar and wind parameters

As discussed in the previous section, the orbital mass ratio and the spectroscopic mass ratio agree within their respective uncertainties. This agreement enabled us to estimate the inclination of the system to be i_orb = 16 ± 1° (see Sect. 3.3.2). This implies that we must be looking at the system pole-on. Given our estimates on the rotational velocities and the fact that the two stars are in contact, they must be oblate and deformed. Hence, the stellar parameters of these stars might no longer be uniform over the surface of the star.

Under the assumption of a star rotating close to break-up velocity, the effective surface gravity at the equator is g_equator ≈ 0 as the gravitational and centrifugal forces cancel each other out. On the other hand, the effective surface gravity on the pole must be rather unperturbed and is mostly affected by gravity. Since we are looking at the pole region of the star, we conclude that our estimate on the surface gravity is reliable within its uncertainties, and hence also our estimate on the inclination.

Due to the effect of gravity darkening, an ellipsoidal star is cooler at the equator and hotter at the poles. Abdul-Masih (2023) studied the effect of rotation and inclination on the determined effective temperature by using synthetic models created with the SPAMMS code on a three-dimensional surface. He concluded that in the case of rapid rotation and low inclination, the star’s temperature can be overestimated by 10%. This means that both of our stars actually might be 3 kK–4 kK cooler. Furthermore, Abdul-Masih (2023) finds that the surface helium abundance of a rapidly rotating pole-on star can be underestimated by as much as 60%, which would account for the discrepancy between the observed helium abundance and the helium abundance from the binary evolutionary models.

Given the ellipsoidal shape and hence the different radii and temperatures, it might be evident that the star’s luminosity also depends on the inclination angle. We calculated the expected luminosity of our target at the pole and the equator and find a difference Δlog L/L_⊙ ≈ 0.15, which is within our given uncertainties. The impact of a somewhat lower luminosity on our final results is too small to impact our final conclusions.

Due to the different effective surface gravities between the poles and the equator, wind parameters vary with inclination. The wind is slower and weaker in the equatorial regions, while it is faster and more powerful at the poles. An average mass-loss rate for the whole star would thus be lower than any measured from a low inclination angle, such as in our case. We would therefore treat our stated mass-loss rates as upper limits.

6.2.2. Discrepancy between observed and predicted rotation velocity

The binary evolutionary models agree with almost all spectroscopically derived stellar parameters, except the rotational velocity. Under the assumption that the inclination determined from the orbital analysis (see Sect. 3) is reliable, the observed rotation velocities of the binary components (v₁ = 490 ± 47km s⁻¹ and v₂ = 544 ± 49km s⁻¹) are a factor of two higher than those predicted by the binary evolutionary models.

Given that we are looking at the binary pole-on, and that both stars are deformed due to their fast rotation, this means that the assumption of spherical symmetry breaks down in the evolutionary models. Zahn et al. (2010) have shown this for a uniformly rapidly rotating star R_equator ≈ 1.5R_pole. Hence, under the assumption that R ≈ R_pole, the rotation rates of the binary model need to be corrected to v_{mod,   1} ≈ 325 ± 45km s⁻¹ and v_{mod,   2} ≈ 425 ± 50km s⁻¹. This brings observations and predictions closer together, but still shows some discrepancy. We note that our target is a contact binary, meaning that the stars are tear-drop shaped rather than ellipsoidal. This makes an estimate of the rotation rates even more complicated.

On the other hand, it is worth mentioning that in our stellar atmosphere calculations we neglect additional broadening mechanisms such as macroturbulence and asynchronous rotation (e.g. introduced by differential rotation). These effects are compensated for in our estimate of the projected rotational velocity v sin i, and hence might lead to an overestimation, which can explain the observed discrepancy.

6.3. Empirical mass-loss rates

From our spectroscopic analysis (Sect. 4.2.5) we derived that the primary dominates most of the wind lines (log(Ṁ₁/(M_⊙ yr⁻¹)) = −5.4), while the secondary still has a significant and noticeable contribution in the wind UV lines (log(Ṁ₂/(M_⊙ yr⁻¹)) = −7.3). In the following we want to determine whether the chosen mass-loss rates agree with observations and theoretical expectations.

Following Rickard et al. (2022), we would expect from a star that has evolved in isolation and has a current luminosity of log(L/L_⊙) = 5.75 a mass-loss rate of only log(Ṁ_Rickard/(M_⊙ yr⁻¹)) = −7.35. This is two orders of magnitude lower than that derived for the primary from empirically spectroscopic analysis. However, the primary has lost a huge amount of mass and has an unusual L/M ratio that is not considered in their relation. On the other hand, theoretical mass-loss recipes take into account several parameters, such as luminosity, mass, temperature and metallicity. For comparison, we calculated the mass-loss rate predicted with the frequently used Vink et al. (2001) to be log(Ṁ_Vink,1/(M_⊙ yr⁻¹)) = −5.9. A more updated recipe for O-type stars is from Björklund et al. (2023). Their mass-loss recipe predicts log(Ṁ_{Bjrklund 2022, 1}/(M_⊙ yr⁻¹)) = −6.6, which is one order of magnitude lower than the observed mass-loss rate. However, these mass-loss rates are calculated for H-rich (X_H > 0.7). The primary has already lost a significant amount of its envelope, and its surface H-abundance already dropped to X_H ≈ 0.6 (see Table 8). For completeness, we calculated the corresponding WR mass-loss rates according to Shenar et al. (2019, 2020b) for stars with X_H > 0.4, yielding log(Ṁ_Shenar, 1/(M_⊙ yr⁻¹)) = −5.5 and being in agreement with the empirically derived values.

The secondary seems to be more well behaved. Its empirically derived mass-loss rates agrees well with the expectations from other O stars in the SMC (Rickard et al. 2022). The theoretical predicted rates from Vink et al. (2001) and Björklund et al. (2023) are log(Ṁ_Vink, 2/(M_⊙ yr⁻¹)) = −6.2 and log(Ṁ_{Bjrklund 2022, 1}/(M_⊙ yr⁻¹)) = −7.1, respectively. Only the latter agrees within the error margins of the observed mass-loss rates. For comparison, we illustrate in Fig. 11 the expected mass-loss rates from other O-type stars in NGC 346 compared to our spectroscopically fit values, alongside theoretical prescriptions.

Fig. 11. Observed mass-loss rates of O stars in NGC 346 (including SSN 7) as a function of stellar luminosity compared to theoretical and empirical relations. The position of both binary components of SSN 7 are marked by black stars. The data is complemented by O stars from Rickard et al. (2022). Apparently single stars are shown as filled blue squares and known muliple systems as open blue squares. In the background the empirical relation of main-sequence O-type stars from Rickard et al. (2022) is included as dashed black line. In addition, we show the theoretical relations of Vink et al. (2001) and Björklund et al. (2023) for main-sequence stars as dashed yellow and dashed red line, respectively.

6.4. SSN 7 and its role in NGC346

In the recent paper of Rickard et al. (2022), it is noted that the massive giant NGC 346 SSN 9 (MPG 355, O2III(f^*), Walborn et al. 2000) is the major ionising source of the core region. In their work it is stated that the star outshines the remaining O-star population and is responsible for 50% of the hydrogen ionising flux (log Q_{H,  SSN 9} = 49.98). However, in their analysis SSN 7 was neglected because of its binary nature. According to our stellar atmosphere models, we find that the two early-type stars in the binary have a combined hydrogen ionising flux of log Q_{H,  SSN7 a + b} = 49.88. This makes SSN 7 one of the two major contributors to the H ionising flux of the core of NGC 346. We stated in Rickard et al. (2022) that the morphology of the ionising front within images of NGC 346 centre upon SSN 9. The location of SSN 7 is directly ‘behind’ SSN 9 from the position of the ionising front (Fig. 1).

We compared the calculated H ionising photon count of both components of SSN 7 to other massive stars within both the SMC and LMC (Fig. 12). The predicted ionising flux is consistent with the spectral types assigned to each component by Walborn et al. (2000), even though the primary has an outstanding L/M ratio and the secondary that is suspected to have accreted several solar masses.

Fig. 12. log QH to spectral type. SMC sample stars analysed by Ramachandran et al. (2019) are included (red diamonds) as well as LMC stars from Ramachandran et al. (2018; blue diamonds) and O stars from the core of NGC 346 from Rickard et al. (2022; black diamonds). The two components of the SSN 7 binary are shown separately.

For helium ionising photons, the model of the primary gives log Q_{He I,  SSN 7a} = 48.98 and log Q_{HeII,  SSN 7a} = 40.81, while the secondary model gives log Q_{He I,  SSN 7b} = 48.58 and log Q_{HeII,  SSN 7b} = 43.14. Combined, the helium ionising fluxes of the two components of the binary combined are log Q_{HeI,  SSN7 a + b} = 49.12 and log Q_{HeI,  SSN7 a + b} = 43.14, or 25% and 0.04% of the total of the O stars analysed in the core of NGC 346 between this paper and Rickard et al. (2022). This is compared to log Q_{HeI,  SSN 9} = 49.41 and log Q _{HeII,  SSN 9} = 46.50, or 31% and 99%. This is notable because the combined He I Ionising flux for the binary is comparable to that of SSN 9, but not for the He II Ionising flux, where SSN 9 dominates.

6.5. Age of the NGC 346 O-star population

During our evolutionary analysis it was revealed that the system of SSN 7 has an age of 4.2 Myr (see Sect. 5.2). Meanwhile, much younger ages for NGC346 have been estimated in previous works, such as 1 − 2.6 Myr (Dufton et al. 2019) and 1 − 2 Myr (Walborn et al. 2000). In Fig. 13 we compare the empirically derived stellar parameters of SSN 7 and other O stars (Rickard et al. 2022) to isochrones from Georgy et al. (2013). According to these isochrones, the age of the primary should be 2 Myr, while the secondary would be 3 Myr old. This mismatch between the age predicted through evolutionary modelling and the age compared to single-star isochrones is not an isolated example and could lead to false conclusions on the age of the cluster.

Fig. 13. HR diagram showing the position of the two components of SSN 7 as determined by our spectrographic analysis. Black stars are the two SSN 7 components. Blue symbols are the other O stars within the central core of NGC 346 (Rickard et al. 2022). Open symbols are noted spectral binaries (Dufton et al. 2019). Dotted lines are isochrones with Z = 0.002 (Georgy et al. 2013) and an initial rotation of 100 km s−1.

When inspecting Fig. 13 in more detail we can see that most of the luminous O stars in NGC 346 are known binaries, or higher order multiple systems, with the one exception of SSN 9. However, this star is unusually luminous and hot and undermassive for its position in the HRD, hinting towards a binary nature. Taking into account only the apparent single O stars in this cluster, we conlcude the age of the O-star population, according to single-star evolution isochromes, to be 4 Myr, which is consistent with the age derived from our binary models.

With evidence that the core of NGC 346 has not yet been affected by any supernova explosion (Danforth et al. 2003), the fate of more massive stars within NGC 346 comes into question. Either none had formed yet, or those that did had imploded into BHs directly.

6.6. Luminosity class and spectral type

In the previous work of Walborn (1978), SSN 7 was classified as O4 III(n)(f), but it was not known to be a binary. In a recent work, Dufton et al. (2019) propose that the presence of He Iλ 4471 absorption and weak He IIλ 4686 emission suggest that the primary should be a O4 If star. With the new insights gained from the detailed multi-epoch observations, we assign the He Iλ 4471 absorption mainly to the secondary component, meaning that the primary has an earlier type.

Here we want to use the newly gained insights from the spectral analysis to reclassify the primary and secondary star. It is evident that the primary in our model shows only a very weak contribution to the He Iλ 4471 line (see Fig. 4), meaning that it must be of very early type (O2-3.5). Following Walborn et al. (2002) it is suggested to use for such early-type stars the ratio of N IVλ 4058 to N IIIλ 4640 emission. The multi-epoch observations have shown that the N IIIλ 4640 emission line is a blend of the primary and secondary components, while the N IVλ 4058 purely originates from the primary (see Fig. 7 and Table A.3). Guided by our stellar atmosphere models we find that N IVλ 4058 > N IIIλ 4640, while we see very weak or no absorption in He Iλ 4471, implying a spectral type of O3.

According to our spectral analysis we find that in the primary nitrogen is enhanced, while carbon and oxygen are deficient. Therefore, we add the qualifier ‘N’ to the spectral type. Following the classification criteria of Sota et al. (2011) and taking into account that the N IVλ 4058 emission feature is greater than the N IIIλ 4640 emission and that the He IIλ 4686 is also in emission we suggest adding the suffix ‘f^*’. In Sota et al. (2014) the luminosity classes are linked to the spectral type as a function of the f phenomenon, suggesting a luminosity class I for the primary. Because of the high temperature of the star, combined with the low Si IV abundance (∼2% solar), no Si IVλλ 1394, 1403 P Cygni lines are observed (Crowther et al. 2002), so the presence or absence of this line cannot contribute to our luminosity classification. The final spectral type we assign to the primary is ON3 If^*.

Observational counterparts with similar spectral type at low metallicity are sparse. Galactic examples are Cyg OB2-7 and Cyg OB2-22A. By comparison of their observed spectra to the synthetic spectrum of our primary we find a good agreement between the morphology of the lines, supporting our choice on the spectral type.

For the secondary, we limit its spectral type to O4-8 as it shows He Iλ 4471 in absorption, while He IIλ 4542  >  He Iλ 4388 and He IIλ 4200  >  He Iλ 4144, and it lacks Si IIIλ4552, which can be seen in the earliest O-type stars. Hence, we follow Conti & Alschuler (1971) and use the logarithmic ratio of the equivalent widths of He Iλ 4471 and He IIλ 4542. Unfortunately, the He IIλ 4542 line is blended by the companion. Therefore, we decided to use the equivalent widths of the synthetic spectrum of the secondary, yielding log(EW(4471)/EW(4542)) = −0.45. This corresponds to a spectral type of O5.5. As mentioned above, the N IIIλ 4640 emission line is a blend from the primary and secondary component. As the secondary does not contribute to the N IVλ 4058 line and shows strong He IIλ 4686 absorption we add the suffix ‘((f))’ (Sota et al. 2011). Following Sota et al. (2014), given the spectral type and the f phenomenon, the star should have a luminosity class V. We use the criterion for the luminosity class estimation introduced by Martins (2018). The criterion is based on the equivalent width of He IIλ 4686 absorption line in the synthetic model of the secondary. With an equivalent with of EW(4868) = 0.74, this star is confirmed to have a luminosity class V. The final spectroscopic classification of the binary is ON3 If^*+O5.5 V((f)).

7. Conclusions

We have conducted the first consistent multi-wavelength multi-epoch spectral analysis of the SB 2 contact binary SSN 7 located in the NGC 346 cluster in the SMC. The detailed analysis revealed that SSN 7 is a system with two binary components that have comparable luminosity (log(L/L_⊙)≈5.75), while having divergent stellar and wind parameters. The primary is the hotter (T_eff,  1 = 43.6 kK) but less massive (M₁ = 32 M_⊙) component. We revealed that the surface hydrogen abundance is X_H = 0.6, and that it has CNO values close to those of the CNO equilibrium. The mass-loss rate of the primary, log(Ṁ₁/(M_⊙ yr⁻¹)) = −5.4, is exceptionally high for its luminosity. The secondary is slightly cooler (T_eff,  2 = 38.7 kK), but more massive (M₂ = 55 M_⊙). The secondary has a mass-loss rate of log(Ṁ₂/(M_⊙ yr⁻¹)) = −7.3, which is in agreement with other O stars within the SMC that have similar luminosity. The findings suggest that the two binary components exchange mass. Based on the newly gained insights we reclassified the system as ON3 If^*+O5.5 V((f)).

The multi-epoch RV analysis of 19 separated observations, including 11 spaced over ten days, allowed us to constrain the orbital parameters. The best fit is achieved with a period of P = 3.07 d and a circular orbit. The orbital analysis confirms a mass ratio of q_orb = 1.5, favouring the secondary to be more massive and supporting the theory of a binary interaction. Furthermore, given the orbital and stellar parameters, we calculated the Roche radius and found that both components are filling their Roche lobes and must be in contact.

From initial abundances consistent with NGC 346, we have modelled possible stellar evolution pathways that explain the empirically derived stellar and orbital parameters of SSN 7.

Our favourite model predicts that the system has an age of ∼4 Myr. In this context, we review the age of the O-star population within the core of NGC 346 and see that no single star (except for the peculiar massive object SSN 9) is to the left of the 4Myr isochrone. We postulate that the age of the O-star population within the core of NGC 346 is older than ≳4 Myr.

The binary evolutionary models are capable of explaining the currently observed contact phase and predict that both stars are still undergoing slow Case A mass transfer.

This makes SSN 7 the most massive Algol-like system known to date, enabling the study of the physics of binary interactions in a way previously unavailable. This will shape our understanding of post-interaction massive binaries, binary BHs, and ultimately, the progenitors of GW events.

Acknowledgments

DP acknowledges financial support by the Deutsches Zentrum für Luft und Raumfahrt (DLR) grant FKZ 50 OR 2005. MJR and DP acknowledge Prof. Wolf-Rainer Hamann, Prof. Lida Oskinova, Dr. Rainer Hainich, and Dr. Heldge Todt of the Institut für Physik und Astronomie, Universität Potsdam. Dr. Andreas Sander and Dr. Varsha Ramachandran of Zentrum für Astronomie der Universität Heidelberg, Astronomisches Rechen-Institut, and Dr. Tomer Shenar of the Institut for Astronomy, University of Amsterdam, and all other contributors to the PoWR code. This publication has benefited from a discussion at a team meeting sponsored by the International Space Science Institute at Bern, Switzerland.

References

Appendix A: Additional tables

Appendix B: Light curve modelling with PHOEBE

Although it is in the survey region of the Optical Gravitational Lensing Experiment (OGLE, Udalski et al. 2015) catalogue, to date no light curve of this object has been reported. In addition, the resolution of the Transiting Exoplanet Survey Satellite (TESS) is 21″ pixel⁻¹, not fine enough to disentangle the light curve of our target from the other O stars in NGC 346. To demonstrate if detectable photometric variations based on our results are expected at all, we employed the PHOEBE code (used in Sect. 3.3.2) to create a synthetic light curve.

Typically light curves are observed in the optical. Hence, we decided to create the synthetic light curve (shown in Figure B.1) in the Johnson V band. In the PHOEBE model we used the stellar parameters of the binary components obtained from our spectroscopic analysis (see Sect. 4.1) as input parameters. Furthermore, we assume that gravitational darkening can be approximated by a power law with a coefficient of β_grav = 1, as typically used for hot stars with radiative envelopes. The atmospheres of both components are approximated by a black body, when compared to the SED of our atmospheric model it is a valid approximation for the V-band flux. Limb darkening is modelled by using a quadratic approximation and using coefficients a₁ = 0.2694 and b₁ = −0.003 for the primary and a₂ = 0.1751 and b₂ = 0.0871 for the secondary, obtained from the emergent flux distribution of our stellar atmosphere models. The effect of reflections and resulting heating in the binary is modelled using the build in the Horvat scheme that includes Lambert scattering as well as the reflection scheme of Wilson (1990). It is assumed that 90% of the incident flux will be reflected.

Fig. B.1. Synthetic predicted light curves of SSN 7 for an inclination angle of i = 16° and the orbital configuration of the binary components in the plane of sky at different phases. Upper panel: Phased synthetic light curve in the Johnson V band. Lower panel: Orbital configuration of the binary components at phases ϕ = −0.5, −0.25, 0.0, 0.25, and 0.5. The models are colour-coded by effective temperature, taking into account the physical effect mentioned in the text.

For a better understanding on how the different photometric variations are linked to the orbital configuration of the binary, we depict the orbital configurations of the stars in the lower part of Fig. B.1. It is evident that the expected maximum photometric variation of the binary components is 18 mmag. The shape of the light curve is only caused by ellipsoidal variability and not by eclipses.

Given the expected low variability, we conclude that only the upcoming powerful telescopes, such as the ELT, will be able to detect this kind of variability. Hence, we can say that our derived low inclination angle is consistent with the lack of published light curve data.

All Tables

All Figures

Fig. 1. Optical (left) and infrared (right) image of NGC346 in the SMC. Each image shows an area on the sky of approximately 3.5′×3.5′ (≈62 pc × 62 pc at d = 61 kpc). The positions of the main ionising sources of NGC 346 (see Sect. 6.4), namely SSN 7 (ON3 If*+O5.5 V((f))) and SSN 9 (O2 III(f*)), are highlighted by red and green squares, respectively. The false-colour optical image was taken with the HST ACS. Credit: A. Nota (ESA/STScI). The false-colour infrared image was taken with the James Webb Space Telescopes (JWST) Near Infrared Camera (NIRCam). Credits: NASA, ESA, CSA, O. Jones (UK ATC), G. De Marchi (ESTEC), and M. Meixner (USRA). Image processing: A. Pagan (STScI), N. Habel (USRA), L. Lenkic (USRA), and L. Chu (NASA/Ames).

In the text

	Fig. 2. Multi-epoch spectral observations of SSN 7. The X-shooter spectrum (ID 09) is shown as a solid line, the GIRAFFE spectrum (ID 17) as a dashed line, the GIRAFFE spectrum (ID 18) as a dot-dashed line, and the UVES spectrum (ID 19) as a dotted line. The shifting emission and absorption features are associated with the primary and secondary, respectively.
In the text

	Fig. 3. Observed (triangles) and synthetic (solid lines) RV curves of the primary (green) and secondary (red) component, as obtained by the PHOEBE code. The dashed black line indicates the barycentric velocity offset v0.
In the text

Fig. 4. Close-ups on the carbon and helium lines used to constrain the luminosity ratio (solid blue) compared to the combined synthetic model (dotted red). Left: observed FUSE spectra (ID 01). The weighted synthetic spectra of the primary and secondary components are shown as a dotted black and green line. Right: observed X-shooter spectra (ID 09 and 31). The unweighted synthetic spectra of the primary and secondary components are shown as a dotted black and green line.

In the text

	Fig. 5. Observed (blue) and combined synthetic (red) SED. The individual synthetic SEDs of the primary and secondary are shown as the dotted black and green line, respectively. The photometry, listed in Table 2, is shown as open squares.
In the text

	Fig. 6. Close-ups on the Balmer lines in the X-shooter spectrum (ID 09; solid blue) compared to the combined synthetic model (dotted red). The unweighted synthetic spectra of the primary and secondary components are shown as dotted black and green line.
In the text

	Fig. 7. Close-ups on the metal lines in the X-shooter spectrum (ID 09 and 31; solid blue) compared to the combined synthetic spectrum (dashed red). The unweighted synthetic spectra of the primary and the secondary component are shown as a dotted black and green line, respectively.
In the text

	Fig. 8. Close-up on the iron forest of the HST FOS (ID 03) spectrum (solid blue) compared to the combined synthetic model flux (dashed red). The weighted synthetic model spectra of the primary and secondary are shown as a dotted black and green line, respectively. The iron forest is used to check the iron abundance and to constrain the microturbulence within the stellar atmosphere.
In the text

Fig. 9. Close-ups on the key wind diagnostic lines (solid blue) in the FUV, UV and optical spectra overplotted by the combined synthetic model (dashed red). In the upper panels the weighted and in the lower panels the unweighted synthetic spectra of the primary and secondary components are shown as dotted black and green line, respectively. Upper left: the O VI line of the FUSE (ID 1) spectrum. This line is sensitive to X-rays. Upper middle: the N V resonance doublet in the HST FOS (ID 3) spectrum. Upper right: the C IV resonance doublet in the HST FOS (ID 3) spectrum. Lower left: the He IIλ 4686 line in the X-shooter spectrum (ID 9), showing strong wind emission features from the primary. Lower right: Hα of the X-shooter spectrum (ID 31), again, the primary dominates the wind emission of this line.

In the text

Fig. 10. Evolutionary tracks of the donor (left) and the accretor (right) colour-coded by their surface H-abundance. The track of the mass gainer is only colour-coded until the primary dies and then continues as a solid white line, assuming that the binary is not disrupted by a supernova explosion. The tracks are overlayed by black dots that show equidistant time steps of 30 000 yr to highlight phases in which the stars spend most of their time. The spectroscopic results are shown as green and red triangles surrounded by error ellipses. Start and end phases of fast and slow Case A, as well as Case AB mass-transfer are indicated by arrows. The iso-contours of equal radii are indicated in the background as dashed grey lines.

In the text

Fig. 11. Observed mass-loss rates of O stars in NGC 346 (including SSN 7) as a function of stellar luminosity compared to theoretical and empirical relations. The position of both binary components of SSN 7 are marked by black stars. The data is complemented by O stars from Rickard et al. (2022). Apparently single stars are shown as filled blue squares and known muliple systems as open blue squares. In the background the empirical relation of main-sequence O-type stars from Rickard et al. (2022) is included as dashed black line. In addition, we show the theoretical relations of Vink et al. (2001) and Björklund et al. (2023) for main-sequence stars as dashed yellow and dashed red line, respectively.

In the text

	Fig. 12. log QH to spectral type. SMC sample stars analysed by Ramachandran et al. (2019) are included (red diamonds) as well as LMC stars from Ramachandran et al. (2018; blue diamonds) and O stars from the core of NGC 346 from Rickard et al. (2022; black diamonds). The two components of the SSN 7 binary are shown separately.
In the text

Fig. 13. HR diagram showing the position of the two components of SSN 7 as determined by our spectrographic analysis. Black stars are the two SSN 7 components. Blue symbols are the other O stars within the central core of NGC 346 (Rickard et al. 2022). Open symbols are noted spectral binaries (Dufton et al. 2019). Dotted lines are isochrones with Z = 0.002 (Georgy et al. 2013) and an initial rotation of 100 km s−1.

In the text

Fig. B.1. Synthetic predicted light curves of SSN 7 for an inclination angle of i = 16° and the orbital configuration of the binary components in the plane of sky at different phases. Upper panel: Phased synthetic light curve in the Johnson V band. Lower panel: Orbital configuration of the binary components at phases ϕ = −0.5, −0.25, 0.0, 0.25, and 0.5. The models are colour-coded by effective temperature, taking into account the physical effect mentioned in the text.

In the text