11institutetext: Institut für Physik und Astronomie, Universität Potsdam,
Karl-Liebknecht-Str. 24/25, 14476 Potsdam, Germany 22institutetext: Armagh
Observatory and Planetarium, College Hill, Armagh BT61 9DG, Northern Ireland,
UK 33institutetext: Astronomisches Rechen-Institut, Zentrum für Astronomie der
Universität Heidelberg, Mönchhofstr. 12-14, 69120 Heidelberg, Germany
44institutetext: Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, 3001
Leuven, Belgium 55institutetext: Department of Physics and Astronomy,
University College London, Gower Street, London WC1E 6BT, UK 66institutetext:
Centro de Astrobiología. CSIC-INTA. Campus ESAC. Camino bajo del castillo s/n.
E-28 692 Villanueva de la Cañada, Madrid, Spain.

# The earliest O-type eclipsing binary in the Small Magellanic Cloud, AzV 476:
A comprehensive analysis reveals surprisingly low stellar masses

D. Pauli1 L. M. Oskinova1 W.-R. Hamann1 V. Ramachandran1,3 H. Todt1 A. A. C.
Sander2,3 T. Shenar4 M. Rickard1,5 J. Maíz Apellániz6 R. Prinja5

(Received ; Accepted)

###### Abstract

Context. Massive stars at low metallicity are among the main feedback agents
in the early Universe and in present-day star forming galaxies. When in
binaries, these stars are potential progenitors of gravitational-wave events.
Knowledge of stellar masses is a prerequisite to understanding evolution and
feedback of low-metallicity massive stars.

Aims. Using abundant spectroscopic and photometric measurements of an
outstandingly bright eclipsing binary, we compare its dynamic, spectroscopic,
and evolutionary mass estimates and develop a binary evolution scenario.

Methods. We comprehensively studied the eclipsing binary system, AzV 476, in
the Small Magellanic Cloud (SMC). The light curve and radial velocities were
analyzed to obtain the orbital parameters. The photometric and spectroscopic
data in the UV and optical were analyzed using the Potsdam Wolf-Rayet (PoWR)
model atmospheres. The obtained results are interpreted using detailed binary-
evolution tracks including mass transfer.

Results. AzV 476 consists of an O4 IV-III((f))p primary and an O9.5: Vn
secondary. Both components have similar current masses (${20\,M_{\odot}}$ and
${18\,M_{\odot}}$) obtained consistently from both the orbital and
spectroscopic analysis. The effective temperatures are
$42\text{\,}\mathrm{k}\mathrm{K}$ and $32\text{\,}\mathrm{k}\mathrm{K}$,
respectively. The wind mass-loss rate of
${\log(\dot{M}/({\,M_{\odot}\,\mbox{yr}^{-1}}))=-6.2}$ of the primary is a
factor of ten higher than a recent empirical prescription for single O stars
in the SMC. Only close-binary evolution with mass transfer can reproduce the
current stellar and orbital parameters, including orbital separation,
eccentricity, and the rapid rotation of the secondary. The binary evolutionary
model reveals that the primary has lost about half of its initial mass and is
already core helium burning.

Conclusions. Our comprehensive analysis of AzV 476 yields a consistent set of
parameters and suggests previous case B mass transfer. The derived stellar
masses agree within their uncertainties. The moderate masses of AzV 476
underline the scarcity of bright massive stars in the SMC. The core helium
burning nature of the primary indicates that stripped stars might be hidden
among OB-type populations.

###### Key Words.:

binaries: eclipsing - binaries: close - binaries: spectroscopy - stars: early-
type - stars: fundamental parameters - stars: individual: AzV 476

##  1 Introduction

The most important parameter defining the evolution of a star is its mass. An
often reported problem in stellar astrophysics is the “mass discrepancy”
problem, which refers to the inconsistent masses derived from spectroscopy,
evolutionary tracks, and, for binaries, from orbital motions (Herrero et al.
1992; Weidner & Vink 2010; Markova & Puls 2015). The analysis of O-type stars
located in the Galaxy and the Large Magellanic Cloud (LMC) performed by
Weidner & Vink (2010) reveal good agreement between the aforementioned mass
estimates, but the study of Markova et al. (2018) suggests a mass discrepancy.
Mahy et al. (2020) investigate a sample of O-type binaries in the LMC and find
good agreement between spectroscopic and dynamic mass estimates while the
evolutionary masses are at odds. These latter authors suggest previous binary
interactions as a possible solution for this discrepancy.

However, the mass discrepancy problem has not yet been studied at
metallicities lower than ${\lesssim 1/2\,Z_{\odot}}$. Sufficiently low
metallicity is offered by the nearby Small Magellanic Cloud (SMC) galaxy
($Z_{\rm SMC}\approx 1/7\,Z_{\odot}$;Hunter et al. (2007); Trundle et al.
(2007)). Precise stellar masses in the SMC allow us to tackle questions about
stellar evolution and feedback at low metallicity.

Stars with spectral types around O2-4 are expected to be very massive with
$M_{\ast}\gtrsim 50\,M_{\odot}$ (Martins & Palacios 2021). However, the true
masses of the early-type O stars are only poorly known, and therefore spectral
types and masses might be falsely mapped. So far, only a couple of the SMC
eclipsing binaries with early spectral types have been studied. Morrell et al.
(2003) investigated the O6V+O4/5III(f) system Hodge 53-47 (alias MOA
J010321.3-720538) and found dynamic masses of ${\approx 26\,M_{\odot}}$ and
${\approx 16\,M_{\odot}}$ for the primary and secondary, respectively. The
eclipsing binary, OGLE SMC-SC10 108086, studied by Abdul-Masih et al. (2021)
contains even less massive stars with ${\approx 17\,M_{\odot}}$ and ${\approx
14\,M_{\odot}}$. Only one eclipsing multiple stellar system in the SMC, HD
5980, appears to contain stars with masses $\gtrsim 30\,M_{\odot}$.
Koenigsberger et al. (2014) and Hillier et al. (2019) studied this system in
detail and found an inner eclipsing binary consisting of an LBV and a Wolf-
Rayet (WR) star, and a third O-type supergiant star with a potential fourth
companion. The orbital masses of the LBV and the WR star are ${\approx
61\,M_{\odot}}$ and ${\approx 66\,M_{\odot}}$, respectively. Koenigsberger et
al. (2014) suggest that there was little or no mass transfer, and that the WR
star has formed via quasi-chemically homogeneous evolution. These examples
highlight the complexity of massive star evolution and the need for advanced
studies on the mass discrepancy problem in the high-mass regime.

Table 1: List of all spectra of AzV 476 used in this work and their associated
orbital phase. RVs and the wavelengths ranges and lines used for their
measurements are also listed. The listed RVs are already corrected for the
barycentric motion and the velocity of the NGC 456 complex with
$\varv_{\mathrm{SMC}}=$152\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$$.

spectral ID | Instrument | Wavelength | MJD(a) | Phase $\phi$ (b) | RV1 | RV2
---|---|---|---|---|---|---
|  | [Å] | [d] |  | [$\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$] | [$\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$]
1 | FUSE |  $9501150\text{\,}\AA$ | 52478.3849 | $\,\,\,\,0.0192$ | — | —
2 | HST/COS | $11781777\text{\,}\AA$ | 56185.2387 | $-0.2305$ | $\,\,\,\,\,194\pm\,\,5\,^{(c)}$ | —
3 | HST/STIS | $11401735\text{\,}\AA$ | 59022.5673 | $-0.3124$ | $\,\,\,\,\,\,\,\,91\pm\,\,5\,^{(c)}$ | $\,\,\,\,\,\,\,\,\,-86\pm 12\,^{(c)}$
4 | HST/STIS | $15742673\text{\,}\AA$ | 57321.6698 | $\,\,\,\,0.0968$ | $\,\,\,\,-55\pm\,\,5\,^{(d)}$ | $\,\,\,\,\,\,\,\,\,\,\,\,\,64\pm 45\,^{(d)}$
5 | X-SHOOTER | $300010\,000\text{\,}\AA$ | 59163.2501 | $-0.2928$ | $112\pm 14$ | $-110\pm 10$
6 | UVES | $300011\,000\text{\,}\AA$ | 57703.2868 | $-0.1610$ | $223\pm\,\,9\,$ | $-171\pm 15$
7 | UVES | $300011\,000\text{\,}\AA$ | 57748.0479 | $-0.3823$ | $\,\,\,60\pm 16$ | $\,\,\,-42\pm 17$
8 | UVES | $300011\,000\text{\,}\AA$ | 57749.0558 | $-0.2747$ | $126\pm\,\,8\,$ | $-104\pm 29$
9 | UVES | $300011\,000\text{\,}\AA$ | 57749.0923 | $-0.2708$ | $141\pm 16$ | $\,\,\,-95\pm 37$
10 | UVES | $300011\,000\text{\,}\AA$ | 57750.1132 | $-0.1618$ | $201\pm 13$ | $-184\pm 34$
11 | UVES | $300011\,000\text{\,}\AA$ | 57750.1492 | $-0.1580$ | $207\pm\,\,9\,$ | $-186\pm 43$
12 | UVES | $300011\,000\text{\,}\AA$ | 57751.0782 | $-0.0588$ | $148\pm 12$ | $\,\,\,-71\pm 37$

(a) Mid-exposure in HJD - 2400000.5 (b) Calculated with Eq. (2). (c) Obtained
from a fit over the range $13601405\text{\,}\AA$. (d) Obtained from a fit over
the range $21002230\text{\,}\AA$.

For single stars, mass estimates rely on spectroscopic diagnostics or
comparison with evolutionary tracks. Currently, standard stellar evolutionary
models predict that the most massive stars, either single or in binary
systems, start their lives on the main sequence (MS) as early-type O stars.
Single stars expand and evolve away from the MS. During this phase, stars
undergo strong mass loss where they might lose their entire hydrogen-rich
envelope, revealing the helium core and becoming WR stars. In the case of
close binaries, evolutionary models predict that the expanding star is likely
to interact with its companion. In the majority of cases, the expanding star
will transfer its envelope to the companion, and become a binary-stripped
helium star, possibly also with a WR-type spectrum (Dionne & Robert 2006;
Shenar et al. 2020; Götberg et al. 2020).

The SMC hosts a handful of WR stars that have relatively high masses ranging
from $10\,M_{\odot}$ to $60\,M_{\odot}$ (Shenar et al. 2016). These stars are
so massive that their hydrogen-burning progenitors must have been early O-type
(or WNL/Of) stars. However, in their recent study of the SMC OB-type
population, Ramachandran et al. (2019) reveal a strong deficiency of massive
stars ($>30\,M_{\odot}$) close to the MS in the upper part of the empiric
Hertzsprung–Russell diagram (HRD) (see Holgado et al. 2020; Ramachandran et
al. 2018, for studies of Galactic and the LMC O star populations). O-type
stars are hot and luminous and are therefore easily detectable. Hence, a
paucity of the earliest O-type stars in the SMC cannot be explained either by
selection effects or stellar evolution scenarios (Schootemeijer et al. 2021).
Thus, the deficiency of most massive O stars strongly questions the formation
process of the WR stars in the SMC, as well as our basic understanding of
stellar evolution in low-metallicity environments. In this context it is
crucial to quantify the budget of massive stars ($\gtrsim 30\,M_{\odot}$) in
the SMC.

Poorly constrained processes that affect the lives of massive stars include
radiatively driven winds, stellar envelope inflation, core overshooting, and
rotationally induced internal mixing. In particular their dependence on mass
and metallicity are not yet fully understood. These effects can drastically
alter the evolution of a star. Studies of objects in low-metallicity
environments are needed to constrain these metallicity-dependent effects and
thus allow the establishment of reliable stellar evolutionary tracks.

To address these outstanding questions, we selected one of the earliest
subtype O stars in the SMC that is an eclipsing binary, allowing us to
estimate its mass by various methods. The subject of this paper, AzV 476, is
located in the cluster NGC 456 in the SMC Wing. The cluster contains an active
star-forming region and hosts young stellar objects (Muraoka et al. 2017). AzV
476 is embedded in a H ii region. The system was identified as an eclipsing
binary by the Optical Gravitational Lensing Experiment (OGLE), which monitors
stellar variability in the SMC and LMC (Pawlak et al. 2016). The primary star
was classified by Massey et al. (2005) as an O2-3 V star based on its optical
spectral appearance and is therefore one of the earliest type stars in the
entire SMC. Here we present the first consistent analysis of AzV 476, based on
photometric and spectroscopic data, and in particular accounting for its
binary nature. Newly obtained high-resolution UV and optical data give us the
opportunity to perform such an analysis, yielding estimates on masses and
stellar and wind parameters.

This paper is structured as follows. In Sect. 2 we describe the observations
and the known stellar and orbital parameters used in the analysis below. The
results of the orbital and spectral analysis as well as those from
evolutionary modeling are presented in Sects. 3, 4, and 5, respectively. Their
implications, similarities, and disagreements on the stellar masses and other
stellar parameters are discussed in Sect. 6, and conclusions are given in
Sect. 7.

##  2 Observations

###  2.1 Spectroscopy

Over the last decade, a handful of AzV 476 spectra covering the UV, optical
(VIS), and infrared (IR) have been obtained, yielding multi-epoch data well
suited for measuring radial velocities (RVs) and spectroscopic analysis. We
used a spectroscopic dataset consisting of 12 spectra and Table 1 gives a
brief overview of these, their covered wavelength ranges, their observed date,
and their associated orbital phase. In the remainder of this paper, we refer
to the individual observed spectra by their ID in Table 1.

The wavelength regime $9501150\text{\,}\AA$ is covered by an archival FUSE
(Oegerle et al. 2000) observations111the shorter wavelengths in the FUSE range
are largely contaminated by interstellar features that yield no information
about the stellar parameters. The FUSE spectrum (ID 1) was taken with a total
exposure time of $21\,567\text{\,}\mathrm{s}$ and a resolving power of
${R\approx 20\,000}$. Incidentally, the FUSE observation was taken close to
the primary eclipse. This FUSE spectrum has a known but unsolved calibration
issue. Therefore, we only used the LiF1A, LiF2A, and SiC2A channels, which
appear to be the least affected by this problem. The spectrum is rectified by
division through the combined continuum flux of our models.

AzV 476 is part of the ULLYSES program222https://ullyses.stsci.edu/. It was
observed with the HST/COS (Hirschauer et al. 2021) using the G130M
($11781472\text{\,}\AA$) and G160M ($13831777\text{\,}\AA$) medium-resolution
gratings (ID 2), ${R\approx 19\,000}$. The two spectra were taken sequentially
with exposure times of $330\text{\,}\mathrm{s}$ and $1100\text{\,}\mathrm{s}$,
respectively. In addition, the star was re-observed in the UV with the
HST/STIS spectrograph (Branton et al. 2021) using the E140M echelle graiting
(ID 3) as part of the HST program 15837 (PI Oskinova). The spectrograph covers
a wavelength regime of $11401735\text{\,}\AA$. The exposure time was
$2707\text{\,}\mathrm{s}$ and the final resolving power is ${R\approx
45\,800}$.

There is another spectrum in the ULLYSES program (ID 4) taken with the
HST/STIS spectrograph using the E230M echelle gratings covering a wavelength
range $15742673\text{\,}\AA$. The exposure time was $2820\text{\,}\mathrm{s}$
and the resolving power is ${R=30\,000}$.

For the optical and near-IR range, we use the publicly available spectra taken
with the X-SHOOTER spectrograph (Vernet et al. 2011) mounted on the ESO Very
Large Telescope (VLT). The spectrum (ID 5) was taken as part of the ESO
106.211Z program, which is part of the XSHOOTU program (their Paper I; Vink
and the XShootU Collaboration, in prep.). The X-SHOOTER spectrograph consists
of three different spectroscopic arms, which are optimized for the wavelength
ranges in the UBV ($30005550\text{\,}\AA$), VIS ($530010\,000\text{\,}\AA$),
and near-IR ($10\,00025\,000\text{\,}\AA$). The resolving powers are
${R\approx 6\,600}$, ${R\approx 11\,000,}$ and ${R\approx 8\,000}$,
respectively. The spectra were obtained with exposure times of
$750\text{\,}\mathrm{s}$, $820\text{\,}\mathrm{s}\mathrm{,}$ and
$300\text{\,}\mathrm{s}$ for the UBV, VIS, and NIR arm, respectively.

The remaining seven spectra (ID $612$) used for our analysis are taken with
the UVES spectrograph (Dekker et al. 2000) mounted on the ESO VLT. Each
spectrum was taken with the DIC1 setup covering the wavelength ranges of
$30004000\text{\,}\AA$ and $500011\,000\text{\,}\AA$. The exposure times are
about $2895\text{\,}\mathrm{s}$ for both spectrographs and have a resolving
power of ${R\approx 65\,000}$ and ${R\approx 75\,000}$, respectively. The UVES
spectra were rectified by hand.

###  2.2 Photometry

Table 2: UBVRIJHK photometry of AzV 476. Band | Apparent magnitude
---|---
| [mag]
U | $12.49\pm 0.04$
B | $13.54\pm 0.07$
V | $13.48\pm 0.01$
VOGLE | $13.49$
R | $13.72\pm 0.06$
I | $13.54\pm 0.20$
IOGLE | $13.56$
J | $13.70\pm 0.03$
H | $13.72\pm 0.04$
K | $13.86\pm 0.06$
G | $13.51\pm 0.06$
GBP | $13.46\pm 0.01$
GRP | $13.57\pm 0.01$

The UBI photometry is adopted from the catalog of the SMC stellar population
(Bonanos et al. 2010). For the VR photometry, we use the magnitudes from the
fourth United States Naval Observatory (USNO) CCD Astrograph Catalog (UCAC4)
(Zacharias et al. 2013). For completeness, we compare the values of the V- and
I-band magnitudes to those listed in the IVth OGLE Collection of Variable
Stars (Pawlak et al. 2016) (VOGLE and IOGLE). Unfortunately, the OGLE
magnitudes are published without error margins. Nonetheless, we find that they
are in agreement with the V- and I-band magnitude from Zacharias et al. (2013)
and Bonanos et al. (2010). JHK photometry is from the 2MASS catalog (Cutri et
al. 2003). Additionally, we use the recent EDR3 Gaia photometry (Gaia
Collaboration et al. 2016, 2021). A total list of the used magnitudes is shown
in Table 2. Massey et al. (2005) estimated the extinction towards AzV 476 to
be $E_{\mathrm{B-V}}=$0.28\text{\,}\mathrm{m}\mathrm{a}\mathrm{g}$$ based on
averaging the color excesses in $B-V$ and $U-B$ based on the spectral type.
From our spectral energy distribution (SED) fit we find better agreement when
using a lower extinction of
$E_{\mathrm{B-V}}=$0.26\text{\,}\mathrm{m}\mathrm{a}\mathrm{g}$$ (see Sect.
4.1). The slight difference can be due to the use of different reddening laws.

For the light curve modeling, we use the OGLE I-band photometry from the IVth
OGLE Collection of Variable Stars (Pawlak et al. 2016). The photometric data
are taken over the period from May 2010 to January 2014. We do not use the
OGLE V-band photometry because it contains too few data points and therefore
cannot be used to resolve the eclipses.

AzV 476 was observed by the TESS space telescope in 2018 (sector 2) and 2020
(sectors 27 and 28). However, the relatively low spatial resolution of TESS
($21\text{\,}\mathrm{{}^{\prime}}\mathrm{{}^{\prime}}\,\mathrm{p}\mathrm{x}^{-1}$)
precludes accurate point source photometry in the crowded region around our
target. Therefore, the TESS data are not used for light-curve modeling.
However, we use the TESS data to improve the ephemeris (see Sect. 3.1.1).

###  2.3 Distance and location

Our target, AzV 476 is part of the NGC 456 cluster which is located in the SMC
Wing. The distance to the SMC Wing was estimated by Cignoni et al. (2009)
${d\approx$55\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$}$ corresponding to a
distance modulus of
$\mathrm{DM}=$18.7\text{\,}\mathrm{m}\mathrm{a}\mathrm{g}$$. This is in
agreement with Nidever et al. (2013) who use red clump stars to study the
structure of the SMC, including the SMC Wing, and find that it follows a
bimodal distribution with a near component at a distance of
${d\approx$55\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$}$ and a far component at
${d\approx$67\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$}$. Tatton et al. (2021)
confirm this distance, but they speculate that young structures, such as the
NGC 456 cluster, do not trace substructures that are associated with the
intermediate-age populations and might be located in front of them (see Sect.
6.1). In this work we adopt a distance of
${d=$55\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$}$.

The observed spectra of AzV 476 are corrected for barycentric motion, which
was calculated with the tool described in Wright & Eastman (2014).
Additionally, we shifted the spectra by the RV of the NGC 456 complex in the
SMC $\varv_{\mathrm{SMC}}=$152\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$$
determined by fitting Gaussians to several interstellar medium (ISM) lines
that we associate with the environment of AzV 476. The RV is not uniform
throughout the different regions in the SMC and our finding is in agreement
with the results of De Propris et al. (2010).

##  3 Analysis of the binary orbit

###  3.1 Method: Eclipse light-curve and RV curve modeling

####  3.1.1 Light curve

Our target, AzV 476, is listed in the IVth OGLE Collection of Variable Stars
(Pawlak et al. 2016) with an orbital period of
${P_{\mathrm{OGLE}}=$9.366\,319\,8\text{\,}\mathrm{d}$}$ and an epoch of the
primary eclipse of ${T_{0\mathrm{,\,\,OGLE}}=2457002.7608}$ in Heliocentric
Julian Date (HJD). Unfortunately, the orbital period and the epoch of the
primary eclipse are given without error margins. The observations that
constitute the OGLE I-band light curve were taken with a cadence of
$\approx$2\text{\,}\mathrm{d}$$, and therefore the individual eclipses are not
well resolved. In contrast, the TESS light curve has a much finer time
coverage of $\approx$5\text{\,}\mathrm{h}\mathrm{,}$$ allowing better
constraint of the orbital period and the epoch of the primary eclipse.

The date at which an eclipse occurs can be expressed as

$T(n)=P\cdot n+{T}_{0},$ (1)

where $T$ is the time of the primary eclipse at orbital cycle $n$. We fitted
Gaussians to the primary eclipses in the TESS light curve; the corresponding
orbital cycles and dates of all eclipses are listed in Table 3.

Table 3: Dates of the primary eclipses in the TESS light curve.

orbital cycle $n$ | MJD(a)
---|---
| [d]
145 | $58360.4633\pm 0.0027$
146 | $58369.8256\pm 0.0030$
147 | $58379.1905\pm 0.0029$
218 | $59044.2165\pm 0.0016$
219 | $59053.5942\pm 0.0016$
220 | $59062.9590\pm 0.0014$
221 | $59072.8303\pm 0.0014$
222 | $59081.6913\pm 0.0015$

(a) ${\mathrm{MJD}=\mathrm{HJD}-2400000.5}$. The TESS data are given in
${\mathrm{TBJD}=\mathrm{BJD}-2457000.0}$ and that we converted the
$\mathrm{BJD}$ to $\mathrm{HJD}$ to be comparable to the OGLE data.

Using this procedure, we obtain ${P=$9.366\,65\pm
0.000\,25\text{\,}\mathrm{d}$}$ and ${{T}_{0}=2457002.7968\pm 0.0052}$.

Figure 1 shows a part of the phased light curve centered on the primary
eclipse for the two different ephemerides, those published in the OGLE catalog
and those we obtain in this work. As can be seen, our newly obtained ephemeris
adequately describes the OGLE as well as the TESS light curves, i.e., the data
set that covers $>$10\text{\,}\mathrm{y}\mathrm{r}$$ of observations.

Figure 1:  Phased OGLE I-band and TESS light curves around the primary
eclipse. Left panel: Light curve phased according to the ephemeris from the
OGLE catalog. Right panel: Light curve phased according to the ephemeris we
derive in this work by fitting the primary’s eclipses of the TESS light curve.

In order to convert a date $t$ at which a spectrum was taken to a phase
$\Phi,$ the following formula is used

$\phi(t)=\cases{}{t-{T}_{0}}{P}\mod{1}&\text{if }\phi<0.5\\\
\rule{0.0pt}{19.91684pt}\dfrac{t-{T}_{0}}{P}\mod{1}-1\text{otherwise.}{}$ (2)

With the improved ephemerides and their small error margins, the uncertainties
on the phases are negligible and are therefore not taken into account here.
The OGLE I-band light curve is shown in the upper panel of Fig. 2. The
conjunctions are at phases ${\phi=0.0}$ and ${\phi=-0.34,}$ implying that the
orbit is eccentric.

Figure 2: Upper panel: Phased OGLE I-band light curve of AzV 476 (blue dots)
and the best fit obtained with the PHOEBE code (red line). Above the light
curve, the spectral IDs of all used spectra (see Table 1) are indicated. Lower
panel: Observed (triangles) and fitted (solid lines) RV curves for the primary
(green) and secondary (red). The fits of the light-curve fit and the RVs are
consistently obtained by the PHOEBE code.

####  3.1.2  RVs

To measure the RVs, we use a Markov chain Monte Carlo (MCMC) method combined
with a least-square fitting method. In the MCMC method, the individual
synthetic spectra of the primary and secondary (see Sect. 4.1) are shifted by
different RVs. The combined synthetic spectrum is then compared to the
observation. Using a least-square likelihood function we estimate the quality
of the used RVs. The MCMC method quickly explores a large parameter space of
different RVs until it converges toward the true solution. In the vicinity of
the true solution, the MCMC method calculates the probabilities of different
combinations of the RVs. This yields the final probability distribution around
the true solution. A more detailed explanation is given in Appendix A.

Because the final probability distribution obtained with the MCMC method is
not necessarily a Gaussian, we quote the error as the 68% confidence interval.
The PHOEBE code cannot handle asymmetric errors, and therefore we only use the
larger margin of the probability distribution as it is the safer choice. The
uncertainties are included in the eclipsing binary modeling as described in
Sect. 3.1.3.

The primary dominates the emission and absorption lines. In order to avoid
uncertainties due to the intrinsic variability of lines formed in the stellar
wind, the RVs listed in Table 1 are the averaged values of the RVs obtained
from selected individual lines. A more detailed list of the RVs obtained from
the individual lines in the different optical spectra is given in Tables 7 and
8.

All the absorption lines that are associated with the secondary show a
contribution from the primary. Furthermore, the depths of these absorption
lines are at the level of the noise, which introduces additional
uncertainties; these are reflected in the larger error margins. Two selected
optical spectra obtained at different phases are depicted in Fig. 3 to
demonstrate how the spectral lines associated with the different binary
components shift.

Figure 3:  X-SHOOTER spectrum (ID 5 in Table 1) displayed in blue and one of
the UVES spectra (ID 10 in Table 1) in red. The spectra are convolved with a
Gaussian with an ${\mathrm{FWHM}=$0.4\text{\,}\AA$}$ to reduce the noise and
to make the RV shifts visible. The region containing $\mathrm{H\gamma}$, He i
$\lambda 4387,$ and He i $\lambda 4471$ lines is shown. In the spectrum shown
by the red line, the primary’s spectrum is redshifted (see
$\mathrm{H\gamma}$), while the secondary’s spectrum (broadened He i lines) is
blueshifted. We note that the primary also partially contributes to the He i
lines.

####  3.1.3 Modeling with PHOEBE

The Physics of Eclipsing Binaries (PHOEBE) v.2.3 modeling software (Prša &
Zwitter 2005; Prša et al. 2016; Horvat et al. 2018; Jones et al. 2020; Conroy
et al. 2020) is employed to derive dynamic masses as well as to obtain
measures of the stellar parameters independently from the spectroscopic model.
The simultaneous fitting of the RV and light curve is done with the `emcee`
sampler (Foreman-Mackey et al. 2013).

To reduce the parameter space, we fix the orbital period and the epoch of the
primary eclipse to the values we derive from the TESS light curve; see Sect.
3.1.1. As we are using only one passband, the temperatures of the components
cannot be obtained reliably. Therefore, the temperatures of the primary and
secondary are fixed to the results from the spectral analysis (see Sect. 4.1),
${T_{\mathrm{eff,\,1}}=$42\text{\,}\mathrm{k}\mathrm{K}$}$ and
${T_{\mathrm{eff,\,2}}=$32\text{\,}\mathrm{k}\mathrm{K}$}$, respectively.

The PHOEBE code assumes synchronous stellar rotation. The actual fast rotation
of the secondary with $\varv\sin
i=$425\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$$ is therefore not
consistently taken into account for modeling the secondary eclipse.
Gravitational darkening is modeled by a power law with a coefficient of
$\beta_{\mathrm{grav}}=1$, as recommended for radiative envelopes in the
PHOEBE documentation333www.phoebe-project.org. Given the period of
$\sim$10\text{\,}\mathrm{d}$$, this renders gravitational darkening
unimportant.

The atmospheres of both components are approximated by a blackbody. We
compared the blackbody flux to the SED of our atmospheric model and find that
it is a valid approximation for the I-band flux. We calculated the emergent
flux distribution using our spectral models (see Sect. 4.1) and fitted
different types of limb-darkening laws (e.g., Diaz-Cordoves & Gimenez 1992).
We find that the limb-darkening law that best describes the primary and
secondary is a quadratic approximation in the form of

$I(\mu)=I(1)\,[1-a_{i}\,(1-\mu)-b_{i}\,(1-\mu)^{2}],$ (3)

where $\mu=\cos\theta$ is the cosine of the directional angle $\theta$, and
$a_{i}$ and $b_{i}$ are the limb-darkening coefficients of each stellar
component $i$. For the primary, the best fit is achieved with coefficients
$a_{1}=0.2032$ and $b_{1}=0.0275$, while for the secondary $a_{2}=0.1668$ and
$b_{2}=0.0802$ are required.

In a binary where both stellar components have large radii, on the order of
$\gtrsim 20\%$ of their separation, and rather similar temperatures, the so-
called “reflection effect” becomes important (Wilson 1990). This effect
accounts for the irradiation of the surface by the other component (Prša 2011;
Prša et al. 2016). As AzV 476 is a close binary with two hot O-stars, the
reflection effect is modeled with two reflections. The effect of ellipsoidal
variability due to tidal interaction, which induces periodic variations in the
light curve, is important in close binary systems with orbital periods on the
order of a few days (Mazeh 2008). However, as none of the binary components
are close to filling their Roche lobe (see Table 5) and the mass ratio is not
extreme ($q_{\mathrm{orb}}=0.89$) this effect is expected to be negligible
(Gomel et al. 2021).

Figure 4:  Same as Fig. 3, but now for the region around the He i $\lambda
3819$ and $\mathrm{H\epsilon}$ lines. Most of the lines, including He i
$\lambda 3926$, are redshifted and associated with the primary, while the He i
$\lambda 3819$ line which is associated with the secondary is blueshifted.

###  3.2 Resulting binary parameters

Table 4: Orbital parameters obtained from the RV and light curves by the
PHOEBE code

Parameter | Value
---|---
$P\,(\mathrm{d})$ | $9.36665\,\,\,\,(\mathrm{fixed})$
HJD0 | $7002.7968\,\,\,\,(\mathrm{fixed})$
$e$ | $0.240^{+0.002}_{-0.002}$
$\omega_{0}\,(\mathrm{{}^{\circ}})$ | $19^{+2}_{-1}$
$q_{\mathrm{orb}}$ | $0.89^{+0.06}_{-0.06}$
$\gamma$ | $12^{+3}_{-3}$
$i\,(\mathrm{{}^{\circ}})$ | $77.9^{+0.3}_{-0.3}$
$a\,({R_{\odot}})$ |  $63^{+2}_{-2}$

Table 5: Stellar parameters obtained from the RV and light curves by the
PHOEBE code

Parameter | Primary | Secondary
---|---|---
$T_{\mathrm{eff}}\,\,[\mathrm{kK}]$ | $42\,\,(\mathrm{input})$ | $32\,\,(\mathrm{input})$
$K\,\,[\mathrm{km\,s^{-1}}]$ | $161^{+13}_{-13}$ | $181^{+14}_{-14}$
$M_{\mathrm{orb}}\,\,[{M_{\odot}}]$ | $20.2^{+2.0}_{-2.0}$ | $18.0^{+1.8}_{-1.8}$
$R\,\,[{R_{\odot}}]$ | $10.7^{+0.4}_{-0.4}$ | $7.0^{+0.4}_{-0.4}$
$R_{\mathrm{RL}}\,\,[{R_{\odot}}]$ | $24.5^{+0.9}_{-0.9}$ | $23.3^{+0.8}_{-0.8}$
$\log g\,\,[$\mathrm{c}\mathrm{m}\,\mathrm{s}^{-2}$]$ | $3.69^{+0.06}_{-0.06}$ | $4.00^{+0.07}_{-0.07}$
$\log L\,\,[L_{\odot}]$ | $5.51^{+0.13}_{-0.13}$ | $4.67^{+0.22}_{-0.22}$
$M_{\mathrm{bol}}\,[\mathrm{mag}]$ | $-8.94^{+0.32}_{-0.32}$ | $-6.83^{+0.56}_{-0.56}$
$f_{2}/f_{1}\,\,\,\,$(OGLE I-band) |  $0.31^{+0.05}_{-0.05}$

Table 6: Summary of the stellar parameters of both stellar components obtained
from the different methods.

| Spectroscopic analysis | Orbital analysis(a) | Single star evolution(b) | Binary evolution(c)
---|---|---|---|---
| star 1 | star 2 | star 1 | star 2 | star 1 | star 2 | star 1 | star 2
$T_{\mathrm{eff}}\,\,[$\mathrm{k}\mathrm{K}$]$ | $42^{+3}_{-3}$ | $32^{+4}_{-4}$ |  $42$ (fix) |  $32$ (fix) | $40.1^{+1.7}_{-1.5}$ | $32.5^{+1.8}_{-1.8}$ | ${42.0}$ | $31.6$
$\log\,g\,\,[$\mathrm{c}\mathrm{m}\,\mathrm{s}^{-2}$]$ | $\,\,3.7^{+0.2}_{-0.1}$ | $\,\,4.0^{+0.2}_{-0.2}$ | $3.69^{+0.06}_{-0.06}$ | $4.00^{+0.07}_{-0.07}$ | $3.80^{+0.09}_{-0.08}$ | $4.03^{+0.10}_{-0.11}$ | $3.56$ | $4.01$
$\log\,L\,\,[L_{\odot}]$ | $5.65^{+0.2}_{-0.2}$ | $4.75^{+0.2}_{-0.2}$ | $5.51^{+0.13}_{-0.13}$ | $4.67^{+0.22}_{-0.22}$ | $5.60^{+0.18}_{-0.13}$ | $4.67^{+0.11}_{-0.14}$ | $5.55$ | $4.66$
$R\,\,[R_{\odot}]$ | $12.6^{+1.0}_{-1.0}\,$ | $\,\,7.8^{+2.0}_{-2.0}$ | $10.7^{+0.4}_{-0.4}$ | $7.0^{+0.4}_{-0.4}$ | $12.94^{+2.2}_{-1.7}$ | $6.57^{+0.8}_{-0.9}$ | $11.3$ | $7.2$
$M\,\,[M_{\odot}]$ | $\,\,29^{+17}_{-11}$ | $\,\,22^{+13}_{-8}$ | $20.2^{+2.0}_{-2.0}$ | $18.0^{+1.8}_{-1.8}$ | $39.2^{+8.8}_{-5.8}$ | $17.4^{+1.4}_{-1.6}$ | $17.8$ | $18.2$
$M_{\mathrm{ini}}\,\,[M_{\odot}]$ | — | — | — | — | $40.2^{+9.1}_{-6.4}$ | $17.4^{+1.4}_{-1.6}$ | $33.0$ | $17.5$
$\log\dot{M}\,\,[{\,M_{\odot}\,\mbox{yr}^{-1}}]$ | $-6.1^{+0.2}_{-0.2}\,$ | $-8.8^{+0.5}_{-0.5}\,$ | — | — | $-6.1^{+0.3}_{-0.3}$ | $-7.9^{+0.3}_{-0.3}$ | $-6.4^{(f)}$ | $-6.46^{(f)}$
$\varv_{\infty}\,\,[$\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$]$ | $2500^{+200}_{-200}$ | $2500$ | — | — | — | — | — | —
$\varv\sin i\,\,[$\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$]$ | $\,\,\,140^{(d)}$ | $\,\,\,\,\,425^{(d)}$ | — | — | $150^{+26}_{-28}$ | $410^{+60}_{-40}$ | $96$ | $575$
$X_{\mathrm{H}}\,\,(\mathrm{by\,\,mass})$ | $0.73$ | $0.73$ | — | — | $0.74^{+0.0}_{-0.1}$ | $0.74^{+0.0}_{-0.1}$ | $0.47$ | $0.74$
$X_{\mathrm{C}}/10^{-5}\,\,(\mathrm{by\,\,mass})$ | $\,\,\,2^{+2}_{-1}$ | $\,\,\,21^{(e)}$ | — | — | $27^{+2}_{-5}$ | $7^{+5}_{-3}$ | $3$ | $20$
$X_{\mathrm{N}}/10^{-5}\,\,(\mathrm{by\,\,mass})$ | $\,\,\,45^{+5}_{-10}$ | $\,\,\,\,\,3^{(e)}$ | — | — | $\,\,\,6^{+14}_{-6}$ | $61^{+40}_{-30}$ | $91$ | $7$
$X_{\mathrm{O}}/10^{-5}\,\,(\mathrm{by\,\,mass})$ | $\,\,\,80^{+10}_{-20}$ | $\,110^{(e)}$ | — | — | $91^{+4}_{-14}$ | $70^{+7}_{-37}$ | $35$ | $110$
$\log\,Q_{H\textsc{\,}}\,[$\,\mathrm{s}^{-1}$]$ | $49.43$ | $47.88$ | — | — | — | — | $49.34^{(g)}$ | $47.72^{(g)}$
$\log\,Q_{{He\textsc{i}}{}}\,\,[$\,\mathrm{s}^{-1}$]$ | $48.72$ | $45.83$ | — | — | — | — | $48.62^{(g)}$ | $45.76^{(g)}$
$\log\,Q_{{He\textsc{ii}}{}}\,\,[$\,\mathrm{s}^{-1}$]$ | $40.80$ | $41.38$ | — | — | — | — | $43.63^{(g)}$ | $35.75^{(g)}$
$\mathrm{age}\,\,[$\mathrm{M}\mathrm{y}\mathrm{r}$]$ | — | — | — | — | $2.8^{+0.5}_{-0.4}$ | $5.9^{+1.5}_{-1.5}$ | $6.0$ | $6.0$
|  |  |  |  |  |  |  |

(a) Results from the PHOEBE code, see Sect. 3.2. (b) Results from the BONNSAI
tool, see Sect. 5.1. (c) Models calculated with MESA, see Sect. 5.2. No
uncertainties are given, as we did no in-depth analysis. (d) Obtained with the
`iacob` broad tool. (e) Initial CNO abundances scaled to SMC metallicity; for
more details see Sect. 4.2. (f) According to the mass-loss recipe implemented
in our evolutionary models. (g) Values taken from the spectroscopic model that
is calculated with the parameters of the binary evolutionary models.

The best fitting model light curve and RV fit obtained with the PHOEBE code
are shown in Fig. 2. The corresponding orbital parameters are listed in Table
5 and the stellar parameters of both components in Table 5. The orbital
solution yields similar masses for both components, while the light curve fit
reveals different $T_{\mathrm{eff}}$ and $L$. This indicates that there was a
prior mass-transfer phase that has stripped away most of the primary’s
envelope, such that now has similar mass to its companion. From the light-
curve fit, the fundamental stellar parameters –stellar radius and the surface
gravity— are determined, giving us the opportunity to cross-check the spectral
analysis. Furthermore, the PHOEBE code calculates the light ratio of the
binary components in the observed band outside conjunctions, which is compared
to the light ratio obtained from the spectroscopic analysis in Sect. 4.

##  4 Spectral analysis

###  4.1 Method: Spectral modeling with PoWR

Synthetic spectra for both stellar components were calculated with the Potsdam
Wolf-Rayet (PoWR) model atmosphere code (Gräfener et al. 2002; Hamann &
Gräfener 2004). In the following, we briefly describe the code. For further
details, we refer to Gräfener et al. (2002), Hamann & Gräfener (2003), Todt et
al. (2015), and Sander et al. (2015).

The PoWR code models stellar atmospheres and winds permitting departures from
local thermodynamic equilibrium (non-LTE). The code has been widely applied to
hot stars at various metallicities (e.g., Hainich et al. 2014, 2015; Oskinova
et al. 2011; Reindl et al. 2014; Shenar et al. 2015). The models are assumed
to be spherically symmetric, stationary, and in radiative equilibrium.
Equations of statistical equilibrium are solved in turn with the radiative
transfer in the co-moving frame. Consistency is achieved iteratively using the
“accelerated lambda operator” technique. This yields the population numbers
within the photosphere and wind.

The emergent spectrum in the observer’s frame is calculated with the formal
integral in which the Doppler velocity is split into the depth-dependent
thermal velocity and a “microturbulence velocity” $\zeta(r)$. The
microturbulence grows from its photospheric value
$\zeta_{\mathrm{ph}}=$10\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$$
linearly with the wind velocity up to $0.1\varv_{\infty}$.

By comparison of the synthetic spectrum with the observations, it is possible
to determine the main stellar parameters. In addition to the chemical
composition, the main stellar parameters that specify the model atmosphere are
the luminosity $L$, stellar temperature $T_{*}$, surface gravity $g_{*}$, wind
mass-loss rate $\dot{M}$, and the wind terminal wind velocity
$\varv_{\infty}$. The stellar temperature is defined as the effective
temperature referring to the stellar radius $R_{*}$ by the Stefan-Boltzmann
equation $L=4\pi\sigma R_{*}^{2}T_{*}^{4}$. The stellar radius is defined at
the Rosseland mean optical depth of $\tau=20$. The differences between $T_{*}$
and the effective temperature $T_{\mathrm{eff}}$ (referring to the radius
where the optical depth $\tau=2/3$) are negligibly small, as the winds of OB-
type stars are optically thin.

In the subsonic regions of the stellar atmosphere, the velocity law is
calculated such that the density, related by the equation of continuity,
approaches the hydrostatic stratification. The hydrostatic equation
consistently accounts for the radiation pressure (Sander et al. 2015). In the
supersonic regions, it is assumed that the wind velocity field can be
described by a so-called double $\beta$-law as was first introduced by Hillier
& Miller (1999):

$\varv(r)=\varv_{\infty}\left[\left(1-f\right)\left(1-\dfrac{r_{0}}{r}\right)^{\beta_{1}}+f\,\left(1-\dfrac{r_{1}}{r}\right)^{\beta_{2}}\right].$
(4)

In this work, we assume $\beta_{1}=0.8$, a typical value for O-stars,
$\beta_{2}=4$, and a contribution of $f=0.4$ of the second $\beta$-term.
$r_{0}$ and $r_{1}$ are close to the stellar radius $R_{*}$ and are determined
such that the quasi-hydrostatic part and the wind are smoothly connected. This
choice of the wind velocity law leads to better agreement between the
$\mathrm{H\alpha}$ absorption line and the C iv P Cygni line profile than the
classical $\beta$-law (Castor et al. 1975).

Inhomogeneities within the wind are accounted for as optically thin clumps
(“microclumping”) which are specified by the “clumping factor” $D,$ which
describes by how much the density within the clumps is enhanced compared to a
homogeneous wind with the same mass-loss rate (Hamann & Koesterke 1998b). In
our analysis, we use a depth-dependent clumping which starts at the sonic
point and increases outward until a clumping factor of $D=20$ is reached at a
radius of $R_{\mathrm{D}}=7\,R_{*}$ for the primary and at a radius
$R_{\mathrm{D}}=10\,R_{*}$ for the secondary. The smaller radius at which the
clumping factor is reached in the primary’s wind profile is needed to model
the observed O v $\lambda 1371$ line properly.

Figure 5: Selected regions of the X-SHOOTER spectrum (ID 5 in Table 1)
compared to different synthetic spectra. The observed spectrum, corrected for
the velocity of the SMC and the barycentric motion, is shown as a solid blue
line. The red dotted line is our best-fitting model which consists of the
primary with ${T_{\mathrm{eff,\,1}}=$42\text{\,}\mathrm{k}\mathrm{K}$}$ and
$\log g_{*,\,1}=3.7$ and the secondary with
${T_{\mathrm{eff,\,2}}=$32\text{\,}\mathrm{k}\mathrm{K}$}$ and $\log
g_{*,\,1}=4.0$. The gray solid line is again a combined synthetic spectrum of
the primary and secondary, but this time the surface gravity of the primary
has been reduced to $\log g_{*,\,1}=3.6$ (and the temperature and mass-loss
rate are slightly adjusted so that the spectrum matches the observations)
while the parameters of the secondary are kept fixed. The line identification
marks correspond to the wavelengths in the rest frame. Panel (a) shows the
$\mathrm{H\gamma}$ line, in which the red wing is dominated by the primary and
the blue wing by the secondary. Panel (b) shows the He ii $\lambda 5412$
absorption line, which is dominated by the primary and is sensitive to surface
gravity. Panel (c) shows the region of He ii $\lambda 6528$ and
$\mathrm{H\alpha}$. While the $\mathrm{H\alpha}$ wings are only barely
affected by the change in $\log g_{*,\,1}$, the He ii $\lambda 6528$ line is
more sensitive to it.

The PoWR model atmospheres used here account for detailed model atoms of H,
He, C, N, O, Mg, Si, P, and S. The iron group elements Sc, Ti, V, Cr, Mn, Fe,
Co, and Ni are combined to one generic element “G” with solar abundance
ratios, and treated in a superlevel approach (Gräfener et al. 2002). The
abundances of Si, Mg, and Fe are based on Hunter et al. (2007, their table 17)
and Trundle et al. (2007, their table 9). For the remaining elements, we
divide the solar abundances of Asplund et al. (2005, their table 1) by seven
to match the previously mentioned metallicity of the SMC. This yields the
following mass fractions: ${X_{\mathrm{H}}=0.73}$,
${X_{\mathrm{Si}}=$1.3\text{\times}{10}^{-4}$}$,
${X_{\mathrm{Mg}}=$9.9\text{\times}{10}^{-5}$}$,
${X_{\mathrm{P}}=$8.32\text{\times}{10}^{-7}$}$,
${X_{\mathrm{S}}=$4.42\text{\times}{10}^{-5}$}$, and
${X_{\mathrm{G}}=$3.52\text{\times}{10}^{-4}$}$. The CNO individual abundances
of both stellar components are not fixed but derived from the analysis. The
complement mass fraction to unity is $X_{\rm He}$.

We determined the color excess $E_{\mathrm{B-V}}$ and the luminosity $L$ of
the binary components by fitting the composite SED to photometry (top panel in
Fig. 9). Reddening is modeled as a combined effect of the Galactic foreground,
for which we adopt the reddening law of Seaton (1979) with
$E_{\mathrm{B-V}}=$0.03\text{\,}\mathrm{m}\mathrm{a}\mathrm{g}$$, and the
reddening law of Bouchet et al. (1985) for the SMC.

We use the `iacob-broad` tool (Simón-Díaz & Herrero 2014) in combination with
the high $\mathrm{S/N}$ X-SHOOTER spectrum ($\mathrm{S/N}\sim 100$) to
determine the rotation rates of the primary and the secondary. The helium
lines are potentially pressure broadened and thus are not optimally suited for
rotation broadening measurements. Therefore, for the primary, we use the few
metal lines that are visible in the optical spectrum, namely the N iv $\lambda
3478$ and the O iv $\lambda 3403$ absorption lines and the N iv $\lambda 4058$
emission line. We determine a projected rotation rate of ${\varv_{1}\sin
i=$140\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$}$ for the primary. The
secondary only contributes to the He i lines. Fitting the He i $\lambda 4387$
and He i $\lambda 4471$ absorption lines —in which the contribution of the
primary is smallest— yields a rotation rate of ${\varv_{2}\sin
i=$425\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$}$ for the secondary.

Figure 6: Part of the observed FUSE spectrum (ID 1 in Table 1) covering the
wavelength range of the LiF2 channel – the same as used by Penny & Gies (2009)
to determine the projected rotational velocity. The observed spectrum,
corrected for the velocity of the SMC, is shown as a solid blue line. The red
dotted line is our best-fitting model with the primary’s projected rotational
velocity ${\varv_{1}\sin
i=$140\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$}$. The gray solid line
is another model where the primary’s projected rotational velocity is reduced
to ${\varv_{1}\sin i=$140\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$}$.
The line identification marks correspond to the wavelengths in the rest frame.

Previously, Penny & Gies (2009) determined the projected rotational velocity
of the primary to ${\varv_{1}\sin
i=$65\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$}$ using a cross-
correlation of the FUSE spectrum in the far-UV with a template spectrum. We
calculate tailored spectral models with the respective projected rotational
velocities and find that also the FUSE spectrum is better reproduced when
using ${\varv_{1}\sin i=$140\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$}$
(see Fig. 6). We assume that the template spectrum used in Penny & Gies (2009)
might not have been perfectly calibrated for our target and that the binary
nature leads to additional uncertainties.

###  4.2 Resulting spectroscopic parameters

####  4.2.1 Temperature and surface gravity of the primary

AzV 476 was previously classified as O2-3 V plus a somewhat later O-type
companion. This implies that nitrogen lines are expected in absorption and
emission in the primary spectrum, while He i absorption lines should be
present in the secondary spectrum.

Indeed, we find that the N iv $\lambda 4057$ (hereafter N iv) emission line,
the N iv absorption lines at $\lambda\lambda\,\,$3463\text{\,}\AA$$,
$3478\text{\,}\AA$, $3483\text{\,}\AA\mathrm{,}$ and
$3485\text{\,}\AA\mathrm{,}$ and the He ii $\lambda 6528$ absorption can be
entirely assigned to the primary. We observe only marginal N iii
$\lambda\lambda 4634,4640$ (hereafter N iii) emission and N v $\lambda\lambda
4603,4619$ (hereafter N v) absorption; see Fig. 14. Thus, the optical spectrum
shows nitrogen only as N iv, while N iii as well as N v are virtually absent.
This restricts the temperature of the primary to a narrow range (Rivero
González et al. 2012). We also find that the O iv multiplets around
$3400\text{\,}\AA$ are purely associated with the primary. Because these lines
are highly contaminated by ISM absorption lines of Ti ii $\lambda 3385$ and Co
i $\lambda 3414,$ they are only used as a crosscheck of the oxygen abundance
applied in our spectral model. A more detailed description of this line
complex is given in Appendix B.2.

The secondary does not contribute noticeably to these weak metal lines.
However, the secondary strongly dominates the He i absorption lines at
$\lambda\lambda$ $3819\text{\,}\AA$, $4387\text{\,}\AA\mathrm{,}$ and
$4471\text{\,}\AA\mathrm{,}$ indicating a lower effective temperature. The
primary star also contributes to these lines, giving additional constraints on
the temperature of the primary.

The surface gravity $g_{*,\,1}$ of the primary is determined by fitting the
wings of the Balmer lines using those UVES spectra with the least wavy
patterns and highest RV shifts as well as the X-SHOOTER spectrum, which has a
higher $\mathrm{S/N}$. Because both stars contribute to the Balmer lines, the
He ii $\lambda 5412$ absorption line that is mostly originating from the
primary and also sensitive to changes in surface gravity is used.
Unfortunately, this line is only contained in the X-SHOOTER spectrum. Surface
gravity affects the density structure and therefore the ionization balance of,
for example, nitrogen and helium. As the temperature is adjusted such that the
nitrogen lines are reproduced, changes in the surface gravity will not only
affect the wings of the He ii lines but also change the depth of the He ii
absorption lines. This gives the opportunity to use the He ii $\lambda 6528$
absorption line in the UVES spectra as second criterion for the surface
gravity estimate in the primary.

We tested different temperature values in the range
${T_{\mathrm{eff,\,1}}=$3950\text{\,}\mathrm{k}\mathrm{K}$}$ and surface
gravities in the range
${\log(g_{*,\,1}/($\mathrm{c}\mathrm{m}\,\mathrm{s}^{-2}$))=$3.54.1$}$ and
found that a temperature of ${T_{\mathrm{eff,\,1}}=$42\pm
3\text{\,}\mathrm{k}\mathrm{K}$}$ and a surface gravity of
$\log(g_{*,\,1}/($\mathrm{c}\mathrm{m}\,\mathrm{s}^{-2}$))=3.7\pm 0.2$ are
most suitable for reproducing the primary’s spectrum. The given error margins
take into account that the measurements of temperature and surface gravity are
not independent.

Different regions of the X-SHOOTER spectrum that are used to determine the
surface gravity are shown in Fig. 5. The accuracy in surface gravity highly
depends on the calibration of the spectrum. One can see that
$\mathrm{H\gamma}$ and He ii $\lambda 5412$ line fits show preference to a
surface gravity for the primary of
$\log(g_{*,\,1}/($\mathrm{c}\mathrm{m}\,\mathrm{s}^{-2}$))=3.7$, while
$\mathrm{H\alpha}$ and He ii $\lambda 6528$ indicate that
$\log(g_{*,\,1}/($\mathrm{c}\mathrm{m}\,\mathrm{s}^{-2}$))=3.6$ might be more
suitable. However, from fitting the UVES data with the highest RV shifts we
find that a surface gravity for the primary of
$\log(g_{*,\,1}/($\mathrm{c}\mathrm{m}\,\mathrm{s}^{-2}$))=3.7$ yields the
best fit.

####  4.2.2 Temperature and surface gravity of the secondary

It is more difficult to determine the temperature of the secondary, as only a
few He i lines are associated with the secondary and all these line have a
contribution from the primary. However, from the spectra with the highest RV
shifts (e.g., UVES spectrum with ID 10), we find that the secondary does not
contribute to the He ii $\lambda 4200$ line, and therefore we have an
additional constraint that can be used as an upper limit on the temperature of
the secondary. In addition to the limited number of lines that are associated
with the secondary, the ionization balance that changes depending on the
surface gravity introduces another uncertainty in our temperature estimation.

Therefore, we first adjust $g_{*,\,2}$ of the secondary such that the wings of
the Balmer lines are well reproduced while keeping the surface gravity and
temperature of the primary fixed. As in the case of the primary, in order to
distinguish the contributions of the binary components to the Balmer wings, we
employ the UVES spectra in the phase with the largest RV shifts. This gives
additional information and a surface gravity of
${\log(g_{*,\,2}/($\mathrm{c}\mathrm{m}\,\mathrm{s}^{-2}$))=4.0}$ can be
determined. Finally, the temperature of the secondary is adjusted such that
the He i absorption lines at $\lambda\lambda$ $3819\text{\,}\AA$,
$4387\text{\,}\AA\mathrm{,}$ and $4471\text{\,}\AA$ match the observation.
Following this procedure, the spectrum of the secondary is best reproduced
with a temperature of
${T_{\mathrm{eff,\,2}}=$32\text{\,}\mathrm{k}\mathrm{K}$}$. This method is
accurate to $\Delta T_{\mathrm{eff,\,2}}=$\pm 4\text{\,}\mathrm{k}\mathrm{K}$$
in temperature and $\Delta\log(g_{*,\,2})=\pm 0.2$ for the surface gravity.
The obtained values of the surface gravity of the primary and the secondary
are in agreement with those obtained from the orbital analysis (see Sect. 3).

####  4.2.3 Luminosity

To estimate the light ratio, the luminosities of both stars are adjusted such
that the shapes of the optical He ii lines match the observations. This is an
iterative process done simultaneously with the temperature and surface gravity
estimate. The total luminosity of the system is calibrated such that the
calculated visual magnitude of the synthetic flux matches the observed one.
This results in luminosities of $\log(L_{1}/L_{\odot})=5.65$ for the primary
and $\log(L_{2}/L_{\odot})=4.75$ for the secondary. This estimate is sensitive
to synthetic spectra, light ratio, and distance modulus. Therefore, we
estimate that our measurements of the luminosity of each stellar component are
accurate to $\Delta\log(L/L_{\odot})=\pm 0.2$.

Taking these uncertainties into account, the model flux ratio in the OGLE
I-band yields $f_{\mathrm{2}}/f_{\mathrm{1}}=0.26\pm 0.1$ which is in
agreement with the flux ratio derived using PHOEBE
${f_{\mathrm{2}}/f_{\mathrm{1}}=0.31\pm 0.05}$ (see Sect. 3).

Figure 7: Observed (blue) and synthetic (red) C iv resonance doublet of our
final model. The individual primary and secondary spectra are shown as dotted
black and green lines, respectively. The sharp interstellar absorptions arise
in the Galactic foreground and in the SMC, and are also modeled with their
respective RV shifts. (a) COS observation (ID 2 in Table 1) (b) STIS
observation (ID 3 in Table 1).  Figure 8: Observed X-SHOOTER spectrum (ID 5 in
Table 1) and synthetic spectra of selected regions that show nitrogen emission
lines. The observed spectrum is shown as a solid blue line, while the dotted
red line is our best-fitting model. The observed spectrum is corrected for the
velocity of the SMC and the barycentric motion. The line identification marks
correspond to the wavelengths in the rest frame. Panel (a) shows the area
around the N iv $\lambda 4057$. Panel (b) shows the region of N iii
$\lambda\lambda 4511-4523$, N v $\lambda\lambda 4604,4620,$ and N iii
$\lambda\lambda 4634,4641$. The N iii and N v lines are not observable in the
spectrum. Panel (c) shows the region of the N iv $\lambda\lambda 7103-7129$
complex.

####  4.2.4 Mass-loss rates

After determining the temperatures, surface gravities, and luminosities of
both binary components, we proceed to measuring the properties of their
stellar winds. The $\mathrm{H\alpha}$ line in the optical as well as the C iv
resonance line in the UV are used as the main diagnostic tools for measuring
the mass-loss rate of the primary. In addition, we also pay attention to the
appearance of the optical He i $\lambda 4686$ and the UV He ii $\lambda 1640$
line; see Fig 16. When calibrating the mass-loss rate of the primary, the best
fit is archived with a mass-loss rate of
$\log(\dot{M}_{1}/({\,M_{\odot}\,\mbox{yr}^{-1}}))=-6.1\pm 0.2$ and a terminal
wind velocity of
$\varv_{\infty,\,1}=$2500\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$$.

We have two spectra covering the C iv resonance line at different phases, the
COS and the STIS spectrum (ID 2 and 3 in Table 1), depicted in Fig. 7. The C
iv P Cygni line profile in both spectra does not show any direct indications
for the contribution from the secondary such as a double emission peak or a
step in the absorption trough. The C iv resonance line in the STIS spectrum
shows slightly weaker emission, which is most likely due to the sensitivity of
the different instruments or weak wind variability. We use the bottom part of
the P Cygni profile of C iv resonance line to confirm the continuum
contribution of the secondary in the UV and therefore the used light ratio
($L_{2}/L_{1}$).

Apparently the secondary does not contribute to the wind lines. Its mass-loss
rate can therefore only be limited to
${\log(\dot{M}_{2}/({\,M_{\odot}\,\mbox{yr}^{-1}}))\leq-8.8\pm 0.5}$. For
higher mass-loss rates, we find that the secondary would noticeably contribute
to the C iv resonance line in the UV. Due to the lack of indications of the
secondary’s wind we adopt the terminal velocity of the primary,
$\varv_{\infty,\,2}=$2500\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$$.

Figure 9: Spectral fit for AzV 476. Top panel: SED. Flux-calibrated
observations (blue lines) are the FUSE (1), COS (2), STIS (4), and X-SHOOTER
(5) spectra as listed in Table 1. The light blue squares are the photometric
UBVRIJHK data as listed in Table 2. The model SED composed of both stellar
components is shown as a dashed red line, while the individual SEDs of the
primary and the secondary are shown as dotted black and green lines,
respectively. Lower panels: Normalized FUSE (1), COS (2), and STIS (4)
spectra. The number in the upper left corner corresponds to the ID given to a
specific spectrum as listed in Table 1. The line styles are the same as in the
top panel. The synthetic spectra are calculated with the model parameters
compiled in Table 6 (“Spectroscopy” columns). Figure 10: Same as Fig. 9, but
now the normalized optical X-SHOOTER spectra (ID 5 in Table 1) are shown.

####  4.2.5 CNO surface abundances

As already explained in Sect. 4.1, the abundances of hydrogen and the iron
group elements are fixed, while the CNO abundances of each stellar component
are adjusted during spectral modeling. The nitrogen abundance is adjusted to
find good agreement with the optical N iv absorption lines at wavelengths
$3463\text{\,}\AA$, $3478\text{\,}\AA$, and $3483\text{\,}\AA$, and the N iv
$\lambda 4057$ and N iv ${\lambda\lambda 7103-7129}$ emission lines, as well
as the UV N v $\lambda\lambda 1238,1242$ resonance doublet. The best agreement
is achieved with a nitrogen mass fraction of
${X_{\mathrm{N}}=(45^{+5}_{-10})\times${10}^{-5}$}$. The nitrogen lines seen
in the optical are displayed in Fig. 8.

A carbon abundance of ${X_{\mathrm{C}}=(2^{+2}_{-1})\times${10}^{-5}$}$
reproduces the optical C iv $\lambda\lambda 5801,5812$ lines best while
maintaining the fit of the C iv resonance doublet, at least for mass-loss
rates that do not conflict with other wind features. For lower abundances, the
latter is no longer saturated. The oxygen abundance of the primary is
calibrated with the oxygen absorption lines in the UV and the optical O iv
multiplets around $3400\text{\,}\AA$, resulting in
${X_{\mathrm{O}}=(80^{+10}_{-20})\times${10}^{-5}$}$.

There are no strong CNO lines in the secondary spectrum, and therefore we
adopt initial CNO abundances scaled to the metallicity of the SMC
($Z_{\mathrm{SMC}}=1/7\,Z_{\odot}$). Alternatively, one could assume that the
surface abundances in the secondary are close to those in the primary as the
accreted material is polluting the surface of the secondary. We explored both
assumptions but could not find diagnostic lines that would allow us to pin
down the composition of the atmosphere and wind of the secondary. We are only
able to fix upper limits for the CNO abundances such that these elements would
not produce features that are inconsistent with the observations:
$X_{\mathrm{N}}\lesssim 50\times${10}^{-5}$$, $X_{\mathrm{C}}\lesssim
21\times${10}^{-5}$$ and $X_{\mathrm{O}}\lesssim 110\times${10}^{-5}$$.

####  4.2.6 Required X-ray flux

Similar to other O-type stars, AzV 476 shows a notoriously strong P Cygni
profile in the O vi $\lambda\lambda 1032,\,1038$ resonance doublet. However,
wind models do not predict a sufficient population of this high ion without
the inclusion of additional physical processes. This phenomenon, termed
“super-ionization”, was described by Cassinelli & Olson (1979) and interpreted
as evidence for the presence of an X-ray field in stellar atmospheres.

In order to model the O vi doublet in the observed spectrum of AzV 476, we add
a hot plasma component. The X-ray-emitting plasma has an adopted temperature
of $T_{\mathrm{X}}=$3\text{\,}\mathrm{M}\mathrm{K}$$ and is distributed
throughout the wind outside a radius of $R_{\mathrm{X}}=1.1\,R_{*}$. Its
constant filling factor is a free parameter; the best reproduction of the
observation is achieved for a model with an emergent X-ray luminosity of
${\log\,L_{\mathrm{X}}=$31.4\text{\,}\mathrm{e}\mathrm{r}\mathrm{g}\,\mathrm{s}^{-1}$}$.

Current X-ray telescopes are not sensitive enough to detect individual O stars
in the SMC. Our final model predicts an X-ray flux at earth of
$1.4\text{\times}{10}^{-17}\text{\,}\mathrm{e}\mathrm{r}\mathrm{g}\,\mathrm{s}^{-1}\,\mathrm{c}\mathrm{m}^{-2}\,\AA^{-1}$
integrated over the X-ray band ($\approx 6-60$ Å). The region on the sky where
AzV 476 is located was observed by both XMM-Newton and Chandra X-ray
observatories. However, despite the detection of AzV 476 in the UV by the XMM-
Newton optical monitor, the star is not detected in X-rays. Hence, we set the
upper limit on its X-ray flux at the median flux of detected sources in the
XMMSSC - XMM-Newton Serendipitous Source Catalog (4XMM-DR10 Version); this is
$$5\text{\times}{10}^{-15}\text{\,}\mathrm{e}\mathrm{r}\mathrm{g}\,\mathrm{s}^{-1}\,\mathrm{c}\mathrm{m}^{-2}$\,\AA^{-1}$
, which is far from being sufficiently sensitive for our predicted source.

The O vi resonance doublet is not the only feature that is sensitive to
X-rays. Other ions which might be populated by Auger ionization are N v and N
iv. We find that the N v $\lambda\lambda 1242,\,1238$ resonance doublet in our
model is not significantly affected by the applied X-ray field, while the N iv
$\lambda\lambda 1718,\,1721$ doublet is very sensitive. However, in the region
around the N iv $\lambda\lambda 1718,\,1721$ doublet, the observations are
very noisy and can be explained by the model with X-rays as well as without.

The composed model spectrum as well as the individual spectra are shown in
Figs. 5 – 10. A summary of the stellar parameters including the abundances is
given in Table 6.

####  4.2.7 Spectroscopic stellar masses

The spectroscopic masses derived from our analysis are
${M_{\mathrm{spec,\,1}}=29^{+17}_{-11}\,M_{\odot}}$ for the primary and
$M_{\mathrm{spec,\,2}}=22^{+13}_{-8}\,M_{\odot}$ for the secondary. The
spectroscopic masses, although a factor of $1.5$ higher than the orbital
masses, agree with these latter within their respective uncertainties. The
large error margins on the spectroscopic masses arise mainly from the large
uncertainties on surface gravity and luminosity (e.g., Fig. 5). The
spectroscopic mass ratio is ${q_{\mathrm{spec}}=0.7}$ which is smaller than
the one obtained from the orbital analysis.

It appears that the spectroscopically derived masses and luminosities are
shifted systematically to higher values compared to the results of the orbital
analysis. We discuss this in more detail in Sect. 6.1.

Nevertheless, a mass ratio close to unity, while the two stars in a close
binary have different spectral types, strongly suggests mass transfer in the
past. From our binary evolutionary models (Sect. 5), we expect that mass
transfer removed most of the hydrogen-rich envelope from the primary. We
calculated model spectra with strongly depleted surface hydrogen
($X_{\mathrm{H}}=0.2$ and $0.5$). However, these models yield poorer spectral
fits as the helium lines become too deep. Still, moderate hydrogen depletion
and helium enrichment cannot fully be excluded.

###  4.3 Revisiting the spectral classification

####  4.3.1 Spectral classification aided by stellar atmosphere modeling

AzV 476 was previously classified on the basis of the appearance of its
optical spectrum. However, during our spectroscopic analysis, we found that
most of the classification lines are blended by the companion. Hence,
including the additional information from our spectroscopic models and the
available UV data allows us to reconsider the spectral classification.

We used the Marxist Ghost Buster (MGB) spectral classification code (Maíz
Apellániz et al. 2012; Maíz Apellániz 2019). MGB compares observed spectra
with the criteria of Walborn et al. (2000, 2002, 2014) and Sota et al. (2011,
2014) while allowing the presence of two components by iteratively adjusting
the spectral types, light ratio, RVs, and rotational indices. This approach
poses two problems. First, the low metallicity of the SMC leads to weaker
nitrogen lines compared to stars with similar type in the MW and LMC. Second,
only the X-SHOOTER spectrum is covering the whole wavelength range of
$39005100\text{\,}\AA$ which is needed for a spectral classification. In
addition we consider those UVES spectra with the highest RV shifts in order to
better separate the individual binary components.

Therefore, we are forced to use the X-SHOOTER spectrum first and then consider
the UVES spectra with the highest RV shifts that cover the range from 3300 to
4500Å.

The MGB code provides the spectral type O4 IV((f)) for the primary. The main
reason for the slightly later subtype assigned to the primary is its
contribution to the He i $\lambda 4471$ line. The suffix “((f))” is assigned
because of the barely detectable N iii $\lambda\lambda 4634,4641$ lines (see
Fig. 14). Regarding the secondary, the absence of the He ii $\lambda 4200$
line suggests the latest O subtype, O9.5: Vn. The suffix “n” refers to the
high rotational broadening of the lines of the secondary.

The MGB classification code only considers the optical range of the spectrum.
If one inspects the UV spectrum, which mostly originates from the primary, and
compares it to the N v and C iv UV wind-profile templates with those presented
in Walborn et al. (2000), one might prefer a luminosity class “III” for the
primary. The difference in the optical spectral appearance is explainable by
the high L/M ratio which leads to very strong mass loss and thus to
exceptionally strong wind lines in the UV. In light of these circumstances and
the unusual chemical composition, we add the suffix “p” for this peculiar
object. In the end, we therefore classify the system as O4 IV-III((f))p +
O9.5: Vn.

Following the work of Weidner & Vink (2010), an isolated star with spectral
type O4 III is expected to have a stellar mass of about
$M_{\mathrm{exp,\,1}}=49^{+7}_{-6}\,M_{\odot},$ which is more than twice as
high as that which we derive using two different methods (Sects. 3 and 4.1).
For the O9.5: V secondary, the mass expected from the spectral type would be
$M_{\mathrm{exp,\,2}}=16^{+7}_{-3}\,M_{\odot},$ which agrees within
uncertainties with the findings of the orbital and spectroscopic analysis.

##  5 Stellar evolution modeling

From orbital and spectroscopic analysis, we find that the mass ratio is close
to unity, although the two binary components have distinctly different
luminosities, effective temperatures, and rotation rates: the secondary is a
fast rotator while the rotation of the primary is only moderate. These facts
strongly indicate that the system has already undergone mass transfer.

###  5.1 Single-star models

As a first approach we investigate whether the stellar parameters of the
individual components can be reproduced by single-star evolutionary models.
For this purpose, we employ the Bayesian statistic tool “The BONN Stellar
Astrophysics Interface” (BONNSAI444www.astro.uni-bonn.de/stars/bonnsai),
(Schneider et al. 2014) in combination with the BONN-SMC tracks (Brott et al.
2011). To find a suitable model, we request the tool to match the current
luminosity, effective temperature, rotation rate, and orbital mass of the two
stars. Indeed, we find that the fast-rotating secondary can be partially
explained by a single-star track. With respect to the primary, the BONNSAI
tool cannot find any track that would explain its current low mass. Only when
ignoring the mass constraint are we able to find a suitable track that
reproduces the remaining stellar parameters. In that case, the predicted ages
of the primary and secondary differ by a factor of two. This finding confirms
our suggestion that the system has already undergone mass transfer. The best-
fitting stellar parameters obtained with the BONNSAI tool —applicable only for
single stars— are listed in Table 6.

###  5.2 Binary evolutionary models

Figure 11: HRD showing the positions of each binary component according to the
spectroscopic and orbital analysis, and the best fitting binary evolutionary
tracks. The results of the spectroscopic analysis are shown as triangles and
the orbital analysis as upside-down triangles. The positions of the primary
and secondary are marked in green and red, respectively. The shaded areas
indicate the respective error-ellipses. Evolutionary tracks (solid lines) of
both binary components are according to the best fitting binary model
calculated with MESA (see Sect. 5.2). The tracks are labeled by their initial
masses. The black dots on the tracks correspond to equidistant time-steps of
$0.1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ to emphasize the most probable
observable phases. The light blue dots mark the best-fitting model.

The binary models are calculated with MESA v. 10398. To mimic the BONN-SMC
tracks we adjusted MESA in a similar way to that described by Marchant (2016).
We follow most of the physical assumptions from the BONN-SMC models (e.g.,
overshooting, thermohaline mixing, etc.) and adopt them from Marchant (2016)
with two exceptions. In our models, we use a more efficient semi-convection
with $\alpha_{\mathrm{sc}}=10$ and calculate the mass transfer according to
the “Kolb” scheme (Kolb & Ritter 1990), as it allows us to include the
eccentricity enhancement mechanism meaning that the evolution of the
eccentricity can be modeled properly.

We want to emphasize that our goal is to check whether binary evolutionary
models can explain the stellar parameters of AzV 476 and especially its
masses. Therefore, we only explore a narrow parameter space and compute a
small set of models with initial primary masses in the range of
${M_{\mathrm{ini,\,1}}=($2538$)\,M_{\odot}}$ and secondary masses
${M_{\mathrm{ini,\,2}}=($1025$)\,M_{\odot}}$, orbital periods in the range of
${P_{\mathrm{\,ini}}=$825\text{\,}\mathrm{d}$}$ and initial eccentricities
${e_{\mathrm{\,ini}}=$0.00.4\text{\,}$}$. The initial parameters are adjusted
such that at some later stage the model binary has properties similar to AzV
476. We assume that the stellar components initially rotate with
${\varv_{\mathrm{rot}}=0.1\,\varv_{\mathrm{crit}}\approx$65\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$}$.
These values are chosen arbitrarily but as we are dealing with a close binary,
the tidal forces nevertheless lead to a tidal synchronization. Further
important assumptions are as follows:

Our models include the effect of inflation, which appears inside a stellar
envelope when the local Eddington limit is reached and exceeded, leading to a
convective region and a density inversion (Sanyal et al. 2015).

The stellar wind prescription is inspired by the work of Brott et al. (2011).
The winds of hot H-rich stars are described according to the Vink et al.
(2001) recipe. For stars with effective temperatures below the bi-stability
jump, where mass-loss rates abruptly increase (see Vink et al. 2001), we use
the maximum $\dot{M}$ from either Vink et al. (2001) or Nieuwenhuijzen & de
Jager (1990). WR mass-loss rates are according to Shenar et al. (2019) but are
assumed to scale with metallicity as $\dot{M}\propto Z^{1.2}$ as recommended
by Hainich et al. (2017). In the transition phase from a hot H-rich star
($X_{\mathrm{H}}\geq 0.7$) to the WR stage ($X_{\mathrm{H}}\leq 0.4$), the
mass-loss rate is interpolated between the prescriptions of Vink et al. (2001)
and Shenar et al. (2019).

Rotational mixing is modeled as a diffusive process including the effects of
dynamical and secular shear instabilities, the Goldreich-Schubert-Fricke
instability, and Eddington-Sweet circulations (Heger et al. 2000). In addition
to the angular momentum transport by rotation, the transport via magnetic
fields from the Tayler-Spruit dynamo (Spruit 2002) is also included.

Convection is described according to the Ledoux criterion and the mixing
length theory (Böhm-Vitense 1958) with a mixing length parameter of
${\alpha_{\mathrm{mlt}}=l/H_{\mathrm{P}}=1.5}$. For hydrogen burning cores, a
steep overshooting is used such that the convective core is extended by
${0.335H_{\mathrm{P}}}$ (Brott et al. 2011) where ${H_{\mathrm{P}}}$ is the
pressure scale height at the boundary of the convective core. Thermohaline
mixing is included with an efficiency parameter of ${\alpha_{\mathrm{th}}=1}$
(Kippenhahn et al. 1980), as well as semiconvection with an efficiency
parameter of ${\alpha_{\mathrm{sc}}=10}$ (Langer et al. 1983).

Mass transfer in a binary is modeled using the “Kolb” scheme (Kolb & Ritter
1990). This allows us to use the Soker eccentricity enhancement (Soker 2000),
which assumes phase-dependent mass loss and calculates the change in
eccentricity due to the mass that is lost from the system, and the mass that
is accreted by the companion. Tidal circularisation is taken into account
throughout the entire evolution of the system.

The remaining mass collapses into a compact object —represented by a point
mass— when helium is depleted in the stellar core. This simplification avoids
numerical problems in the latest burning stages and the unknowns that come
with a supernova explosion and a possible kick.

As the luminosities obtained from the spectroscopic and orbital analysis
differ to some extent, we put additional focus on reproducing the orbital
masses, which are our most reliable estimate. The best agreement with the
empirically derived parameters is found for a system with initial masses
${M_{\mathrm{ini,\,1}}=33\,M_{\odot}}$ and
$M_{\mathrm{ini,\,2}}=17.5\,M_{\odot}$, initial orbital period
${P_{\mathrm{\,ini}}=$12.4\text{\,}\mathrm{d}$}$, and initial eccentricity
${e_{\mathrm{\,ini}}=0.14}$. The model closest to the current stellar
parameters gives similar masses for the primary and secondary:
${M_{\mathrm{evol,\,1}}=17.8\,M_{\odot}}$ and
${M_{\mathrm{evol,\,2}}=18.2\,M_{\odot}}$,. The remaining stellar parameters
of our favorite models are listed in Table 6.

The Hertzsprung-Russell diagram (HRD) displaying the best fitting stellar
evolutionary tracks for both binary components is illustrated in Fig. 11. The
binary evolution model is able to reproduce almost all empirically derived
stellar parameters including the current orbital period
(${P_{\mathrm{model}}=$9.3\text{\,}\mathrm{d}$}$) and the eccentricity of the
system (${e_{\mathrm{model}}=0.25}$). However, the evolutionary models over-
predict the rotation rate of the secondary and the surface abundances of the
primary; see Table 6. We calculated PoWR atmosphere models with the best
fitting parameters obtained with our MESA models. The synthetic spectra are
shown in Fig. 17.

A more fine-tuned model would likely solve some of these problems. For
example, we did not include the initial rotation as a free parameter. This
might solve the faster rotation of the primary; however, the secondary would
still be expected to be rotating too fast as it spins up to criticality during
mass transfer. According to our favorite evolutionary model, the primary in
AzV 476 is currently evolving towards the helium zero age MS (ZAMS) and will
possibly spend the rest of its life as a hot helium or WR-type star.

##  6 Discussion

###  6.1 Comparison of the orbital, spectroscopic, and evolutionary mass
estimates

The orbital masses and the masses predicted by the binary evolutionary models
agree well within their error margins and deviate only by $10\%$. On the other
hand, the spectroscopic mass and luminosity of both components appear to be
higher by a factor of $1.5$ compared to the orbital solutions. The question
therefore arises as to whether the discrepancy between the orbital and
spectroscopic solutions is significant.

####  6.1.1 Orbital versus spectroscopic mass

The simplest way to resolve this discrepancy is to assume that the distance to
the system is lower than the canonical SMC distance of
(${d=$55\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$}$). As the luminosity depends
on the assumed distance as $L\propto d^{-2}$, the radius
(${R\propto\sqrt{L}\propto d^{-1}}$) and the spectroscopic mass also depend on
distance (${M_{\mathrm{spec}}\propto R^{2}\propto L\propto d^{-2}}$). The SMC
galaxy is extended, and the distances to its various structural parts are not
well established. Tatton et al. (2021) suggest that young structures, such as
the NGC 456 cluster, are not well represented by the intermediate age stellar
populations, which are usually used for distance estimates. Instead, young
clusters might be located in front of these latter. This is in agreement with
the argumentation by Hammer et al. (2015), who showed that the interactions
between the LMC and SMC could lead to multiple tidal and ram-pressure stripped
structures, which spatially separate young star-forming regions from older
stellar populations.

In order to bring the spectroscopic and orbital parameters of AzV 476 into
accordance, the system needs to be shifted to a distance of
$49\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$. The resulting spectroscopic
masses are then ${M_{\mathrm{shift,\,1}}}\approx 23\,M_{\odot}$ and
${M_{\mathrm{shift,\,2}}\approx 18\,M_{\odot}}$. However,
$d=$49\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$$ would imply that our target is
located at the same distance as the LMC.

The LMC has higher metallcity (${Z_{\mathrm{LMC}}=1/2\,Z_{\odot}}$) than the
SMC (${Z_{\mathrm{SMC}}=1/7\,Z_{\odot}}$), and therefore we speculate that
stars formed in the interaction regions may have higher metallicity than the
SMC. To test this, we increased the metallicity content of our spectroscopic
models to the LMC metallicity while keeping the remaining stellar parameters
unchanged. We find that the lines in the iron forest, especially in the far-
UV, are too deep compared to the observations. Therefore, it cannot be
confirmed that our target has LMC metallicity, but at the same time this does
not provide additional constraints on its distance.

Another option we need to consider is the possibility of a third light
contribution. AzV 476 is located in a very crowded region. Therefore, it is
possible that the observed light is contaminated by a nearby object and hence
the true luminosity (and mass) is lower. We inspected the HST acquisition
image of our target and find no other nearby UV-bright stars. Alternatively
our target could be a multiple system rather than a binary. If there were a
third component, it would need to be of similar luminosity to the secondary in
order to have an impact on our analysis. However, we do not see any
indications of a third component in our spectra.

A third option is that the surface gravities determined from the spectral
analysis are overestimated. The accuracy of the spectroscopic mass strongly
depends on surface gravity, $\log\,g$, which is mainly determined from fitting
the pressure-broadened line wings. Such estimates suffer from limited $S/N$
ratio in the spectra. In cases where the Balmer lines include contributions
from both components, further uncertainties are induced. As explained in Sect.
4.2 and shown in Fig. 5, the He ii lines help to improve the estimates, but
still do not allow us to pin point $\log\,g$ for each component separately.
This leaves room for various combinations of $\log\,g$ able to reproduce the
Balmer lines and is reflected in the large uncertainties of the spectroscopic
masses.

####  6.1.2 Orbital versus evolutionary mass

The evolutionary modeling yields a slightly lower mass for the primary than
the spectroscopic and orbital analyses. As can be seen in Table 6, the
evolutionary mass is $\leavevmode\nobreak\ 2\,M_{\odot}$ lower than the
orbital one.

For a given initial mass, the main parameter defining how much mass is lost
during mass transfer by the primary is the initial binary orbital period. As a
rule of thumb, the longer the initial orbital period, the less mass is removed
from the donor star (Marchant 2016). However, several effects (e.g.,
eccentricity) influence the response of the orbit to mass transfer. For
instance, an eccentric orbit leads to phase-dependent mass transfer, which
potentially increases the eccentricity and widens the orbit. Thus, the donor
can expand more, and less mass is removed via Roche-lobe overflow (Soker 2000;
Vos et al. 2015).

The amount of mass lost during mass transfer is also reduced when the radial
expansion of the star is minimized such that the star can stay underneath the
Roche lobe. The radial expansion is affected by the mixing processes. For
instance, a star with more efficient semi-convective mixing is expected to be
more compact (e.g., Klencki et al. 2020; Gilkis et al. 2021). We computed a
binary evolutionary model with a less efficient semi-convection of
${\alpha_{\mathrm{sc}}=1}$ and find that the post-mass-transfer (after
detachment) mass of the primary is lowered by $1\,M_{\odot}$. This is only one
example and there are a handful of other processes that impact the mixing
efficiency, energy transport, and the expansion of the envelope.

One more process that affects the radial expansion is the mass-loss rate which
is described by the adopted mass-loss recipe (see Sect. 5). Over- or
underestimating the mass-loss rate influences stellar evolution before,
during, and after mass transfer. Interestingly, we find approximate agreement
between the empirically derived mass-loss rate of the primary and the
prediction of the wind recipe (see Sect. 6.2.1).

The evolutionary mass of a star after mass transfer is sensitive to various
assumptions. Therefore, a difference of $2\,M_{\odot}$ between evolutionary
and orbital mass appears to be within the range of the evolutionary model
uncertainties. While a more fine-tuned binary evolutionary model is beyond the
scope of this paper, our modeling confirms that the orbital masses of AzV 476
are well explained by a post-mass-transfer binary evolutionary model.
Therefore, this unique system is a useful benchmark for improving the
understanding of mass transfer and mixing processes in close massive binaries.
We conclude that the orbital, spectroscopic, and evolutionary mass estimates
agree within their uncertainties.

###  6.2 The empirical mass-loss rate

Stellar evolution depends on the mass-loss rate. Stellar evolution models rely
on recipes that are only valid for a given evolutionary phase. In transition-
phases one typically interpolates between corresponding mass-loss recipes
(e.g., Sect. 5). In the following, we compare the mass-loss rates we
empirically derive from spectroscopy with other estimates as well as with the
various recipes.

In this work, the mass-loss rates are primarily derived from resonance lines
in the UV and the $\mathrm{H\alpha}$ line while taking into account the
morphology of the He ii $\lambda 1640$ and He ii $\lambda 4686$ lines. The
winds of all hot stars are clumpy (e.g., Hamann et al. 2008) which enhances
the emission lines fed by recombination cascades (e.g., Hamann & Koesterke
1998a) such as the $\mathrm{H\alpha}$ line. Resonance lines in the UV are
mostly formed by line scattering, a process that scales linearly with density
and is thus independent from microclumping. Considering both
$\mathrm{H\alpha}$ and the resonance lines allows us to constrain the clumping
factor $D$ and the mass-loss rate $\dot{M}$ consistently.

####  6.2.1 The primary

The spectroscopically derived mass-loss rate of the primary is
${\log(\dot{M}/({\,M_{\odot}\,\mbox{yr}^{-1}}))=-6.1\pm 0.2}$. Massa et al.
(2017) derive a mass-loss rate of
$\log(\dot{M}/({\,M_{\odot}\,\mbox{yr}^{-1}}))=-5.54^{+0.12}_{-0.10}$ from the
IR excess. However, the clumping effects were not included in their analysis.
As the IR excess originates from free-free emission, it scales with
$\sqrt{D}$. Correcting the mass-loss rate with a clumping factor of $D=20$, as
assumed in our spectroscopic models, yields
$\log(\dot{M}/({\,M_{\odot}\,\mbox{yr}^{-1}}))=-6.19^{+0.12}_{-0.10}$, which
is in agreement with the optical and UV analyses.

In the massive star evolutionary models, the terminal wind velocity is
prescribed as $\varv_{\infty}/\varv_{\mathrm{esc}}=2.6$ and the mass-loss
rates of OB stars are prescribed according to the Vink et al. (2001) recipe.
Using the stellar parameters obtained with spectroscopic analysis (see Table
6) and the prescribed
$\varv_{\infty}=$1911\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$$, the
“Vink’s recipe” yields ${\log(\dot{M}/({\,M_{\odot}\,\mbox{yr}^{-1}}))=-6.1}$.
Alternatively, applying the actual wind velocity measured from the UV spectra,
$\varv_{\infty}=$2500\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$$, the
theoretical prediction changes to
${\log(\dot{M}/({\,M_{\odot}\,\mbox{yr}^{-1}}))=-6.2}$, which is in agreement
with our empirical result within the uncertainties.

However, the primary’s current mass-loss rate adopted in our favored binary
evolutionary model (${\log(\dot{M}/({\,M_{\odot}\,\mbox{yr}^{-1}}))=-6.4}$) is
lower. This is mainly caused by the interpolation between the mass-loss rate
prescriptions for OB stars by Vink et al. (2001) and for WR stars by Shenar et
al. (2019) which is employed when the surface hydrogen mass-fraction of the
primary’s evolutionary model already dropped below $X_{\mathrm{H}}<0.7$.
Additionally, the stellar parameters (e.g., mass and luminosity) that enter
the different mass-loss recipes employed by the evolutionary code are slightly
different than those obtained from our spectral analysis.

Based on dynamically consistent Monte Carlo wind models, Vink & Sander (2021)
suggest an updated scaling of metallicity for their mass-loss rates. This
scaling yields higher mass-loss rates of
${\log(\dot{M}/({\,M_{\odot}\,\mbox{yr}^{-1}}))=-5.7}$, or, using the measured
$\varv_{\infty}=$2500\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$$,
${\log(\dot{M}/({\,M_{\odot}\,\mbox{yr}^{-1}}))=-5.9}$. Correcting in a simple
way the former theoretical mass-loss rate for clumping results in
${\log(\dot{M}/({\,M_{\odot}\,\mbox{yr}^{-1}}))=-6.4}$, which is somewhat
lower than our empirical measurement.

Figure 12: Abundance profiles of the primary at different evolutionary stages
as predicted by the evolutionary models. The legend at the top of each plot
explains the colors indicating the dominating mixing processes. The H, He, C,
N, and O abundances of the model are shown as solid lines while the observed
surface CNO abundances are shown as dotted lines of the same respective color.
The vertical black dotted line indicates the mass of the primary after mass
transfer. Left panel: Primary has depleted hydrogen in the core and is at the
terminal age MS. Right panel: Primary fills its Roche lobe and starts
transferring mass.

The spectroscopic measurements of mass-loss rates of O-dwarfs in the SMC
deliver useful empiric prescriptions (e.g., Bouret et al. 2003; Ramachandran
et al. 2019). According to Ramachandran et al. (2019, see their figure 15),
the mass-loss rate of the primary star in AzV 476 would be only
${\log(\dot{M}/({\,M_{\odot}\,\mbox{yr}^{-1}}))=-7.4}$, which is much lower
than obtained from our analysis. We explain this discrepancy by the advanced
evolutionary status: the primary has already a substantially reduced mass due
to its close-binary evolution, leading to a high L/M ratio. This precludes the
comparison with the recipes which use the mass–luminosity relations for single
stars.

####  6.2.2 The secondary

There are no indications of wind lines of the secondary in the observed
spectra. Therefore, we can only set an upper limit on the mass-loss rate
${\log(\dot{M}/({\,M_{\odot}\,\mbox{yr}^{-1}}))\leq-8.8\pm 0.5}$. The Vink et
al. (2001) recipe predicts a higher value of
${\log(\dot{M}/({\,M_{\odot}\,\mbox{yr}^{-1}}))=-8.2}$. Correcting for
clumping with $D=20$ yields
${\log(\dot{M}/({\,M_{\odot}\,\mbox{yr}^{-1}}))=-8.8}$ which is in agreement
with the empirical result.

According to Vink & Sander (2021), the predicted mass-loss rate of the
secondary is ${\log(\dot{M}/({\,M_{\odot}\,\mbox{yr}^{-1}}))=-7.8}$;
correcting it for clumping yields
${\log(\dot{M}/({\,M_{\odot}\,\mbox{yr}^{-1}}))=-8.4}$, which is higher than
the empirical upper limit. The empirical relation of Ramachandran et al.
(2019, their figure 15,) suggests
${\log(\dot{M}/({\,M_{\odot}\,\mbox{yr}^{-1}}))=-8.6\pm 0.2}$, which is
somewhat higher than our result but consistent with its error margins.

As can be seen in Table 6, the mass-loss rate used in the binary evolution
calculation is higher than predicted by the recipes as well as empirically
measured. This is because in the MESA models, the mass-loss rate of a star
rotating close to critical is enhanced compared to slowly rotating stars
(Paxton et al. 2013). The synthetic spectrum calculated with the parameters
obtained by the binary evolutionary models displayed in Fig. 17 shows strong P
Cygni profiles in the spectrum of the secondary (e.g., Si iv $\lambda\lambda
1393.8,1402.8$) that are not observed. We conclude that the mass loss of the
secondary is not rotationally enhanced above the observed limit.

###  6.3 Comparison of the stellar parameters as predicted by MESA
evolutionary tracks versus derived spectroscopically

In this work, AzV 476 is studied using two different approaches. (1) The
empiric approach is to model the observed light curves and spectra using the
PHOEBE code and the stellar atmosphere model PoWR in order to derive orbital,
stellar, and wind parameters. (2) The evolutionary modeling approach is to use
the MESA code to compute tracks reproducing the current stellar parameters of
AzV 476 components (Fig. 11) simultaneously with its orbital parameters. Below
we compare the outcomes of these two approaches.

In binary evolutionary models, during the mass-transfer phase, most of the
hydrogen-rich envelope is stripped off such that products of the CNO burning
process are exposed at the surface. The predicted surface abundances in our
favorite model correspond to an intermediate stage between the initial CNO
abundances and the CNO equilibrium at metallicity
${Z_{\mathrm{SMC}}=1/7\,Z_{\odot}}$. This means that most of the initial C, N,
and O become N, and the CO abundances are depleted.

Indeed, the abundances derived spectroscopically are intermediate between the
initial CNO abundances and the CNO equilibrium. The C and O abundances are
reduced while the N abundance increased drastically. However, we find a factor
of two difference between the predicted and observed N and O abundances (see
Table 1). We note that the _surface_ CNO abundances predicted by the MESA
models are subject to various assumptions; for example, on mixing efficiencies
or the mass removal by Roche-lobe overflow. Because of the scarcity of
suitable photospheric absorption lines, the CNO abundances of the primary are
partially derived from emission lines. Those lines are sensitive to many
effects, such as the mass-loss rate, temperature, clumping, and turbulence,
introducing additional uncertainties in the abundance measurement.

Furthermore, we find disagreement between the observed and predicted H mass-
fraction. In evolutionary models, a significant amount of the H-rich envelope
(i.e., with $X_{\rm H}\,>\,0.7$) is removed during the mass-transfer phase,
revealing the products of the established chemical gradient between the core
and the envelope, leading to the depletion of surface hydrogen. However, from
the spectroscopic analysis, we do not find such strong H depletion. As in the
case of the CNO measurements, this indicates some deficiencies in the
evolutionary models that are likely related to the mass-transfer phase and/or
to the mixing.

For a better understanding, we are showing in Fig. 12 abundance profiles at
two different stages. The left panel in Fig. 12 shows the abundances at core
hydrogen depletion. At this stage, the star is hydrogen shell burning, which
leads to the formation of an intermediate convection zone (ICZ, Langer
(1987)), which ranges from the mass coordinates
$2124\text{\,}\mathrm{M}_{\odot}$. One can see that the H, He, O, and N
abundances in this region are the same as the abundances that we observe; only
the predicted C abundance is higher. However, the spectroscopically derived C
abundance is the least accurate one, as the observed spectra have no
photospheric C absorption lines and the C iv resonance line in the UV highly
depends on the mass-loss rate. At this evolutionary stage, the envelope is not
well mixed, which can be seen in the step-like structure between the different
mixing regions.

The right panel in Fig. 12 shows the abundances at the onset of mass transfer.
Mass transfer happens on the dynamical timescale, while the mixing and burning
processes take place on longer timescales. Therefore, this plot can be used to
decipher the surface composition of the final model if a specific amount of
mass is removed from the star. Compared to the left panel, the stellar model
has had more time to mix the CNO elements through the envelope. The ICZ is now
largely extended and the CNO abundance ratios have changed; for instance N is
now the most abundant element in this region. We would like to emphasize that
this highly depends on the assumed mixing efficiencies, as a less efficient
semi-convection (e.g., $\alpha_{\mathrm{sc}}=1$) would lead to less mixing and
might preserve the step-like structure seen the left panel. Semi-convection
does not only affect the mixing efficiency but also impacts the change of the
stellar radius (Klencki et al. 2020), and therefore a star with more effective
semi-convective mixing spends more time hydrogen shell burning before it fills
the Roche lobe and initiates mass transfer. This additionally impacts the
final chemical composition.

By now it should be evident that accurate prediction of the surface CNO
abundances in binary evolutionary models is a nontrivial task and is beyond
the scope of this paper. However, this underlines the importance of AzV 476,
which can be used to get a better understanding of the different mixing
efficiencies, as it reveals the abundances formed by burning stages that would
normally be hidden from the observer, making them difficult to calibrate.

In order to understand how the predicted stellar parameters of the binary
evolution models — including the surface abundances — compare to the observed
spectrum, synthetic spectra of the primary and secondary are calculated using
these parameters (see Table 6). The resulting spectrum is shown in Fig. 17.
When comparing the synthetic spectrum to the observed one, it is evident that
the N emission lines are too strong and the helium absorption lines are too
deep.

The abundances are not the only difference between the evolutionary model
predictions and the observations. Evolutionary model predicts much higher
rotation rates than observed. In the evolutionary model, the secondary has
spun up to critical rotation, which is not observed, and also the rotation
rate of the primary is somewhat slower than observed. Vanbeveren et al. (2018)
and Shara et al. (2020) studied the rotation rates of O stars in WR+O star
binaries at different metallicities and compared them to predictions from
evolutionary models. This comparison revealed that the observed rotation rates
of the O stars are lower than the predictions from the evolutionary models.
These latter authors suggest that there is a process during the fast case A
mass transfer that removes the angular momentum of the star; their favored
explanation is an angular momentum loss induced by a magnetic field that only
develops during the mass-transfer phase.

The disagreements on the chemical abundances and the rotation rates between
evolutionary models and the empiric analysis suggest some unconsidered physics
during the mass-transfer phase. A more fine-tuned model would likely solve
some of these problems. For example, we did not include the initial rotation
as a free parameter. This might solve the faster rotation of the primary;
however, the secondary would still be expected to be rotating too fast as it
spins up to criticality during mass transfer.

###  6.4 AzV 476 and its future evolution in the context of the most massive
star population in the SMC

The binary components of AzV 476 are among the most luminous O stars in the
entire SMC. Nevertheless, our analysis reveals that the mass of the O4 IV-
III((f))p primary star is only $M_{\mathrm{orb,\,1}}\approx 20M_{\odot}$. We
estimate its initial mass to be $M_{\mathrm{ini,\,1}}\approx 33\,M_{\odot}$
which is significantly lower than the most massive O stars in the Galaxy and
the LMC (e.g., Ramírez-Agudelo et al. 2017; Gruner et al. 2019). This system
therefore further exacerbates the problem with the deficiency of very massive
O stars in the SMC (Ramachandran et al. 2019; Schootemeijer et al. 2021).

In the future, the primary in AzV 476 will evolve bluewards on the HRD, but it
will never be able to lose its entire H-rich envelope. As its luminosity is on
the edge of the least luminous observed WR star of the SMC (Shenar et al.
2016, 2020), it is unclear as to whether or not it will be able to develop an
optically thick wind and become a WR star. Our evolutionary model predicts
that the primary will have a final mass of
$M_{\mathrm{final,\,1}}=16\,M_{\odot}$ at core helium exhaustion. The primary
will probably collapse and form a BH (Heger et al. 2003). To our knowledge, no
BH has yet been identified in the SMC.

After the primary collapse and the formation of a compact object, the system
will consist of a MS star with a compact object. The MESA model predicts that
the secondary will have $T_{\mathrm{eff}}=$28\text{\,}\mathrm{k}\mathrm{K}$$
and $\log(L/L_{\odot})=4.85$, at the halfway point of its remaining MS
lifetime. We calculated a spectroscopic model of the secondary in this
evolutionary stage and find that the secondary would be classified as a B0 V
star. Its critical rotation makes the system a likely progenitor of a high-
mass X-ray binary (HMXB) with a Be donor star. After the stellar model evolves
beyond the MS, it undergoes a mass-transfer phase, offering an additional
opportunity to show up as a HMXB. After large fractions of the secondary’s
envelope are stripped off, it will remain as a helium star with some hydrogen
left on the surface with a luminosity around $\log(L/L_{\odot})\approx 4.9$.
Using the recent estimates from Sander & Vink (2020) for massive hydrogen-free
stars, the L/M-ratio will be significantly too low to yield a sufficiently
strong stellar wind for it to become optically thick at this metallicity (even
when accounting for additional hydrogen and larger radii). Thus, our secondary
star will most likely never appear as a WR star. With a final mass of
$M_{\mathrm{final,\,2}}=7\,M_{\odot}$ at core helium depletion, the secondary
is expected to explode as a supernova leaving a NS. With a final orbital
period of about $\approx$100\text{\,}\mathrm{d}\mathrm{,}$$ this system is not
expected to merge within a Hubble time. The impact of a supernova kick is
neglected in our binary models; nevertheless, in principle, it could change
the final orbital period of the system drastically.

##  7 Summary and conclusions

In this work we study the earliest known O-type eclipsing binary in the SMC,
AzV 476. We derived the masses of its companions using a spectroscopic
analysis and by modeling its light- and RV curves. We compare the results with
binary evolutionary models.

We conducted a quantitative analysis of the multi-phase spectra in the optical
and UV obtained by the ESO VLT and the HST and supplemented by the FUSE
spectroscopy in the far-UV. To this end we used the non-LTE stellar atmosphere
model code PoWR. The spectra of both components were disentangled, allowing us
to determine their masses, stellar, and wind parameters. Independently, we
used the PHOEBE code to derive the stellar and the orbital parameters from the
light curve in the I-band as well as RV curves. We employed the MESA code to
compute detailed binary evolutionary models that are able to reproduce the
observed stellar parameters in significant detail (see Table 6 and Fig. 11).

Our conclusions can be summarized as follows:

  1. (1)

The eclipsing binary AzV 476 harbors one of the most luminous O stars in the
SMC. It consists of an O4 IV-III((f))p type primary and an O9.5: Vn type
secondary. Single stars with similar spectral types typically have masses of
${M_{\mathrm{exp,\,1}}}=49M_{\odot}$ and ${M_{\mathrm{exp,\,2}}}=16M_{\odot}$.

  2. (2)

By analyzing the light curve and the RV curve consistently, we derive orbital
masses of ${M_{\mathrm{orb,\,1}}=20M_{\odot}}$ and
${M_{\mathrm{orb,\,2}}=18M_{\odot}}$ for the binary components. We find that
the primary has less than half of the mass that is expected from its spectral
type. The modest mass of one of the most luminous stars in the SMC highlights
the conspicuous deficiency of very massive O-type stars in this metal-poor
galaxy.

  3. (3)

The spectroscopic masses agree within the uncertainties with the more reliable
orbital masses, and they confirm a mass ratio of close to unity. The observed
spectra reveal that the surface N abundance of the primary is enhanced while C
and O abundances are reduced. These two aspects suggest a prior mass-transfer
phase.

  4. (4)

The spectroscopic analysis uncovers that the mass-loss rate of the primary
with ${\log(\dot{M}/({\,M_{\odot}\,\mbox{yr}^{-1}}))=-6.1}$ is ten times
stronger compared to recent empirical prescriptions for single O stars in the
SMC. Likely, this is due to its high $L/M$ ratio. Our empirical mass-loss rate
is a factor of two higher compared to the one used in the evolutionary models.
While for our primary the result as such is in line with standard theoretical
predictions, this underlines that the present treatment of wind mass loss in
stellar evolution models needs to be improved to properly account for the
products of binary evolution.

  5. (5)

The current moderate mass of the primary can only be explained by binary
evolutionary models. According to our favorite binary model, the initial mass
of the primary was ${M_{\mathrm{ini,\,1}}=33\,M_{\odot}}$, while the secondary
formed with ${M_{\mathrm{ini,\,2}}=18\,M_{\odot}}$. The system is
$\sim$6\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$$ old.

  6. (6)

The binary evolutionary model confirms that this binary has undergone case B
mass transfer and is now in a detached phase. The binary evolutionary model
reveals that the primary must be core helium burning and that the observed CNO
abundances correspond to those of the hydrogen shell burning layers.

  7. (7)

According to our binary evolution model, the primary will become a helium star
(or maybe a WR star) with a portion of the hydrogen remaining in the envelope
($X_{\mathrm{H}}<50\%$) and finally collapse to a BH or NS. If the binary
system stays bound after core-collapse, it might show up as a HMXB with a
rapidly rotating Be-type donor. After core hydrogen depletion, the secondary
is expected to expand and transfer mass onto the compact object, stripping off
most of its hydrogen-rich envelope. The secondary will spend most of its late
evolutionary phases in the blue part of the HRD as a helium star with a small
fraction of hydrogen left in the envelope. However, due to its low $L/M$
ratio, it will never be able to display a WR-like spectrum. Finally, the
secondary will collapse, and there is a small chance it will form a binary of
compact objects that could potentially merge within a Hubble time.

The three different methods that we used to derive the stellar masses of the
binary components show that the different mass estimates agree within their
respective uncertainties (e.g., Table 6 and Fig. 11). Nonetheless, the
different methods yield somewhat different results. We point out several
difficulties that come along with the different methods applied, for instance
when measuring surface gravity (see Fig. 5). Finally, we conclude that the
earliest type eclipsing binary in the SMC, AzV 476, provides a unique
laboratory for studying massive binaries at low metallicity.

###### Acknowledgements.

We dedicate this publication to the late Dr. Rodolfo H. Barba, whose passing
is a tragic loss to the astrophysics community. As the original referee of
this paper, he provided an outstandingly detailed and substantial report,
which was extremely helpful and led to significant improvements of this work.
Dr. Barba also provided us with the extracted TESS light curve of AzV 476
which resulted in the update in the orbital period. The authors thank Dr. Alex
Fullerton for useful discussion and advise related to the FUSE spectroscopy
and Soetkin Janssens for very helpful advises regarding the PHOEBE code. We
gratefully thank the staff of the FUSE, HST, ESO telescopes for making the
useful spectroscopic data publicly available. The results presented in this
paper are based on observations obtained with the NASA/ESA Hubble Space
Telescope, retrieved from the Mikulski Archive for Space Telescopes (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. Furthermore, its conclusions are based on observations
collected at the European Southern Observatory (ESO) under the program
098.A-0049. The authors thank Andrea Mehner for preparing the OBs of the
XSHOOTU project and the people on the management committee of XSHOOTU. 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. DP acknowledges financial
support by the Deutsches Zentrum für Luft und Raumfahrt (DLR) grant FKZ 50 OR
2005. AACS is Öpik Research Fellow at Armagh Observatory & Planetarium. AACS
and VR acknowledge support by the Deutsche Forschungsgemeinschaft (DFG -
German Research Foundation) in the form of an Emmy Noether Research Group
(grant number SA4064/1-1, PI Sander). JMA acknowledges support from the
Spanish Government Ministerio de Ciencia through grant PGC2018-095 049-B-C22.

## References

  * Abdul-Masih et al. (2021) Abdul-Masih, M., Sana, H., Hawcroft, C., et al. 2021, arXiv e-prints, arXiv:2104.07621
  * Asplund et al. (2005) Asplund, M., Grevesse, N., & Sauval, A. J. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 336, Cosmic Abundances as Records of Stellar Evolution and Nucleosynthesis, ed. T. G. Barnes III & F. N. Bash (San Francisco: Astronomical Society of the Pacific), 25
  * Böhm-Vitense (1958) Böhm-Vitense, E. 1958, ZAp, 46, 108
  * Bonanos et al. (2010) Bonanos, A. Z., Lennon, D. J., Köhlinger, F., et al. 2010, AJ, 140, 416
  * Bouchet et al. (1985) Bouchet, P., Lequeux, J., Maurice, E., Prevot, L., & Prevot-Burnichon, M. L. 1985, A&A, 149, 330
  * Bouret et al. (2003) Bouret, J.-C., Lanz, T., Hillier, D. J., et al. 2003, ApJ, 595, 1182
  * Branton et al. (2021) Branton, D., Riley, A., et al. 2021, STIS Instrument Handbook, Version 20.0, Baltimore: STScI
  * Brott et al. (2011) Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115
  * Cassinelli & Olson (1979) Cassinelli, J. P. & Olson, G. L. 1979, ApJ, 229, 304
  * Castor et al. (1975) Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157
  * Cignoni et al. (2009) Cignoni, M., Sabbi, E., Nota, A., et al. 2009, AJ, 137, 3668
  * Conroy et al. (2020) Conroy, K. E., Kochoska, A., Hey, D., et al. 2020, ApJS, 250, 34
  * Cutri et al. (2003) Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog, II/246
  * De Propris et al. (2010) De Propris, R., Rich, R. M., Mallery, R. C., & Howard, C. D. 2010, ApJ, 714, L249
  * Dekker et al. (2000) Dekker, H., D’Odorico, S., Kaufer, A., Delabre, B., & Kotzlowski, H. 2000, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4008, Optical and IR Telescope Instrumentation and Detectors, ed. M. Iye & A. F. Moorwood, 534–545
  * Diaz-Cordoves & Gimenez (1992) Diaz-Cordoves, J. & Gimenez, A. 1992, A&A, 259, 227
  * Dionne & Robert (2006) Dionne, D. & Robert, C. 2006, ApJ, 641, 252
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * Gilkis et al. (2021) Gilkis, A., Shenar, T., Ramachandran, V., et al. 2021, MNRAS, 503, 1884
  * Gomel et al. (2021) Gomel, R., Faigler, S., & Mazeh, T. 2021, MNRAS, 501, 2822
  * Götberg et al. (2020) Götberg, Y., Korol, V., Lamberts, A., et al. 2020, ApJ, 904, 56
  * Gräfener et al. (2002) Gräfener, G., Koesterke, L., & Hamann, W.-R. 2002, A&A, 387, 244
  * Gruner et al. (2019) Gruner, D., Hainich, R., Sander, A. A. C., et al. 2019, A&A, 621, A63
  * Hainich et al. (2015) Hainich, R., Pasemann, D., Todt, H., et al. 2015, A&A, 581, A21
  * Hainich et al. (2014) Hainich, R., Rühling, U., Todt, H., et al. 2014, A&A, 565, A27
  * Hainich et al. (2017) Hainich, R., Shenar, T., Sander, A., Hamann, W. R., & Todt, H. 2017, in The Lives and Death-Throes of Massive Stars, ed. J. J. Eldridge, J. C. Bray, L. A. S. McClelland, & L. Xiao, Vol. 329, 171–175
  * Hamann & Gräfener (2003) Hamann, W.-R. & Gräfener, G. 2003, A&A, 410, 993
  * Hamann & Gräfener (2004) Hamann, W.-R. & Gräfener, G. 2004, A&A, 427, 697
  * Hamann et al. (2008) Hamann, W.-R., Gräfener, G., Oskinova, L., & Liermann, A. 2008, in Astronomical Society of the Pacific Conference Series, Vol. 388, Mass Loss from Stars and the Evolution of Stellar Clusters, ed. A. de Koter, L. J. Smith, & L. B. F. M. Waters, 171
  * Hamann & Koesterke (1998a) Hamann, W.-R. & Koesterke, L. 1998a, A&A, 335, 1003
  * Hamann & Koesterke (1998b) Hamann, W.-R. & Koesterke, L. 1998b, A&A, 333, 251
  * Hammer et al. (2015) Hammer, F., Yang, Y. B., Flores, H., Puech, M., & Fouquet, S. 2015, ApJ, 813, 110
  * Heger et al. (2003) Heger, A., Fryer, C. L., Woosley, S. E., Langer, N., & Hartmann, D. H. 2003, ApJ, 591, 288
  * Heger et al. (2000) Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368
  * Herrero et al. (1992) Herrero, A., Kudritzki, R. P., Vilchez, J. M., et al. 1992, A&A, 261, 209
  * Hillier et al. (2019) Hillier, D. J., Koenigsberger, G., Nazé, Y., et al. 2019, MNRAS, 486, 725
  * Hillier & Miller (1999) Hillier, D. J. & Miller, D. L. 1999, ApJ, 519, 354
  * Hirschauer et al. (2021) Hirschauer, A. S. et al. 2021, Cosmic Origins Spectrograph Instrument Handbook, Version 13.0, Baltimore: STScI
  * Holgado et al. (2020) Holgado, G., Simón-Díaz, S., Haemmerlé, L., et al. 2020, A&A, 638, A157
  * Horvat et al. (2018) Horvat, M., Conroy, K. E., Pablo, H., et al. 2018, ApJS, 237, 26
  * Hunter et al. (2007) Hunter, I., Dufton, P. L., Smartt, S. J., et al. 2007, A&A, 466, 277
  * Jones et al. (2020) Jones, D., Conroy, K. E., Horvat, M., et al. 2020, ApJS, 247, 63
  * Kippenhahn et al. (1980) Kippenhahn, R., Ruschenplatt, G., & Thomas, H. C. 1980, A&A, 91, 175
  * Klencki et al. (2020) Klencki, J., Nelemans, G., Istrate, A. G., & Pols, O. 2020, A&A, 638, A55
  * Koenigsberger et al. (2014) Koenigsberger, G., Morrell, N., Hillier, D. J., et al. 2014, AJ, 148, 62
  * Kolb & Ritter (1990) Kolb, U. & Ritter, H. 1990, A&A, 236, 385
  * Langer (1987) Langer, N. 1987, A&A, 171, L1
  * Langer et al. (1983) Langer, N., Fricke, K. J., & Sugimoto, D. 1983, A&A, 126, 207
  * Mahy et al. (2020) Mahy, L., Almeida, L. A., Sana, H., et al. 2020, A&A, 634, A119
  * Maíz Apellániz (2019) Maíz Apellániz, J. 2019, MGB: Interactive spectral classification code
  * Maíz Apellániz et al. (2012) Maíz Apellániz, J., Pellerin, A., Barbá, R. H., et al. 2012, in Astronomical Society of the Pacific Conference Series, Vol. 465, Proceedings of a Scientific Meeting in Honor of Anthony F. J. Moffat, ed. L. Drissen, C. Robert, N. St-Louis, & A. F. J. Moffat, 484
  * Marchant (2016) Marchant, P. 2016, PhD thesis, Rheinischen Friedrich-Wilhelms-Universität Bonn
  * Markova & Puls (2015) Markova, N. & Puls, J. 2015, in New Windows on Massive Stars, ed. G. Meynet, C. Georgy, J. Groh, & P. Stee, Vol. 307, 117–118
  * Markova et al. (2018) Markova, N., Puls, J., & Langer, N. 2018, A&A, 613, A12
  * Martins & Palacios (2021) Martins, F. & Palacios, A. 2021, A&A, 645, A67
  * Massa et al. (2017) Massa, D., Fullerton, A. W., & Prinja, R. K. 2017, Monthly Notices of the Royal Astronomical Society, 470, 3765
  * Massey et al. (2005) Massey, P., Puls, J., Pauldrach, A. W. A., et al. 2005, ApJ, 627, 477
  * Mazeh (2008) Mazeh, T. 2008, in EAS Publications Series, Vol. 29, EAS Publications Series, ed. M. J. Goupil & J. P. Zahn, 1–65
  * Morrell et al. (2003) Morrell, N., Ostrov, P., Massey, P., & Gamen, R. 2003, MNRAS, 341, 583
  * Muraoka et al. (2017) Muraoka, K., Homma, A., Onishi, T., et al. 2017, ApJ, 844, 98
  * Nidever et al. (2013) Nidever, D. L., Monachesi, A., Bell, E. F., et al. 2013, ApJ, 779, 145
  * Nieuwenhuijzen & de Jager (1990) Nieuwenhuijzen, H. & de Jager, C. 1990, A&A, 231, 134
  * Oegerle et al. (2000) Oegerle, W., Murphy, E., & Kriss, J. 2000, The FUSE Data Handbook, Version 1.1
  * Oskinova et al. (2011) Oskinova, L. M., Todt, H., Ignace, R., et al. 2011, MNRAS, 416, 1456
  * Pawlak et al. (2016) Pawlak, M., Soszyński, I., Udalski, A., et al. 2016, Acta Astron., 66, 421
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4
  * Penny & Gies (2009) Penny, L. R. & Gies, D. R. 2009, ApJ, 700, 844
  * Prša (2011) Prša, A. 2011, PHOEBE Scientific Reference, Version 0.30
  * Prša et al. (2016) Prša, A., Conroy, K. E., Horvat, M., et al. 2016, ApJS, 227, 29
  * Prša & Zwitter (2005) Prša, A. & Zwitter, T. 2005, ApJ, 628, 426
  * Ramachandran et al. (2018) Ramachandran, V., Hamann, W. R., Hainich, R., et al. 2018, A&A, 615, A40
  * Ramachandran et al. (2019) Ramachandran, V., Hamann, W. R., Oskinova, L. M., et al. 2019, A&A, 625, A104
  * Ramírez-Agudelo et al. (2017) Ramírez-Agudelo, O. H., Sana, H., de Koter, A., et al. 2017, A&A, 600, A81
  * Reindl et al. (2014) Reindl, N., Rauch, T., Parthasarathy, M., et al. 2014, A&A, 565, A40
  * Rivero González et al. (2012) Rivero González, J. G., Puls, J., Massey, P., & Najarro, F. 2012, A&A, 543, A95
  * Sander et al. (2015) Sander, A., Shenar, T., Hainich, R., et al. 2015, A&A, 577, A13
  * Sander & Vink (2020) Sander, A. A. C. & Vink, J. S. 2020, MNRAS, 499, 873
  * Sanyal et al. (2015) Sanyal, D., Grassitelli, L., Langer, N., & Bestenlehner, J. M. 2015, A&A, 580, A20
  * Schneider et al. (2014) Schneider, F. R. N., Langer, N., de Koter, A., et al. 2014, A&A, 570, A66
  * Schootemeijer et al. (2021) Schootemeijer, A., Langer, N., Lennon, D., et al. 2021, A&A, 646, A106
  * Seaton (1979) Seaton, M. J. 1979, MNRAS, 187, 73P
  * Shara et al. (2020) Shara, M. M., Crawford, S. M., Vanbeveren, D., et al. 2020, Monthly Notices of the Royal Astronomical Society, 492, 4430
  * Shenar et al. (2020) Shenar, T., Gilkis, A., Vink, J. S., Sana, H., & Sand er, A. A. C. 2020, A&A, 634, A79
  * Shenar et al. (2016) Shenar, T., Hainich, R., Todt, H., et al. 2016, A&A, 591, A22
  * Shenar et al. (2015) Shenar, T., Oskinova, L., Hamann, W. R., et al. 2015, ApJ, 809, 135
  * Shenar et al. (2019) Shenar, T., Sablowski, D. P., Hainich, R., et al. 2019, A&A, 627, A151
  * Simón-Díaz & Herrero (2014) Simón-Díaz, S. & Herrero, A. 2014, A&A, 562, A135
  * Soker (2000) Soker, N. 2000, A&A, 357, 557
  * Sota et al. (2014) Sota, A., Maíz Apellániz, J., Morrell, N. I., et al. 2014, ApJS, 211, 10
  * Sota et al. (2011) Sota, A., Maíz Apellániz, J., Walborn, N. R., et al. 2011, ApJS, 193, 24
  * Spruit (2002) Spruit, H. C. 2002, A&A, 381, 923
  * Tatton et al. (2021) Tatton, B. L., van Loon, J. T., Cioni, M. R. L., et al. 2021, MNRAS, 504, 2983
  * Todt et al. (2015) Todt, H., Sander, A., Hainich, R., et al. 2015, A&A, 579, A75
  * Trundle et al. (2007) Trundle, C., Dufton, P. L., Hunter, I., et al. 2007, A&A, 471, 625
  * Vanbeveren et al. (2018) Vanbeveren, D., Mennekens, N., Shara, M. M., & Moffat, A. F. J. 2018, A&A, 615, A65
  * Vernet et al. (2011) Vernet, J., Dekker, H., D’Odorico, S., et al. 2011, A&A, 536, A105
  * Vink et al. (2001) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574
  * Vink & Sander (2021) Vink, J. S. & Sander, A. A. C. 2021, Monthly Notices of the Royal Astronomical Society, 504, 2051
  * Vos et al. (2015) Vos, J., Østensen, R. H., Marchant, P., & Van Winckel, H. 2015, A&A, 579, A49
  * Walborn et al. (2002) Walborn, N. R., Howarth, I. D., Lennon, D. J., et al. 2002, AJ, 123, 2754
  * Walborn et al. (2000) Walborn, N. R., Lennon, D. J., Heap, S. R., et al. 2000, PASP, 112, 1243
  * Walborn et al. (2004) Walborn, N. R., Morrell, N. I., Howarth, I. D., et al. 2004, ApJ, 608, 1028
  * Walborn et al. (2014) Walborn, N. R., Sana, H., Simón-Díaz, S., et al. 2014, A&A, 564, A40
  * Weidner & Vink (2010) Weidner, C. & Vink, J. S. 2010, A&A, 524, A98
  * Wilson (1990) Wilson, R. E. 1990, ApJ, 356, 613
  * Wright & Eastman (2014) Wright, J. T. & Eastman, J. D. 2014, PASP, 126, 838
  * Zacharias et al. (2013) Zacharias, N., Finch, C. T., Girard, T. M., et al. 2013, AJ, 145, 44

##  Appendix A RV determination

To determine the RV shift of a star during a specific phase we use a MCMC
approach. In our code we use the synthetic spectra obtained from the spectral
analysis and shift them until the shifted synthetic spectrum best matches the
observed spectrum. In our procedure, we use the method described in the
`emcee` documentation555https://emcee.readthedocs.io with small adjustments
such that it fulfills our needs.

The probability function, a measure of how well the synthetic spectrum fits
the observed one, consists of the combination of a prior function (i.e., a
function that limits the parameter space to reasonable numbers) and a
likelihood function. Uniform distributions are used as prior functions to
avoid biases towards specific RV shifts,

$\log p_{\mathrm{prior}}(\varv_{1},\varv_{2},f)=\cases{&}\text{ if
}$-50\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$<\varv_{1}<$600\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$\\\
\,\,\,\,0&\text{ and
}$-600\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$<\varv_{2}<$50\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$\\\
\text{ and }-\infty<\log f<1.0\\\
\rule{0.0pt}{19.91684pt}-\infty\text{otherwise.}{}$ (5)

For the likelihood function we use the least square function and assume that
the variance is underestimated by some fractional amount $f$,

$\log
p_{\mathrm{likelihood}}=-\dfrac{1}{2}\sum_{n}\dfrac{(y_{n}-y_{\mathrm{model,\,}n})^{2}}{s_{n}^{2}}+\log
2\pi s_{n}^{2},$ (6)

where

$s_{n}^{2}=\sigma_{n}^{2}+f(y_{\mathrm{model,\,}n})^{2}.$ (7)

Using the above-mentioned functions, the probability function can be expressed
as

$\log P=\log p_{\mathrm{prior}}+\log p_{\mathrm{likelihood}}.$ (8)

As recommended in the `emcee` documentation${}^{\ref{emcee_footnote}}$ the
starting values are initialized as tiny Gaussian balls around
${\varv_{1}=$300\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$}$,
${\varv_{2}=$-300\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$}$ and ${\log
f=0.0}$. To ensure a good converged posterior distribution, 32 walkers (i.e.,
individual Makrov chains) are used and each one is iterated for 5000 steps. We
cut off the first 1000 steps to ensure that the chains are “burned-in” such
that the remaining distribution resembles the posterior distribution. The
uncertainties are based on the 16th, 50th, and 84th percentiles of the
samples, corresponding to approximately $1\sigma$ deviation. To make these
values readable for the PHOEBE code, which requires symmetric error margins,
only the larger uncertainty is used.

For each of the components, we take the mean value of the RVs obtained from
several absorption lines. The primary is causing most of the He ii lines seen
in the composite spectrum. Only the He ii $\lambda 4025$ line shows
contributions from both components. The secondary contributes the majority of
the observed He i absorption lines, but in many of them also the primary can
be seen. We use the He ii $\lambda 4025$ line to fit both components and to
obtain their RVs simultaneously, such that this line can be used as a cross-
check to see if these fits are in agreement with the other lines that are
purely from the primary, like He ii $\lambda 4200$, or the lines that are
associated with the secondary. This gives us confidence that the RVs of the
secondary obtained with this method are trustworthy. The obtained RVs of the
primary and secondary can be found in Table 7 and 8.

This fitting procedure turns out to work well even for line complexes with
contributions from both stars. Figure 13 shows one of our fits of the He i
$\lambda 3819$ region. The He i $\lambda 3819$ absorption line is present in
both spectra. The absorption line of the primary is shifted redwards and the
rotationally broadened absorption line of the secondary is shifted bluewards.
This complex is blended by the blueshifted He ii $\lambda 3813$ absorption
line of the primary.

Figure 13: Radial-velocity fit of the synthetic spectrum to the observed He ii
$\lambda 3813$ and He i $\lambda 3819$ absorption lines of the UVES spectrum
with ID 10 (Table 1). The observed spectrum is shown as a solid blue line, our
best fitting model as a dotted red line, the unweighted synthetic spectrum of
the primary as a dashed black line, and the unweighted spectrum of the
secondary as a dotted green line. The observed spectrum is convolved with a
Gaussian with a $\mathrm{FWHM}=$0.2\text{\,}\AA$$ and is corrected for the
velocity of the SMC and the barycentric motion. The line identification marks
correspond to the wavelengths of the absorption lines in the rest frame.

This method requires sufficiently high $\mathrm{S/N}$ observations. The MCMC
code struggles with some of the He i lines that are associated with the
secondary. In most cases, the problems emerge when the depth of the absorption
lines is of the same order as the noise (i.e., the UVES spectra have a
$\mathrm{S/N}\sim 25$). Moreover, some of the UVES spectra show a wavy pattern
from the echelle orders. In some cases, these waves overlap with absorption
lines, making an accurate RV determination almost impossible. All fits are
inspected by eye and the observations that show the aforementioned issues are
discarded from our analysis.

Table 7:  RVs of the primary determined from different spectral lines.

spectral ID | He i $\lambda 3819\,^{(a)}$ | He i $\lambda 3935$ | He ii $\lambda 4025\,^{(b)}$ | N iv $\lambda 4058$ | He ii $\lambda 4200$ | mean value(c)
---|---|---|---|---|---|---
| [$\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$] | [$\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$] | [$\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$] | [$\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$] | [$\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$] | [$\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$]
5 | $114\pm 63$ | $117\pm 22$ | $120\pm 8$ | $103\pm 7$ | $105\pm 4$ | $112\pm 14$
6 | — | $229\pm 26$ | $236\pm 8$ | $217\pm 21$ | $210\pm 4$ | $223\pm 9$
7 | $48\pm 63$ | $36\pm 47$ | $53\pm 8$ | $46\pm 9$ | $55\pm 8$ | $47\pm 16$
8 | $140\pm 13$ | $115\pm 34$ | $134\pm 11$ | $117\pm 18$ | $122\pm 4$ | $126\pm 8$
9 | $145\pm 59$ | — | $153\pm 10$ | $134\pm 22$ | $131\pm 5$ | $141\pm 16$
10 | $201\pm 53$ | $215\pm 25$ | $209\pm 11$ | $195\pm 18$ | $183\pm 8$ | $201\pm 13$
11 | $223\pm 32$ | $212\pm 26$ | $214\pm 9$ | $192\pm 13$ | $196\pm 4$ | $207\pm 9$
12 | $162\pm 18$ | $167\pm 32$ | $131\pm 14$ | $137\pm 45$ | $144\pm 8$ | $148\pm 12$
|  |  |  |  |  |

(a) The He i $\lambda 3819$ is present in the primary and secondary spectrum.
We note that the observed He i $\lambda 3819$ line also contains the He ii
$\lambda 3813$ line from the primary. (b) The He ii $\lambda 4025$ line has a
strong contribution from the secondary. (c) The mean value is calculated as
$\sum_{i=1}^{i=n}\mathrm{RV}_{i}/n$ where $\mathrm{RV}_{i}$ is the RV of line
$i$ and $n$ is the total number of lines used. The errors of the mean values
are calculated via Gaussian error propagation.

Table 8:  RVs of the secondary determined from different spectral lines.

spectral ID | He i $\lambda 3819\,^{(a)}$ | He ii $\lambda 4025\,^{(b)}$ | He i $\lambda 4143$ | He i $\lambda 4387$ | He i $\lambda 4471$ | mean value(c)
---|---|---|---|---|---|---
| [$\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$] | [$\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$] | [$\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$] | [$\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$] | [$\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$] | [$\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$]
5 | $-88\pm 29$ | $-122\pm 31$ | — | $-117\pm 20$ | $-113\pm 18$ | $-110\pm 13$
6 | — | $-164\pm 22$ | — | $-198\pm 52$ | $-152\pm 24$ | $-171\pm 21$
7 | $-41\pm 83$ | $-33\pm 23$ | — | $-32\pm 35$ | $-35\pm 31$ | $-35\pm 25$
8 | $-137\pm 89$ | $-79\pm 41$ | $-110\pm 55$ | — | $-91\pm 24$ | $-104\pm 29$
9 | $-108\pm 112$ | — | $-95\pm 88$ | $-81\pm 37$ | $-97\pm 24$ | $-95\pm 37$
10 | $-169\pm 52$ | $-199\pm 43$ | — | — | — | $-184\pm 34$
11 | $-194\pm 34$ | $-151\pm 26$ | $-213\pm 121$ | — | — | $-186\pm 43$
12 | $-149\pm 42$ | $-125\pm 79$ | $-146\pm 68$ | $-130\pm 66$ | $132\pm 42$ | $-132\pm 27$
|  |  |  |  |  |

Same footnotes as in Table 7.

##  Appendix B Additional figures

###  B.1 Additional optical N iii and N v lines

Figure 14 shows N iii and N v lines in our spectrum with the highest
$\mathrm{S/N}$, the X-SHOOTER spectrum (ID 5 in Table 1, ${\mathrm{S/N}\sim
100}$). These lines are very weak and therefore barely visible in the
spectrum. N v $\lambda 4604$ is marginally in absorption, while N iii
$\lambda\lambda 4634,4641$ appears in slight emission. Although nearly hidden
in the noise, these weak features are reproduced by our final model (red-
dotted).

###  B.2 Optical O iv and N iv lines

Walborn et al. (2004) suggest to use the optical O iv multiplets around
$3400\text{\,}\AA$ and the N iv absorption lines at wavelengths
$3463\text{\,}\AA$, $3478\text{\,}\AA$, $3483\text{\,}\AA$ and
$3485\text{\,}\AA$ to derive the N and O abundances in hot massive stars.
Luckily, the X-SHOOTER spectrum extends to such short wavelengths. Figure 15
shows that the N iv lines of our final model are in agreement with the
observation. However, for the O iv multiplet at wavelengths $3403\text{\,}\AA$
and $3414\text{\,}\AA$ the synthetic spectrum is in emission and not in
absorption as observed. For the O iv $\lambda$3414\text{\,}\AA$$ line, there
is an ISM contribution from Co I that explains the observed absorption.
Regarding the O iv $3403\text{\,}\AA$ line, we are not aware of any blending
lines. However, this particular line is formed close to the onset of the wind
and is sensitive to different details within the wind (e.g., temperature
stratification, wind velocity law). As the other O iv multiplet
$\lambda\lambda$33813389\text{\,}\AA$$ can be well reproduced, we trust our
determined O abundance, which is supported by other absorption lines in the
UV.

###  B.3 Additional wind lines

As main diagnostic wind lines, the C iv resonance doublet in the UV and the
$\mathrm{H\alpha}$ line in the optical are used. However, the He ii $\lambda
1640$ in the UV and the He ii $\lambda 4868$ in the optical are also sensitive
to the stellar wind and are taken into account when adjusting the stellar
parameters of our synthetic model. These lines are shown in Fig. 16.

Figure 14: Same as Fig. 8, but now zoomed on optical N iii and N v lines.
Figure 15: Observed (X-SHOOTER spectrum; ID 5 in Table 1) and synthetic
spectra of the region around $3400\text{\,}\AA$ showing the O iv and N iv
lines. The observed spectrum (blue) shows clear indications of ISM absorption
lines. The black line is the synthetic spectrum without a contribution of ISM
lines. The spectrum with the modeled ISM lines is shown as red dotted line.