11institutetext: Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, 3001
Leuven, Belgium; 11email<EMAIL_ADDRESS>
American Association of Variable Star Observers, 49 Bay State Road, Cambridge,
MA 02138, USA
Department of Physics and Astronomy, Embry-Riddle Aeronautical University,
3700 Willow Creek Road, Prescott, AZ 86301, USA
Département de physique and Centre de Recherche en Astrophysique du Québec
(CRAQ), Université de Montréal, C.P. 6128, Succ. Centre-Ville, Montréal,
Québec, H3C 3J7, Canada
Departamento de Astronomía, Universidad de La Serena, Av. Juan Cisternas 1200
Norte, La Serena, Chile
Department of Astronomy & Astrophysics, 525 Davey Laboratory, Pennsylvania
State University, University Park, PA 16802, USA
Department of Physics & Astronomy, University of Sheffield, Hounsfield Road,
Sheffield, S3 7RH, UK
School of Physical Sciences, The Open University, Milton Keynes, UK
Anton Pannekoek Institute for Astronomy, University of Amsterdam, Postbus
94249, 1090 GE Amsterdam, The Netherlands
Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, 85741
Garching, Germany
Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA
02138, USA
Argelander-Institut fur Astronomie, Universit ¨ at Bonn, Auf dem Hügel 71,
53121 Bonn, Germany
Department of Physics and Astronomy, University College London, Gower Street,
London WC1E 6BT, UK
Royal Observatory of Belgium, Avenue Circulaire 3, B-1180 Brussel, Belgium
Centro de Astrobiología, CSIC-INTA, Campus ESAC, Camino bajo del castillo,
E-28 692 Villanueva de la Canada, Madrid, Spain
Heidelberger Institut fur Theoretische Studien, Schloss-Wolfsbrunnenweg 35,
69118 Heidelberg, Germany
Astronomisches Rechen-Institut, Zentrum fur Astronomie der Universität
Heidelberg, Mönchhofstr. 12-14, 69120 Heidelberg, Germany
Armagh Observatory and Planetarium, College Hill, Armagh, BT61 9DG, Northern
Ireland

# The Tarantula Massive Binary Monitoring
V. R 144 – a wind-eclipsing binary with a total mass $\gtrsim 140$
$M_{\odot}$††thanks: Based on observations collected at the European Southern
Observatory under program IDs 085.D-0704(A), 085.D-0704(B), 086.D-0446(A),
086.D-0446(B), 087.C-0442(A), 087.D-0946, 090.D-0212(A), 090.D-0323(A),
092.D-0136(A), 292.D-5016(A).

T. Shenar 11 H. Sana 11 P. Marchant 11 B. Pablo 22 N. Richardson 33 A. F. J.
Moffat 44 T. Van Reeth 11 R. H. Barbá 55 D. M. Bowman 11 P. Broos 66 P. A.
Crowther 77 J. S. Clark We are saddened to report the death of Dr Simon Clark,
who passed away during the preparation of the manuscript.88 A. de Koter 99 S.
E. de Mink 1010 9 9 11 11 K. Dsilva 11 G. Gräfener 1212 I. D. Howarth 1313 N.
Langer 1212 L. Mahy 11 14 14 J. Maíz Apellániz 1515 A. M. T. Pollock 77 F. R.
N. Schneider 1616 17 17 L. Townsley 66 J. S. Vink
1818112233445566778899101011111212131314141515161617171818

(Received March 02, 2021; accepted April 06, 2021)

###### Abstract

Context. The evolution of the most massive stars and their upper-mass limit
remain insufficiently constrained. Very massive stars are characterized by
powerful winds and spectroscopically appear as hydrogen-rich Wolf-Rayet (WR)
stars on the main sequence. R 144 is the visually brightest WR star in the
Large Magellanic Cloud (LMC). R 144 was reported to be a binary, making it
potentially the most massive binary thus observed. However, the orbit and
properties of R 144 are yet to be established.

Aims. Our aim is to derive the physical, atmospheric, and orbital parameters
of R 144 and interpret its evolutionary status.

Methods. We perform a comprehensive spectral, photometric, orbital, and
polarimetric analysis of R 144. Radial velocities are measured via cross-
correlation. Spectral disentangling is performed using the shift-and-add
technique. We use the Potsdam Wolf-Rayet (PoWR) code for the spectral
analysis. We further present X-ray and optical light-curves of R 144, and
analyse the latter using a hybrid model combining wind eclipses and colliding
winds to constrain the orbital inclination $i$.

Results. R 144 is an eccentric ($e=0.51$) $74.2-$d binary comprising two
relatively evolved (age $\approx 2\,$Myr), H-rich WR stars (surface mass
fraction $X_{\rm H}\approx 0.4$). The hotter primary (WN5/6h, $T_{*}=50\,$kK)
and the cooler secondary (WN6/7h, $T_{*}=45\,$kK) have nearly equal masses,
with $M\sin^{3}i=48.3\pm 1.8$ $M_{\odot}$ and $45.5\pm 1.9$ $M_{\odot}$,
respectively. The combination of low rotation and H-depletion observed in the
system is well reproduced by contemporary evolution models that include
boosted mass-loss at the upper-mass end. The systemic velocity of R 144 and
its relative isolation suggest that it was ejected as a runaway from the
neighboring R 136 cluster. The optical light-curve shows a clear orbital
modulation that can be well explained as a combination of two processes:
excess emission stemming from wind-wind collisions and double wind eclipses.
Our light-curve model implies an orbital inclination of $i=60.4\pm
1.5^{\circ}$, resulting in accurately constrained dynamical masses of $M_{\rm
1,dyn}=74\pm 4\,M_{\odot}$ and $M_{\rm 2,dyn}=69\pm 4\,M_{\odot}$. Assuming
that both binary components are core H-burning, these masses are difficult to
reconcile with the derived luminosities ($\log L_{1,2}/L_{\odot}=6.44,6.39$),
which correspond to evolutionary masses of the order of $M_{\rm 1,ev}\approx
110\,M_{\odot}$ and $M_{\rm 2,ev}\approx 100\,M_{\odot}$. Taken at face value,
our results imply that both stars have high classical Eddington factors of
$\Gamma_{\rm e}=0.78\pm 0.10$. If the stars are on the main sequence, their
derived radii ($R_{*}\approx 25\,R_{\odot}$) suggest that they are only
slightly inflated, even at this high Eddington factor. Alternatively, the
stars could be core-He burning, strongly inflated from the regular size of
classical Wolf-Rayet stars ($\approx 1\,R_{\odot}$), a scenario that could
help resolve the observed mass discrepancy.

Conclusions. R 144 is one of the few very massive extragalactic binaries ever
weighed without usage of evolution models, but poses several challenges in
terms of the measured masses of its components. To advance, we strongly
advocate for future polarimetric, photometric, and spectroscopic monitoring of
R 144 and other very massive binaries.

###### Key Words.:

stars: massive – stars: Wolf-Rayet – binaries: close – binaries: spectroscopic
– Magellanic Clouds – Stars: individual: RMC 144, BAT99 118, HD 38282

##  1 Introduction

The upper-mass limit throughout cosmic times remains a fundamental uncertainty
in models of star formation and feedback (Figer, 2005). Theoretical
estimations of this parameter widely vary from $\approx 120$ $M_{\odot}$ to a
few thousands of solar masses, depending on the metallicity and modeling
assumptions (e.g., Larson & Starrfield, 1971; Oey & Clarke, 2005). The most
massive stars largely dictate the radiative and mechanical energy budget of
their host galaxies (Doran et al., 2013; Bestenlehner et al., 2014;
Ramachandran et al., 2019). They are invoked as chemical factories to explain
the presence of multiple populations in globular clusters (Gieles et al.,
2018; Bastian & Lardo, 2018; Vink, 2018), and are the presumed progenitors of
the most massive black holes (BH) observed with the LIGO-VIRGO detectors (e.g.
Abbott et al., 2020) and pair-instability supernovae (Fryer et al., 2001;
Woosley et al., 2007; Marchant et al., 2019).

Much effort has been dedicated to finding the most massive stars in the Local
Group. It was realised by de Koter et al. (1997) that stars with current
masses in the excess of $M\approx 100$ $M_{\odot}$ spectroscopically appear as
Wolf-Rayet (WR) stars: hot stars with emission-line dominated spectra that
form in their powerful stellar wind. While the WR phase is typically
associated with evolved, core He-burning massive stars (dubbed classical WR
stars), very massive stars already retain a WR-like appearance on the zero-age
main sequence (see also Maíz Apellániz et al. 2017 for stars just below this
limit). Being relatively hydrogen and nitrogen rich, they are usually
classified as WNh stars.

In terms of stellar mass, the current record-breaker is the WN5h star R 136a1
(alias BAT99~108), followed by the WN5h star R136c (alias BAT99~112), with
reported masses in the range of 200 to 300 $M_{\odot}$ (Crowther et al., 2010;
Bestenlehner et al., 2020). These objects are found in the core of the
Tarantula nebula in the sub-solar metallicity environment ($Z\approx
0.5\,Z_{\odot}$) of the Large Magellanic Cloud (LMC). Other reported examples
include the LMC WN5h star VFTS 682 (Bestenlehner et al., 2014) and several
very massive WNh and Of stars in the Galactic clusters Arches (e.g. Figer et
al., 2002; Najarro et al., 2004; Martins et al., 2008; Lohr et al., 2018) and
NGC~3603 (Crowther et al., 2010). Masses for these stars were inferred through
luminosity calibrations under the assumption that the stars are single. Given
the difficulty in disproving stellar multiplicity, the question of whether
these most massive stars are truly single remains open.

A more reliable way to measure the masses of stars is offered by binary
systems. With the measurement of both radial-velocity (RV) amplitudes $K_{1}$
and $K_{2}$, the minimum masses $M_{1}\sin^{3}i$ and $M_{2}\sin^{3}i$ can be
derived from Newtonian mechanics, noting the strong dependence on the orbital
inclination $i$. If the inclination can be derived, the masses then follow.
This is particularly achievable in eclipsing binaries such as WR 43a (alias
NGC 3603-A1), which enabled a Keplerian mass estimate of $M_{1}=116\pm 31$
$M_{\odot}$ for the primary component (Schnurr et al., 2008a). The inclination
can also be constrained by other means, such as polarimetry, as performed by
Shenar et al. (2017b) for R 145 in the LMC ( BAT99~119), or through an
interferometric orbit determination (e.g., WR~137 and WR~138 in the Galaxy,
Richardson et al. 2016).

In cases where the inclination is not known, mass calibrations of one or both
of the binary components to their spectral classes and/or luminosities can
convert the minimum masses to actual masses. While the calibration makes this
method model-dependent, performing this in a binary provides a much more solid
basis for the mass estimate compared to single-star estimates. The reason is
twofold. First, by resolving the orbit of the two components, the danger of
contamination by additional stellar components, as is the case for apparently-
single stars, decreases. Second, a calibration of one of the components
automatically results in the mass of the second component, which allows to
test the consistency of the calibration for both components simultaneously.
Examples include the Galactic binaries WR~21a (Tramper et al., 2016) and
WR~20a (Bonanos et al., 2004; Rauw et al., 2004). This method has recently
been applied by Tehrani et al. (2019) to the LMC WN5h+WN5h binary Melnick~34
(alias BAT99~116) to establish the most massive stellar components ever
weighed in a binary: $M_{1}=139^{+21}_{-18}$ $M_{\odot}$ and $M_{2}=127\pm 17$
$M_{\odot}$.

With a Smith visual magnitude of $\varv=11.15\,$mag (Breysacher et al., 1999),
R 144 (_a.k.a._ RMC~144, BAT99~118, Brey~89, HD~38282) is the visually
brightest WR star in the LMC. Early attempts to study R 144 by Moffat (1989)
lacked the sensitivity and time coverage necessary to uncover its
multiplicity. Through multi-epoch spectroscopy in the near IR, Schnurr et al.
(2008b, 2009) showed that R 144 is a binary candidate, but could not establish
a periodic behaviour. More recently, Sana et al. (2013b) identified R 144 as a
WNh+WNh double-lined spectroscopic binary (SB2) with a large RV amplitude
(hence possibly high stellar masses), but the authors lacked the necessary
time coverage to establish an orbital solution for this system. Given the
total luminosity $\log(L/L_{\odot})\approx 6.8$ and large RV variations
exhibited by the system, Sana et al. (2013b) suggested that R 144 might be the
most massive binary known in the Local Group.

This paper is the fifth paper in the framework of the Tarantula Massive Binary
Monitoring (TMBM) project, a multi-epoch observing campaign design to
constrain the orbital properties of about 100 massive binary candidates
detected by the VLT-FLAMES Tarantula survey (VFTS, Evans et al., 2011; Sana et
al., 2013a). While R 144 was not part of VFTS as it was used as a guiding star
for the secondary guiding of the FLAMES plate in the VFTS, R 144 was manually
inserted in TMBM in view of its science interest.

The first paper of the TMBM series provided orbital solutions for 82 O-type
systems in the sample (Almeida et al., 2017), the second focused on the very
massive WR binary R 145 (Shenar et al., 2017b), while the third and fourth
papers provided a spectroscopic and photometric analysis of double-lined
spectroscopic binaries in the sample (Mahy et al., 2020a, b). In this study,
we provide an orbital, spectroscopic, photometric, and polarimetric analysis
of R 144. We do so relying on multi-epoch data acquired with the FLAMES-UVES
and X-SHOOTER instruments mounted on the European Southern Observatory’s (ESO)
Very Large Telescope (VLT), described in Sect. 2. Our analysis is presented in
Sect. 3, with the goal of deriving its orbit, disentangling its composite
spectra, assessing the dynamical masses of both binary components, and
deriving their physical parameters. In Sect. 4, we discuss the implications of
our results on our understanding of stellar evolution at the upper-mass end
and of the upper-mass limit. We conclude our findings in Sect. 5.

##  2 Data and reduction

###  2.1 Spectroscopic data

Figure 1: X-SHOOTER spectrum of R 144 taken on 20 Feb 2014 (MJD = 56708.04),
showing an overview of the optical part of the spectrum.

Our main spectroscopic dataset includes three different instruments spanning
roughly 4 yr. The first set of data was obtained with X-SHOOTER through the
Dutch program for guaranteed observing time from Jan 2011 to Feb 2013,
allowing the detection of the system as an SB2 binary and of large RV
variations. The details of the data reduction and the early scientific results
were presented in Sana et al. (2013b). R 144 was then monitored over a 1.5-yr
period from Oct 2012 to Mar 2014 as the sole FLAMES-UVES target of the TMBM
campaign, providing 27 useful observational epochs. The FLAMES-UVES data were
obtained with the RED520 setup, providing a wavelength coverage from 4200 to
6800 $\AA$ at a spectral resolving power of 47000. Each observation was
composed of three consecutive 980s exposures for optimal cosmic removal. The
data were reduced using the ESO FLAMES-UVES CPL pipeline under the _esorex_
environment. The sky subtraction was performed using the median value of three
simultaneous sky spectra obtained through three individual FLAMES-UVES fibers
located at empty locations in the 20’-diameter field-of-view of the FLAMES
fiber positioner.

A near real-time preliminary analysis of the FLAMES-UVES data allowed us to
estimate the orbital period and time of periastron passage with sufficient
accuracy to trigger a DDT monitoring of the 2014.2 periastron passage with the
X-SHOOTER spectrograph. In this higher-cadence campaign, we obtained 13
X-SHOOTER spectra over a 10-d time base. The data were obtained in nodding-
mode, with two nodding cycles completed to guarantee optimal sky subtraction.
The nod-throw was 5” and a jitter of 1” around each nodding position was
included. NDITxDITs of 1x150s and 2x85sec were performed at each nodding
position for, respectively, the UVB/VIS arms and the NIR arm. Narrow slit
widths of 0.5” / 0.4” / 0.4” for the UVB / VIS / NIR spectra were adopted,
delivering spectral resolving powers of 9700 / 18400 / 11600. The data were
reduced following the same procedure as the X-SHOOTER GTO observations of R
144 (see Sana et al., 2013b) and normalized through a polynomial fit to the
continuum (see example in Fig. 1). The overall S/N after combining all nodding
exposures of a given arm is in the range of 250 to 350.

An additional spectrum obtained with The Fiber-fed Extended Range Optical
Spectrograph (FEROS) instrument ($R=48000,S/N\approx 60$, Kaufer et al. 1997)
mounted on the MPG/ESO-2.20m telescope in La Silla and reduced with the MIDAS
FEROS pipeline. The data were acquired on 2011-05-16 in the framework of the
Galactic O- and WN-star survey (OWN, Barbá et al., 2010), and were retrieved
from the Library of Libraries of Massive-Star High-Resolution Spectra
(LiLiMaRlin, Maíz Apellániz et al. 2019). Finally, we retrieved two reduced
UVES spectra from the ESO science portal. These were obtained in Nov 2003 with
the DIC1 390+564 instrument setup and an entrance slit of 0.7”, delivering a
resolving power of $\sim$55 000. The data provide continuous wavelength
coverage from 3260 to 6680 Å, but for a 40 Å gap around 4560 Å. The data are
described in Crowther & Walborn (2011).

In addition, we retrieved three spectra covering the far UV from the Mikulski
Archive for Space Telescopes (MAST). One spectrum covers the spectral range
$900-1200\,\AA$ and was obtained with the Far Ultraviolet Spectroscopic
Explorer (FUSE), presented by Willis et al. (2004). The spectrum was acquired
on 1999 Dec 16 (PI: Sembach, ID: P1031802000), and has a resolving power of
$R\approx 20000$ and a signal-to-noise ratio of S/N $\approx 100$. We further
retrieved two high-resolution spectra acquired with the International
Ultraviolet Explorer (IUE). The first was acquired on 1978 Sep 29 (PI: Savage,
ID: SWP02798), and covers the spectral range $1200-2000\,\AA$ with $R\approx
10000$ and S/N $\approx 15$. The second was acquired on 1979 Feb 14
(PI:Savage, ID: LWR03766), and covers $2000-3000\,\AA$ with $R\approx 10000$
and $S/N\approx 10$. All three spectra are flux-calibrated, and are thus
useful for fitting the spectral energy distribution (SED) of R 144, as well
for deriving the wind parameters of both components.

###  2.2 Photometric data

We obtained photometric data using the All-Sky Automated Survey for Supernovae
(ASAS-SN) automated pipeline (Shappee et al., 2014; Kochanek et al., 2017).
These data are available in the V, and more recently, $g$ filters. We chose to
use the V dataset as it offers a broader bandwidth and the signal in this
filter is more pronounced, covering a baseline of 1609 d (roughly 22 orbital
orbital cycles, see Sect. 3.5). We phase-folded the data on the orbital period
derived in Sect. 3.5 and removed this signal from the light-curve. We then
performed basic sigma clipping and removed clear alias/instrumental peaks
which were present in the Fourier transform. This was done by a process, of
phase-folding on the signal, templating this signal using a LOWESS filter
(Cleveland, 1979), and removing this signal from the light-curve. After
reductions were finished the initial orbital signal was added back to the
light-curve.

R 144 was also observed in several sectors of cycle 1 of the Transiting
Exoplanet Survey Satellite (TESS) mission as part of the 30-min cadence full
frame image (FFI) data (Ricker et al., 2015). More recently, R 144 was
included as a target for 2-min cadence data as part of a TESS guest
investigator proposal (PI: Bowman; GO3059). To date, 2-min cadence TESS data
in sectors 27–30 are publicly available from the Mikulski Archive for Space
Telescopes (MAST), which span a total of 108 d. We extracted a customised
light-curve of R 144 from the 2-min TESS postage stamp pixel data using a
1-$\sigma$ threshold binary aperture mask, as defined in the lightkurve python
package (Lightkurve Collaboration et al., 2018). The background flux estimates
and barycentric timestamp corrections provided by the SPOC pipeline (Jenkins
et al., 2016) were also taken into account, and the light curve was normalised
per sector by dividing through the median observed flux. To validate the
extracted light-curve, we repeated the data reduction using the TESS FFI data
with different aperture masks, and visually compared the results.

In addition, we retrieved photometry using
VizieR111https://vizier.u-strasbg.fr/viz-bin/VizieR: UBV photometry was
retrieved from Parker (1992), R magnitude from Monet et al. (1998), I
magnitude from Pojmanski (2002), JHK and IRAC magnitudes from a compilation by
Bonanos et al. (2009), and WISE magnitudes from Cutri & et al. (2014).

###  2.3 X-ray data

X-ray observations of R 144 in the range $\approx$ 0.2 - 9 keV were acquired
with the Chandra X-ray Observatory in the framework of the Chandra Visionary
Program T-ReX (PI: Townsley). The observations were analysed with the ACIS
Extract package (Broos et al., 2010). We refer the reader to Pollock et al.
(2018) and Townsley et al. (2011) for a detailed account of the observing
program and data reduction, respectively.

###  2.4 Polarimetric data

We collected linear polarimetric data for R 144 acquired and presented by
Schnurr et al. (2009). The data were obtained between 1988 October and 1990
May at the 2.2-m telescope of the Complejo Astronomico El Leoncito (CASLEO)
near San Juan, Argentina, and with Vatican Observatory Polarimeter (VATPOL,
Magalhaes et al., 1984).

##  3 Analysis

###  3.1 Radial-velocity measurements and orbital solution

Figure 2: Plotted are three X-SHOOTER spectra (MJD = 55452.89, 55585.65,
56708.04, see legend) corresponding to phases close to the radial velocity
extremes ($\phi=0.05,0.84$ with ephemeris derived in this study) and the
systemic velocity ($\phi=0.96$), and focusing on the N iv $\lambda 4058$ and N
v $\lambda 4945$ lines. Both binary components are seen in the N iv $\lambda
4058$ line, but only the hotter primary is seen in the N v $\lambda 4945$
line. Conversion to Doppler space is with reference to rest wavelengths.

The optical spectrum of R 144 is dominated by emission lines belonging
primarily to He and N, as is typical for WN stars (Fig. 1). A closer
inspection of the N iv $\lambda 4058$ and N v $\lambda 4945$ lines reveals the
presence of two WR stars in the system (Fig. 2). The two stars exhibit similar
spectra, although one star appears to be slightly hotter. This is evident when
examining, for example, the N v $\lambda 4945$ line, which forms predominantly
in one of the components. Another example can be seen in Fig. 3 for the N v
$\lambda\lambda 4604,4620$ doublet and the N iii $\lambda\lambda 4634,4641$
triplet. Both stars appear to contribute in both line blends, but the
indicative anti-phase behaviour of these lines implies that each of them is
dominated by a different stellar component. The dynamical spectra of the N v
$\lambda 4945$, N iii $\lambda\lambda 4634,4641$, and N v $\lambda\lambda
4604,4620$ lines point towards similar conclusions (Fig. 4). Throughout this
paper, we refer to the hotter component in R 144 as the primary star.

Figure 3: As Fig. 2, but showing three FLAMES-UVES spectra (MJD = 56645.57, 56645.57, 56697.69), and focusing on the N v $\lambda\lambda 4604,4620$ and N iii $\lambda\lambda 4634,4641$ lines. Both components contribute to the N v $\lambda\lambda 4604,4620$ lines, while the cooler secondary dominates the N iii $\lambda\lambda 4634,4641$ lines.  |  |
---|---|---
Figure 4: Dynamical spectra for the N v $\lambda 4945$, N iii $\lambda\lambda
4634,4641$, and N v $\lambda\lambda 4604,4620$ lines. Red and green curves
depict the orbital solution derived in Sect. 3.1 for the primary and secondary
components, respectively.

The RVs are measured via cross-correlation (Zucker & Mazeh, 1994). The
formalism is outlined in Shenar et al. (2019) and Dsilva et al. (2020). As a
first step, a single observation is used as a template against which the
remaining observations are cross-correlated to obtain the relative RVs. In the
second step, a new template is formed by co-adding all observations in the
same frame. This template is then used to re-measure the RVs, and the process
is repeated until convergence is attained (typically 2-3 iterations).

For the (hotter) primary, we use the clearly isolated N v $\lambda 4945$ line
(see Fig. 2). For the (cooler) secondary, we use the peak of the N iii
$\lambda\lambda 4634,4641$ triplet (see Fig. 3), as it is strongly dominated
by the secondary. The template of the primary, focused on the N v $\lambda
4945$ line, is calibrated to the rest-frame by comparing it with our model
atmosphere (Sect. 3.3), hence enabling us to convert relative RVs to absolute
values for the primary. As the positions of the line centroids depend on the
wind parameters, providing absolute RVs for WR stars is model-dependent, and
we estimate an accuracy of $\approx 20$ km ${\rm s}^{-1}$ on the absolute
value from comparison to models. The final measurements are given in Table 2.

Relying on the RV measurements, we derive an orbital solution, constraining
the time of periastron $T_{0}$, eccentricity $e$, argument of periastron
$\omega$, systemic velocity $V_{0}$, and RV amplitudes $K_{1}$ and $K_{2}$.
The orbital period is fixed to $P=74.2074$-d based on our light-curve
analysis, which is described in Sect. 3.5. The remaining orbital parameters
are fitting parameters. For the minimization, we use a self-written Python
script that relies on the minimization package
lmfit222https://lmfit.github.io/lmfit-py (Shenar et al., 2019), applying the
differential evolution algorithm (Newville et al., 2014). Since the RVs of the
secondary are relative, we allow for an additional constant offset between the
two RV sets to allow an orbital solution with a single systemic velocity. We
fit the RV curves simultaneously with our light-curve model, which is
described in Sect. 3.5. However, no notable differences are obtained when
minimizing both separately.

The derived orbital parameters are given in Table 3, and the best-fitting
orbital solution is shown in Fig. 5. We obtain a root mean square error (rms)
of 15 km ${\rm s}^{-1}$ for both components. In contrast, the formal
measurement errors have a mean of $\sigma\approx 10\,$km ${\rm s}^{-1}$. It
is, however, not uncommon for formal errors of RV measurements to be
underestimated in the case of WR stars, primarily due to spectral line
variability that stems from either intrinsic wind variability or wind-wind
collisions in the system. As the Keplerian model should provide a virtually
perfect description to time variability of the RVs, we multiply the formal
errors by the ratio of the measured rms value and the mean of the formal
errors, since these scaled errors should reflect more realistically the true
measurement errors.

Figure 5: RV measurements of the N v $\lambda 4945$ line (primary) and the N
iii $\lambda\lambda 4634,4641$ triplet (secondary), along with the RV curves
corresponding to the orbital parameters given in Table 3.

The inferred orbital parameters are given in Table 3. The two stars are found
to have almost identical RV-amplitudes and hence masses, with an indication
for the hotter companion being slightly more massive. The system is eccentric
and seems to strongly resemble the properties of other known WNh+WNh binaries
in the LMC such as R 145 (Shenar et al., 2017b) and Mk 34 (Tehrani et al.,
2019). The minimum masses of the primary and secondary are
$M\sin^{3}i=48.3\,M_{\odot}$ and $45.5\,M_{\odot}$, respectively.

###  3.2 Spectral disentangling

Spectral disentangling is a mathematical technique that separates a series of
phase-dependent spectra of a multiple system to the individual spectra of the
stellar constituents (e.g., Hadrava, 1995). The appearance of the separated
spectra depends on the orbital parameters (or individual RVs of the
components) as an input. Here, we use the shift-and-add technique to separate
the spectra (Marchenko et al., 1998; González & Levato, 2006). The technique
has been used in the past for the study of WR binaries (e.g., Demers et al.,
2002; David-Uraz et al., 2012; Shenar et al., 2017b, 2018, 2019). As all
orbital parameters were derived in Sect. 3.1, we can in principle fix them to
separate the spectra. However, by allowing $K_{1}$ and $K_{2}$ to vary, we can
use additional line blends to verify the results obtained in Sect. 3.1. We do
this by calculating the $\chi^{2}$ of the residuals obtained when shifting-
and-adding the disentangled spectra and subtracting the resulting spectrum
from the individual observations iteratively (see e.g., Shenar et al., 2020).

In Fig. 6, we show for the N iv $\lambda 4058$ line the match between two
observations close to RV extremes and the sum of the disentangled spectra
shifted according to the best-fitting $K_{1}$ and $K_{2}$ values. The analysis
relies on the X-SHOOTER data, which cover the relevant wavelength regime. The
lower panel shows the reduced $\chi^{2}$ map, where a clear minimum consistent
with the values given in Table 3 is shown. Figure 7 shows the same as Fig. 6,
but for the N v $\lambda\lambda 4604,4620$ lines and using both the X-SHOOTER
and FLAMES-UVES spectra. Again, the results agree very well (within 1$\sigma$)
with our measurements from the RV curves. Overall, our three measurements of
$K_{1}$ and $K_{2}$ are consistent within their respective 1$\sigma$ error
measurements and indicate that the hotter primary is slightly more massive
than its cooler companion.

---
Figure 6: Upper panels: comparison between the disentangled spectra obtained
for the best-fitting $K_{1}$ and $K_{2}$ values, their sum, and observations
at two phases close to RV extremes for the N iv $\lambda 4058$. Lower panel:
reduced $\chi^{2}$ as a function of $K_{1}$ and $K_{2}$ for the N iv $\lambda
4058$ line. The analysis is performed using the X-SHOOTER data. The minimum is
at $K_{1}=129\pm 6\,$km ${\rm s}^{-1}$, $K_{2}=134\pm 4$ km ${\rm s}^{-1}$.
---
Figure 7: As Fig. 6, but for N v $\lambda\lambda 4604,4620$. The analysis is
performed using the FLAMES-UVES data and X-SHOOTER data. The minimum is at
$K_{1}=130\pm 3\,$km ${\rm s}^{-1}$, $K_{2}=134\pm 4$ km ${\rm s}^{-1}$.
Figure 8: Disentangled spectra of R 144 obtained with the shift-and-add
technique for the hotter primary (WN5/6h, black line) and cooler secondary
(WN6/7h, green line).

The disentangled X-SHOOTER spectra are shown in Fig. 8. The spectra can be
disentangled up to a scaling factor which depends on the optical light ratio
of the two stars. For reasons that will be outlined in Sect. 3.3, we adopt the
optical light ratio of $f_{1}/f_{2}(V)=0.79$. The spectra shown in Fig. 8 are
scaled accordingly. The process of disentangling, especially of two emission-
line stars and in the presence of non-Keplerian variability (e.g., wind-wind
collisions), can introduce spurious features that should not be over-
interpreted. For example, the shape of the electron-scattering wings (notably
H and He ii lines) and P-Cygni lines (notably He i lines), as well as the
absolute strengths of the lines, can be quite degenerate in the solution.
Hence, we primarily rely on the disentangled spectra to classify the stars and
compare them qualitatively.

While the two spectra are similar, a few differences are apparent. Both
components show all three ionisation stages N iii, iv, v, but the hotter
primary shows N v more prominently, while the cooler secondary dominates in N
iii lines, as anticipated. The spectral features of the cooler secondary tend
to be sharper, whereas those of the primary more smeared, pointing towards
slightly larger terminal velocity $\varv_{\infty}$ for the hotter primary.

Following the quantitative classification schemes by Smith et al. (1996), the
primary falls in-between the WN5h and WN6h classes, while the secondary falls
in-between the WN6h and WN7h classes. This is also confirmed by morphological
comparisons with spectra presented by Crowther & Smith (1997) and Crowther &
Dessart (1998). We therefore classify the primary as WN5/6h and the secondary
as WN6/7h. The ’-h’ suffix stems from the ratio of Balmer to He ii lines,
which suggests that both still have a significant amount of hydrogen in their
envelopes.

###  3.3 Spectral analysis

Table 1: Inferred parameters for R 144 from the orbital/light-
curve/polarimetric (upper part, Sects. 3.1, 3.5, 3.6.) spectral (middle part,
Sect. 3.3), and evolutionary (lower part, Sect.4.1) analyses.

Parameter | primary | secondary
---|---|---
Spectral type | WN5/6h | WN6/7h
orbital analysis
$P_{\mathrm{orb}}$ [d] a | $74.2074\pm 0.0043$
$T_{0}$ [MJD] | $58268.98\pm 0.07$
$e$ | $0.506\pm 0.004$
$V_{0}$ [km ${\rm s}^{-1}$] b | 210$\pm 20$
$\omega$ | $304.6\pm 0.8$ | $124.6\pm 0.8$
$\Omega$ c | $114\pm 9$
$K$ [km ${\rm s}^{-1}$] | $129.5\pm 2.3$ | $137.4\pm 2.1$
$M\,\sin^{3}i$ [$M_{\odot}$] | $48.3\pm 1.8$ | $45.5\pm 1.9$
$a\sin i$ [$R_{\odot}$] | $164.0\pm 3.0$ | $174.3\pm 2.7$
$q(\frac{M_{2}}{M_{1}})$ | $0.94\pm 0.02$
$i$ [deg] d | 60.4$\pm$1.5
$a$ [$R_{\odot}$] | $189\pm 5$ | $200\pm 4$
$M_{\rm dyn}$ [$M_{\odot}$] | $74\pm 4$ | $69\pm 4$
spectral analysis
$f/f_{\rm tot}(V)$ e | 0.44 (fixed) | 0.56 (fixed)
$T_{*}$ [kK] | $50\pm 2$ | $45\pm 2$
$T_{2/3}$ [kK] | $45\pm 2$ | $40\pm 2$
$\log R_{\rm t}$ [R⊙] | $1.10\pm 0.05$ | $1.15\pm 0.05$
$\log\,L$ [L⊙] | $6.44\pm 0.05$ | $6.39\pm 0.05$
$M_{V,{\rm Smith}}$ [mag] | $-7.3\pm 0.2$ | $-7.6\pm 0.2$
$E_{B-V}$ [mag] | $0.20\pm 0.01$
$R_{*}$ [R⊙] | $22\pm 3$ | $26\pm 3$
$R_{2/3}$ [R⊙] | $27\pm 3$ | $31\pm 3$
$D$ |  $10$ (fixed) |  $10$ (fixed)
$\xi(R_{*})$ |  $30$ (fixed) |  $30$ (fixed)
$\log\dot{M}$ [$M_{\odot}\,{\rm yr}^{-1}$] | $-4.38\pm 0.15$ | $-4.34\pm 0.15$
$\varv_{\infty}$ [km ${\rm s}^{-1}$] | $1400\pm 100$ | $1200\pm 100$
$X_{\rm H}$ | $0.35\pm 0.05$ | $0.40\pm 0.05$
$X_{\rm C}/10^{-4}$ | $7\pm 3$ | $7\pm 3$
$X_{\rm N}/10^{-4}$ | $40_{-10}^{+20}$ | $40_{-10}^{+20}$
$X_{\rm O}/10^{-4}$ | $\lesssim 0.1$ | $\lesssim 0.1$
$\varv\sin i$ [km ${\rm s}^{-1}$] | $\lesssim 80$ | $\lesssim 100$
$\varv_{\rm eq}$ [km ${\rm s}^{-1}$] f | $\lesssim 100$ | $\lesssim 120$
evolution models g
$M_{\rm ev,ini}$ [$M_{\odot}$] | $148\pm 14$ | $132\pm 20$
$\varv_{\rm eq,ini}$ [km ${\rm s}^{-1}$] | $90^{+315}_{-62}$ | $10^{+36}_{-15}$
$M_{\rm ev,cur}$ [$M_{\odot}$] | $111\pm 12$ | $100\pm 12$
$\varv_{\rm eq,cur}$ [km ${\rm s}^{-1}$] | $\leq 20$ | $\leq 10$
age [Myr] | $1.96\pm 0.21$ | $2.10\pm 0.20$
333 a Obtained with the ASAS-SN and TESS data using the pwkit pdm algorithm
(Sect. 3.5). b The absolute calibration of the systematic velocity is model-
dependent and is obtained through calibration with model spectra. The 20 km
${\rm s}^{-1}$ uncertainty is a rough estimate. c Obtained from polarimetry
(Sect. 3.6). d Obtained from light-curve analysis (Sect. 3.5). e Fixed based
on mass-luminosity calibrations for homogeneous stars relying on the derived
mass ratio. f Assuming the orbital and rotational axes are aligned. g Inferred
with BONNSAI using tracks calculated by Gräfener (2021).

Figure 9: Comparison between observed SED (blue squares and lines, upper
panel) and normalized FUSE, IUE, and X-SHOOTER ($\phi=0.76$) spectra (lower
panel) and the synthetic composite spectrum (red dotted line). The composite
spectrum is the sum of the hotter primary (black solid line) and cooler
secondary (green dashed line). Prominent telluric and molecular bands are
marked in red.

To derive the physical parameters of both stellar components, we perform a
spectroscopic analysis of the system while relaxing the assumption of local
thermodynamic equilibrium (non-LTE). We use the Potsdam Wolf-Rayet (PoWR)
model atmosphere code (Hamann & Gräfener, 2003; Gräfener et al., 2002; Sander
et al., 2015). PoWR is a 1D code that solves the radiative transfer problem in
spherical geometry. It was originally developed for the analysis of WR stars,
but is applicable to any hot star with an expanding atmosphere (e.g., Shenar
et al., 2015; Giménez-García et al., 2016; Ramachandran et al., 2019). While
PoWR is designed for the analysis of single stars, the analysis of binaries is
possible by relying on the disentangled spectra (Sect. 3.2) and on the sum of
individual models calculated for the two binary components. Even though this
neglects the impact of wind-wind collisions (WWC), mutual irradiation, or
other non-spherical effects, these effects are globally negligible, although
they can dominate in the case of individual lines (e.g., the He i lines, see
Sect. 3.4).

We calculate tailored models for our analysis, but we strongly rely on PoWR
grids444http://www.astro.physik.uni-potsdam.de/PoWR calculated for LMC
metallicities (Hamann & Gräfener, 2004; Todt et al., 2015) for error
estimations. We include model atoms for H, He, C, N, O, P, Si, S, and the iron
group elements (dominated by Fe). The mass fractions of P, Si, S, and Fe are
fixed similarly to Hainich et al. (2014) and Shenar et al. (2019): $X_{\rm
P}=2.91\times 10^{-6}$, $X_{\rm Si}=3.21\times 10^{-4}$, $X_{\rm S}=1.55\times
10^{-4}$, $X_{\rm Fe}=7.02\times 10^{-4}$; the rest are fitting parameters.
The goodness of the fit is judged through visual inspection; producing
$\chi^{2}$-like estimates in the framework of full non-LTE calculations,
especially in a binary, is not feasible and would anyhow neglect systematic
errors originating in the uncertain structure of the winds, which are the
dominant factor of uncertainty here.

The main parameters defining a model atmosphere are the effective temperature
$T_{*}$, the luminosity $L$, the mass-loss rate $\dot{M}$, and the chemical
abundances, expressed here as mass fractions. Because of the optically-thick
winds of WR stars, their photospheres (defined at a mean optical depth of
$\tau_{\rm ross}\approx 2/3$) are typically located well above their
hydrostatic stellar surfaces. The inner boundary of the model, referred to as
the stellar radius $R_{*}$, is defined at $\tau_{\rm Ross}=20$.555$\tau=20$ is
a standard convention; any value above $\approx 10$ would be appropriate due
to the exponential increase of $\tau(r)$ The stellar radius $R_{*}$ is related
to $T_{*}$ and $L$ via the Stefan-Boltzmann equation
$L/L_{\odot}=(R/R_{\odot})^{2}\,(T/T_{\odot})^{4}$. The wind velocity
$\varv(r)$ assumes the functional form of a $\beta$-law (Castor et al., 1975),
characterised by a unitless parameter $\beta$ of the order of unity and the
terminal velocity $\varv_{\infty}$. We explored various $\beta$ values in the
range 1-10, and can conclude that $\beta\lesssim 3$, which is consistent with
previous literature for late-type WN stars (e.g., Chené et al., 2008). We
therefore adopt the standard value of $\beta=1$, lacking evidence to assume
otherwise. The Doppler widths of the opacity and emissivity profiles are
defined as $\varv_{\rm Dop}=\sqrt{\varv_{\rm th}^{2}+\xi^{2}}$, where
$\varv_{\rm th}(r)$ is the thermal motion of the chemical species, and
$\xi(r)$ is a prespecified microturbulence (Shenar et al., 2015). Here, we
assume $\xi(R_{*})=30\,$km ${\rm s}^{-1}$ and that $\xi(r)$ grows with the
wind velocity with a proportionality factor of 10%.

Winds of massive stars are not homogeneous but are rather clumped (Moffat et
al., 1988; Hillier, 1991; Lépine & Moffat, 1999; Puls et al., 2006; Oskinova
et al., 2007; Sundqvist et al., 2010). Optically thin clumping (microclumping)
is accounted for here by introducing the clumping factor $D$, which describes
the density factor of clumped material compared to the equivalent smooth wind
($D=1/f$, where $f$ is the filling factor). In ideal cases, $D$ can be derived
by simultaneously considering spectral features whose strength is proportional
to the density $\rho$ (P-Cygni lines, electron-scattering wings) and
$\rho^{2}$ (recombination lines, free-free emission). However, given the
uncertainties in the disentangling procedure (see Sect. 3.2) and the ambiguity
in associating UV features with the individual binary components, this is not
feasible here. While we explored various clumping parameters up to $D=100$, we
cannot derive $D$ unambiguously, and hence fix the clumping factor to the
typical value of $D=10$ (e.g., Hainich et al., 2014; Shenar et al., 2019). The
mass-loss rates derived can be easily scaled to other clumping factors by
preserving the product $\dot{M}\sqrt{D}$.

For a given value of $T_{*}$ and chemical abundances, the strength of emission
lines in visual wavelengths depends not only on $\dot{M}$, but also on the
size of the stellar surface ($\propto R_{*}^{2}$). The transformed radius
$R_{\rm t}$ (Schmutz et al., 1989), defined as

$R_{\text{t}}=R_{*}\left[\frac{v_{\infty}}{2500\,{\rm km}\,{\rm
s}^{-1}\,}\middle/\frac{\dot{M}\sqrt{D}}{10^{-4}\,M_{\odot}\,{\rm
yr}^{-1}}\right]^{2/3},$ (1)

provides a useful parameter that preserves the strength of emission
recombination lines in the model. Models with the same values of $T_{*}$,
$R_{\rm t}$, and chemical abundances would exhibit emission lines of
comparable equivalent width, irrespective of $D$, $\dot{M}$, $\varv_{\infty}$,
and $R_{*}$.

To estimate the optical light ratio in the visual between the two components,
we follow Tehrani et al. (2019) and rely on the mass ratio derived in Sect.
3.1. Gräfener et al. (2011) provided mass-luminosity relations for homogeneous
stars. Given the large masses of the components of R 144 and their large
convective cores, the assumption of homogeneity should be reasonable. The
luminosities further depend on the hydrogen mass fractions $X_{\rm H}$ of the
stars. Fixing $X_{\rm H}$ to the values derived in our analysis (see below),
we find that the luminosity of the hotter primary should be larger by about
0.05 dex. Alternatively, if we assume that both stars comprise entirely of
helium (i.e., the mass of the outer hydrogen layer is negligible), the
luminosity difference should be 0.03 dex, which is comparable. Such luminosity
differences are obtained when assuming an optical light contribution of $44\%$
for the hotter primary and $56\%$ for the cooler secondary, which we adopt
here.

The effective temperatures are determined from the balance of lines belonging
to N iii, iv, v. The wind parameters $R_{\rm t}$ and $v_{\infty}$ are
determined from the strengths and shapes of recombination lines and P-Cygni
lines. The mass fractions $X_{\rm N},X_{\rm H},$ and $X_{\rm C}$ are
determined from the overall strengths of H, He, N, and C lines. The
luminosities and reddening are derived by comparing the sum of the SEDs of
both components to available photometry and flux-calibrated FUSE and IUE
spectra. The reddening comprises two laws: a Galactic foreground contribution
following Cardelli et al. (1989) with a reddening of $E_{B-V,{\rm
Gal}}=0.03\,$mag and $R_{V}=3.1\,$, combined with the reddening law derived by
Howarth (1983) for the LMC extinction. Our final provided reddening value
represents the sum of both contributions: $E_{B-V}=E_{B-V,{\rm
Gal}}+E_{B-V,{\rm LMC}}=0.03+0.17=0.20\,$mag,

Figure 10: Schematic of R 144 as seen during primary inferior conjunction
($\phi\approx-0.026$). The primary (filled blue circle) and secondary (filled
teal circle) are plotted to scale with their relative separation. The thick
red line indicates the wind-wind collision (WWC) cone. Marked are the line-
forming regions of several diagnostic lines for the primary star, as
calculated from our tailored model atmosphere (Sect. 3.3). The line-forming
regions of the secondary star are comparable.
Figure 11: Comparison between an observed spectrum of the N v $\lambda 4945$
multiplet and model spectra for the primary calculated with. the parameters in
Table 3, but with various $v\sin i$ values (see legend). The calculation is
performed using a 3D integration of the formal integral (Shenar et al., 2014).

A comparison between the final models and observations is shown in Fig. 9, and
our derived parameters are given in Table 3. The overall global fit generally
reproduces the observations well, with a few exceptions. First, the strength
of the He i lines is not reproduced. Reproducing them consistently appears to
come at a cost of reproducing diagnostic lines such as N v $\lambda\lambda
4604,4620$. Since He i lines form far out beyond the WWC zone, as illustrated
in Fig. 10, the strengths of He i lines may strongly deviate from our model
spectra. We therefore do not consider He i lines in our fit. Second, assuming
the baseline LMC abundance adopted here is adequate, we are not capable of
reproducing the desaturated line profiles of the P v $\lambda\lambda
1118,1128$ resonance doublet, which also form in regions that could be highly
impacted by WWCs. This is a known problem that is also observed in the
analysis of OB-type stars (Fullerton et al., 2006), and is typically thought
to be related to clumpiness and porosity in the wind. Increasing the value of
the clumping factor $D$ could improve the discrepancy. However, having
calculated models with values up to $D=100$, we find that these lines remain
saturated in both the primary and the secondary. We were able to reproduce the
de-saturation when accounting for optically thick clumps in the formal
integration using the macroclumping formalism (Oskinova et al. 2007; Šurlan et
al. 2013, Hawcroft et al. 2021, submitted). However, this formalism introduces
three additional free parameters, and degrades the quality of the fit in other
spectral lines. We therefore avoid including macroclumping in our final fit.

Our derived values for the physical parameters of both components are in broad
agreement with similar objects (e.g., Shenar et al., 2017b; Tehrani et al.,
2019; Bestenlehner et al., 2020). The lower derived value of $X_{\rm H}$ in
the primary is a consequence of the ratio of pure helium lines to H+He lines,
which is slightly larger for the primary in the disentangled spectra (Fig. 8),
but should be taken with caution in light of potential contamination of line-
profile variability on the disentangled spectra.

Despite the lack of absorption lines, we can provide realistic estimates for
the projected rotational velocity $\varv\sin i$ using emission lines that form
close to the stellar surface, such as the N v $\lambda 4945$ and N iv $\lambda
4058$ lines (see Fig. 10). To account for rotation in an expanding atmosphere,
we utilise the 3D integration module in PoWR when calculating the formal
integral (Shenar et al., 2014). A comparison between various $\varv\sin i$
calculations for the primary component and observations is shown in Fig. 11,
focusing on the N v $\lambda 4945$ line. Evidently, $\varv\sin i\approx
60\,$km ${\rm s}^{-1}$ best reproduces the observed profile. However, there is
a significant degeneracy between the adopted microturbulence $\xi$ and
$\varv\sin i$, which prevents us from providing actual estimates for
$\varv\sin i$. Nevertheless, upper limits can be obtained. For the secondary
component, the N iv $\lambda 4058$ line was used. These upper limits are
provided in Table 3.

The uncertainties on $T_{*},R_{\rm t},\log L,E_{B-V},X_{\rm H},X_{\rm
C},X_{\rm N},X_{\rm O}$ are order-of-magnitude estimates based on the
sensitivity of the fit to these parameters, as explored by the few hundreds of
models calculated here. The errors on the dependent variables (e.g., the
stellar radii) follow from error propagation.

###  3.4 Wind-wind collision

When two dense winds in a massive binary collide, they form a wind-wind
collision (WWC) cone, with its tip at the region where the dynamical pressures
of both outflows equalize (Stevens et al., 1992). As the material flows along
the WWC cone, it cools down and emits light that can be seen as photometric
excess, from the X-ray regime (e.g., Cherepashchuk, 1976; Pollock, 1987) down
to the infrared and radio (Williams et al., 1997). Furthermore, as the plasma
recombines, WWC excess emission can be seen specifically in recombination
lines (e.g. Rauw et al., 1999; Marchenko et al., 2003; Shenar et al., 2017b).

To study whether excess WWC emission may influence the data, we show in Fig.
12 the equivalent widths (EW) of several diagnostic spectral lines as a
function of orbital phase. It is apparent that almost all lines show excess of
flux close to periastron, with the exception of the N v $\lambda\lambda
4604,4620$ doublet. In fact, the WWC emission may increase the line flux by up
to $\approx 100\%$, but it is absorbed by the stellar winds close to
periastron. This can be seen, for example, in the EW change of the He ii
$\lambda 4686$, line (Fig. 12), which increases rapidly towards periastron,
but then drops at periastron. Since the emission lines contribute 5-10% to the
visual flux, such a variation in their strengths can lead to changes in the
visual light-curve of the order of $\approx 10\,$%, as we indeed observe
(Sect. 3.5).

|
---|---
|
|
Figure 12: EW variations with phase for He i $\lambda 4472$ (4465–4490 Å), He
i $\lambda 5876$ (5875–5900 Å), He ii $\lambda 4686$ (4670–4720 Å), H$\beta$
(4840–4890 Å), N iii $\lambda\lambda 4634,4641$ (4620–4650 Å), and N v
$\lambda\lambda 4604,4620$ (4595–4615 Å). Only the red part of the He i lines
is integrated to avoid the excessive blue-shifted absorption, which we
interpret as line-of-sight absorption by the shock cone

###  3.5 Visual light-curve

Figure 13: ASAS-SN light-curve of R 144 phased with the ephemeris given in
Table 3. Marked in the inset are the phases of inferior and superior
conjunction (hotter primary in front and behind the cooler secondary,
respectively) and periastron passage.

The ASAS-SN optical light-curve of R 144 shows a clear periodic signature (see
Fig. 13), which covers roughly 22 orbital cycles. The same signature is seen
in our extracted TESS light-curve, which covers 1.5 orbital cycles (Appendix
B). To derive the period, we combined the ASAS-SN and TESS data, which
together cover a time base of 6.5 yr. We then conducted a Period search with
Phase Dispersion Minimization (pwkit pdm, Stellingwerf, 1978; Schwarzenberg-
Czerny, 1997) to derive the period of R 144 of $P=74.2074\pm 0.0043$ d (Table
3). The period is consistent with that found from fitting the RV curves, but
given the sharp minimum observed in the light-curve, we fix the period to that
found from the light-curve.

Figure 13 shows the ASAS-SN light-curve of R 144 phased with the ephemeris
given in Table 3 and binned at $\Delta\phi=0.002$. The light-curve mimics that
of heartbeat stars (Thompson et al., 2012) – a class of eccentric binaries
that show rapid periodic variations in their light-curve, typically associated
with tidal distortion of the star(s) during periastron passage. We attempted
to model the ASAS-SN light-curve using the Physics of Eclipsing Binaries
(PHOEBE) light-curve modelling tool (Prša et al., 2016), using our derived
stellar and orbital parameters as input. Our PHOEBE model is shown in Appendix
C. However, the observed amplitude of the ”heartbeat” signal is roughly 100
times larger than obtained in the model, and the shape was not well
reproduced, suggesting that a different mechanism is responsible for the
observed behaviour. As we show below, the light-curve can be explained by a
combination of phase-dependent excess emission and eclipses.

To model the light-curve, we construct a hybrid model containing two physical
ingredients. First, as discussed in Sect. 3.4, recombination photons stemming
from the WWC region can cause flux changes of the order of several percent.
This is confirmed from the variation of integrated line-flux in our spectra,
and is observed in other WR (e.g., $γ^2$ Vel – Richardson et al. 2017). For
large separations where adiabatic conditions prevail, this emission is
expected to grow roughly as $E_{\rm WWC}\propto 1/\mathcal{D}(\phi)$ reaching
a maximum at periastron ($\phi=0$), where $\mathcal{D}(\phi)$ is the
instantaneous separation between the stars (Usov, 1992). We model the excess
emission stemming in WWC as

$E_{\rm WWC}(\phi)=A_{\rm WWC}\cdot\left(\frac{\mathcal{D}_{\rm
min}}{\mathcal{D}(\phi)}\right)^{\gamma_{\rm WWC}}=A_{\rm
WWC}\cdot\left(\frac{1+e}{1+e\,\cos\nu}\right)^{\gamma_{\rm WWC}},$ (2)

where $\nu(\phi)$ is the true anomaly, $A_{\rm WWC}$ is the maximum excess at
periastron, and $\gamma_{\rm WWC}$ is a number of the order of unity. $A_{\rm
WWC}$ and $\gamma_{\rm WWC}$ are treated as free fitting parameters.

Second, we account for eclipses in the system. Given its shape and strength,
the flux minimum observed precisely at inferior conjunction ($\phi=-0.026$,
primary in front of secondary) is not an ordinary eclipse, but rather a ”wind
eclipse”. As one star is behind the other, the light from the eclipsed star is
scattered off the free electrons in the wind of the eclipsing star (Lamontagne
et al., 1996; St-Louis et al., 2005). Lamontagne et al. (1996) developed a
model for the case of O+WR binaries in circular orbits involving a single
eclipse of the O star by the WR wind. Here, we extend this model to elliptical
orbits by substituting the angle $2\pi\phi$ with the angle666Sometimes also
denoted as $\lambda$ $\varpi=\nu+\omega+\pi/2$ (defined to be 0 when the
primary eclipses the secondary, i.e., at $\nu=-\omega-\pi/2$), and
substituting the constant separation $a$ in the case of a circular orbit with
$\mathcal{D}(\varpi)$. We further consider the fact that both stars have
significant winds and are eclipsing each other by adding the corresponding
terms for the companion. We adopt the terminal velocities and stellar radii
$R_{*}$ given in Table 3. Aside from the orbital parameters
$P,e,\omega,T_{0}$, the hybrid model contains six fitting parameters: $A_{\rm
WWC}$, $\gamma_{\rm WWC}$, $\dot{M}_{1}$, $\dot{M}_{2}$, $i$, and the visual
light ratio $f_{1}/f_{2}(V)$. However, we cannot fit the light ratio
simultaneously to both $\dot{M}_{1},\dot{M}_{2}$ as the model becomes
degenerate. Therefore, we fix the light ratio to the value estimated in Sect.
3.3 ($f_{1}/f_{2}=0.79$), but discuss the impact of changing this value in
Appendix D.

The model relies on the following assumptions: (i) The stars can be
approximated as point sources. With the derived geometry of R 144, this should
be a fairly good approximation (however, see discussion in Appendix D). (ii)
The optical depth along the line-of-sight is dominated by Thomson scattering
off free electrons. This assumption is verified by our model atmospheres,
where electron scattering contributes $\approx 80\%$ to the total opacity in
the outer layers. (iii) The scattering is optically thin. This is verified by
our model, where the optical depths typically remain well below unity. (iv)
The wind velocity field takes the form of a $\beta$-law. Analytical solutions
are provided by Lamontagne et al. (1996) for the $\beta=0$ (constant velocity)
and $\beta=1$, but arbitrary values of $\beta$ may be assumed if the
integration is performed numerically. For now, we adopt $\beta=1$, but explore
$\beta=2$ in Appendices D and E. (v) The wind is fully ionized, which yields
the number of electrons per baryon $\alpha\simeq(1+X_{\rm H})/2$, resulting in
$\alpha_{1}=0.68$ and $\alpha_{2}=0.7$. Our model atmospheres suggest that the
assumption of fully-ionized winds holds reasonably well for the primary, but
breaks down at a few stellar radii above the stellar surface for the
secondary. For now, we fix $\alpha$ as constant, but discuss this assumption
in Appendix D. (vi) Both stellar winds are spherically symmetric and
unperturbed by each other, and no additional opacity sources are present. We
discuss the validity of this assumption in Appendix D.

Figure 14: Upper panel: Comparison between the observed ASAS-SN light-curve
and our best-fitting hybrid light-curve model (red solid line) consisting of
double wind-eclipse (blue dashed line) + WWC excess emission (green dotted
line), which are multiplied by 0.5 for clarity. Lower panel: Residuals between
observation and model (O-C).

The ASAS-SN light-curve is fitted simultaneously with the RV curves derived in
Sect. 3.1. We refrain from combining the TESS data with the ASAS-SN data in
the fit. This is both because the time span of the TESS light-curve only
covers 1.5 orbital periods (i.e., only one wind-eclipse is seen), and also
because the amplitude of the eclipse is sensitive to the choice of aperture
mask and background subtraction in the TESS reduction procedure due to the
relatively large pixel size of TESS. However, the TESS data allow for a
refinement of the orbital period, and provide an encouraging sanity check for
our model. A comparison between our model and the TESS data is presented in
the Appendix B. A comparison between our model and the ASAS-SN data is shown
in Fig. 14. The derived parameters are $A_{\rm WWC}=0.12\pm 0.01$,
$\gamma_{\rm WWC}=1.5\pm 0.1$, $\log\dot{M}_{1}=-4.29\pm
0.05\,$$M_{\odot}\,{\rm yr}^{-1}$, $\log\dot{M}_{2}=-4.42\pm
0.05\,$$M_{\odot}\,{\rm yr}^{-1}$, and $i=60.4^{\circ}$ (given also in Table
10).

Figure 15: Same as Fig. 14, but showing the impact of changing the inclination
by $\pm 5^{\circ}$ while fixing all other parameters (green and blue lines,
see legend).

The agreement between the ASAS-SN data and our analytical model is excellent,
with a formal reduced $\chi^{2}$ for the RV curves + light-curve combined of
$\chi^{2}=0.7$. As this indicates that the formal errors on the photometry are
overestimated, we rescale them with a constant factor (0.82) to obtain
$\chi^{2}$ values of the order of unity, which ensures that the relative
weighting between the photometric and RV measurements is adequate777The
conclusion that the noise is overestimated assumes that the model is adequate
and the noise is Gaussian. Nevertheless, we note that we get very similar
results when no correction is applied to the errors.. Figure 14 also shows the
residuals between observation and computation (O-C). No clear discrepancies
are seen. However, the O-C plot reveals a low-frequency modulation, also
observed in other WR stars (e.g. Ramiaramanantsoa et al., 2019). However, a
periodogram of the residuals did not reveal any signficant frequency peaks.

While it is not immediately evident, the light-curve is impacted by both
eclipses. The eclipse at primary superior conjunction ($\phi=0.28$, secondary
is in front) is responsible for the rapid decrease of flux after periastron
and to the slight discontinuity observed at $\phi=0.28$. It is wider than the
eclipse at $\phi=-0.026$ due to the longer time spent at apastron, and it is
weaker primarily due to the larger separation between the two components,
resulting in a larger impact parameter and hence lower column density.

The derived parameters are in very good agreement with expectation. The mass-
loss rates, $\log\dot{M}_{1}=-4.29\pm 0.05\,$[$M_{\odot}\,{\rm yr}^{-1}$] and
$\log\dot{M}_{2}=-4.42\pm 0.05\,$[$M_{\odot}\,{\rm yr}^{-1}$], agree within
$1\sigma$ with those derived from spectroscopy. The WWC excess reaches a
maximum of $\approx\pm 10$%, which is roughly consistent with the flux changes
in recombination lines (Fig. 12, Sect. 3.4). The exponent of $\gamma_{\rm
WWC}=1.5\pm 0.1$ is in the accepted range of 1-2 (e.g. Usov, 1992; Schnurr et
al., 2009; Shenar et al., 2017b).

Importantly, the modelling of the light-curve enables us to constrain the
orbital inclination. In fact, the model is extremely sensitive to the
inclination. This is illustrated in Fig. 15, where we show the impact of
changing $i$ by $\pm 5^{\circ}$ when fixing all other parameters. Our fit
implies an inclination of $i=60.4\pm 1.5^{\circ}$. This, in turns, implies
dynamical masses of $74\pm 4\,M_{\odot}$ and $69\pm 4\,M_{\odot}$ for the
primary and secondary, respectively. The validity of our model and an
investigation for possible systematic errors are thoroughly discussed in the
Appendix D.

###  3.6 Polarimetric analysis

---
Figure 16: Comparison between the observed and modelled Stokes parameters $Q$
(blue line) and $U$ (green line) for two distinct exponents $\gamma$. For
$\gamma=1$ (upper panel), we obtain $\Omega=120\pm 11^{\circ}$,
$\tau_{*}=0.22$, $Q_{0}=-0.24\pm 0.02$, $U_{0}=0.16\pm 0.02$. For $\gamma=-1$
(lower panel), we obtain $\Omega=125\pm 11^{\circ}$, $\tau_{*}=0.10\pm 0.02$,
$Q_{0}=-0.23\pm 0.02$, $U_{0}=0.15\pm 0.02$. Dashed lines bounding the curves
show the impact of changing the inclination by $\pm 10^{\circ}$. See text for
details.

An additional observational constraint that is sensitive to the inclination is
provided by polarimetry. R 144 was studied using spectropolarimetry by Vink &
Harries (2017), who reported mild depolarization in recombination lines (”the
line effect”), presumably due to the binary nature of R 144. Time-dependent
polarimetry can yield important constraints on the geometry of the system,
including the orbital inclination. When the light of a star is scattered off
the free electrons in the wind of its companion, a net polarization can be
observed in the Stokes parameters $Q$ and $U$ that varies with time. Brown et
al. (1978, 1982) developed analytical equations that describe the change in
$Q$ and $U$ in elliptical systems where one star possesses a wind888The
equations were later corrected by Simmons & Boyle (1984).. This model has
since been implemented in several cases of WR binaries (St. -Louis et al.,
1988; Moffat et al., 1998), including the WR+WR binary R 145 (Shenar et al.,
2017b). The phase change of the Stokes parameters depends not only on the
orbital parameters $P,T_{0},e,\omega$, and $i$, but also on the longitude of
the ascending node $\Omega$. Moreover, the model depends on the density and
structure of the scattering medium, which is typically parametrised with two
free fitting parameters: an effective optical depth $\tau_{*}$, and an
exponent related to the structure of the wind, $\gamma$. For brevity, we refer
to Shenar et al. (2017b) for a full account of the relevant equations.

Unfortunately, only a few polarimetric measurements are available for R 144,
and they are obtained relatively far from periastron. However, we can still
fix most of the parameters based on our previous analyses (cf. Table 3),
leaving only $\Omega$, $\tau_{*}$, $\gamma$ and two arbitrary offset constants
$Q_{0}$ and $U_{0}$ as fitting parameters. The fitting procedure, performed
using the Python lmfit package, reveals that $\gamma$ and $\tau_{*}$ cannot be
well constrained with the available data, while $\Omega$ can be. Therefore, we
provide our preliminary value of $\Omega$ in Table 3. In Fig. 16, we show the
comparison between the $Q/U$ data and our-best-fitting models for two typical
values of $\gamma=-1$ and $1$. The impact of changing the inclination is also
shown.

The few data points available, along with their relatively large errors, can
accommodate a wide range of inclination angles, and we therefore cannot
provide an independent measurement of $i$, which is critical for the
estimation of the dynamical masses. Future polarimetric monitoring of the
system should is therefore strongly encouraged.

###  3.7 X-ray light-curve

Figure 17: X-ray light-curve of R 144, also shown aggregated into adaptive
wider phase intervals. See text for details.

As discussed in Sect. 3.4, with shock velocities of the order of $\approx
1000\,$km ${\rm s}^{-1}$, WWCs are expected to lead to substantial X-ray
emission (e.g., Cherepashchuk, 1976; Stevens et al., 1992; Usov, 1992). X-ray
excess is indeed often, though not always, observed in binaries with two
massive components (e.g., Pollock, 1987; Corcoran et al., 1996; Nazé et al.,
2007; Guerrero & Chu, 2008). Similarly to the optical emission excess in
recombination lines, one would expect the X-ray emission to increase at
periastron, where the densities are highest, and decrease inversely
proportional to the separation.

Interestingly, the exact opposite is observed in R 144. In Fig. 17, we show
the T-ReX X-ray light-curve of R 144, folded with the ephemeris given in Table
3. Tehrani (2019) reported R 144 to be a fairly hard X-ray source ($kT=4.2\pm
0.4\,$keV) and reported with an X-ray luminosity of $\log L_{\rm
X}=33.4\,$${\rm erg}\,{\rm s}^{-1}$. For such a luminous binary, this is
relatively modest. In comparison, the massive WNh+WNh binary Mk 34, despite
sharing many of the characteristics of R 144 and despite portraying a
comparable X-ray hardness ($kT=3.16\pm 0.03\,$keV), exceeds this X-ray
luminosity by almost two orders of magnitudes at peak luminosity ($\log L_{\rm
X}=35.3\,$${\rm erg}\,{\rm s}^{-1}$, Pollock et al. (2018)). Moreover, while
the X-ray flux of Mk 34 increases shortly before periastron (followed by a
rapid decrease), it decreases at periastron in the case of R 144. This is in
contrast to the pattern observed in the optical light-curve.

The difference in the X-ray behaviour between Mk 34 and R 144 may be related
to their different wind properties. The components of Mk 34 have reported
terminal velocities that are roughly twice as large as reported here for R
144, but they also have mass-loss rates that are about three times lower than
reported here. Additionally, with an orbital period of $P=155\,$d, the
separation between the components in Mk 34 is significantly larger than found
for R 144. The smaller semi-major axis in R 144 in combination with the
stronger winds (compared to Mk 34) could mean that a substantial amount of
X-rays is absorbed by the stellar winds. Moreover, at the winds likely do not
reach their terminal value at the collision region. Assuming $\beta=1$ for the
velocity law, the wind velocities are roughly 75% (i.e., $\approx 1000\,$km
${\rm s}^{-1}$) at collision during periastron, and even lower for larger
$\beta$ values. Hydrodynamical simulations will be needed to assess whether
these facts are sufficient to explain the two orders-of-magnitude difference
observed in their X-ray luminosities. Importantly, R 144 shows that, while
detection of strong X-rays is a helpful indicator of multiplicity, the
converse does not hold.

##  4 Discussion

###  4.1 Evolutionary status

We estimate the evolutionary masses and ages of both components using the
Bayesian statistics tool BONNSAI999 The BONNSAI web-service is available at
www.astro.uni-bonn.de/stars/bonnsai (Schneider et al., 2014). Using a set of
input stellar parameters and their corresponding errors, the tool estimates
the current and initial masses along with the age by interpolating between
evolutionary tracks calculated at LMC metallicity by Brott et al. (2011) and
Köhler et al. (2015) for stars with initial masses up to 500 $M_{\odot}$ and
over a wide range of initial rotation velocities. As input parameters, we use
the derived values of $T_{*}$, $\log L$, $X_{\rm H}$, and upper limits on
$\varv_{\rm eq}$. However, the BONNSAI tool cannot find a solution that
satisfies these three constraints. We therefore advanced by omitting the
condition on the bounded rotational velocity.

Figure 18: HRD positions of the primary and secondary components of R 144,
along with evolution tracks calculated by Köhler et al. (2015) for $M_{\rm
ini}=100,125,$, and 150 $M_{\odot}$ with an initial rotation of $\varv_{\rm
eq,ini}=200\,$km ${\rm s}^{-1}$ and 350 km ${\rm s}^{-1}$. The color
corresponds to the surface hydrogen mass fraction.

Our final results are shown in Table 3. BONNSAI infers an age of roughly 2 Myr
for both components. The consistent ages suggest that the components of R 144
have not interacted via mass-transfer in the past, as expected from the
orbital separation and the significant eccentricity. We derive current masses
of $M_{\rm ev,cur}=110\pm 9$ and $100\pm 10$ $M_{\odot}$ and initial masses of
$M_{\rm ini}=130\pm 13$ and $119\pm 12$ $M_{\odot}$ for the primary and
secondary, respectively (Fig. 18).

####  4.1.1 Rotation discrepancy

While the ages of both components are consistent within their errors, the Bonn
models are only capable of reproducing the observed parameters when assuming
high initial rotation rates of $\varv_{\rm eq,ini}\approx 350-400\,$km ${\rm
s}^{-1}$, which tend to chemically homogenize the stars. This in turn results
in current rotation velocities of $\varv_{\rm eq,cur}\approx 250\,$km ${\rm
s}^{-1}$, which are substantially larger than the conservative upper limits
derived here ($\varv_{\rm eq}\lesssim 100$ km ${\rm s}^{-1}$). As is
illustrated in Fig. 18, models without very high initial rotation are not
capable of reproducing the strong He enrichment that we observe.

The mass-loss rates implemented by Köhler et al. (2015) rely on extrapolations
from Vink et al. (2001). However, recent empirical (Bestenlehner et al., 2014,
2020) and theoretical studies (Gräfener & Hamann, 2008; Vink et al., 2011;
Schneider et al., 2018a; Bestenlehner, 2020; Gräfener, 2021) suggest that
mass-loss rates at the upper-mass end are substantially larger than the
original Vink et al. (2001) prescriptions. Indeed, the current mass-loss rates
retrieved by the BONNSAI tool are of the order of
$\log\dot{M}=-4.8\,$[$M_{\odot}\,{\rm yr}^{-1}$], more than a factor two lower
than obtained empirically in our study. Higher mass-loss rates would result in
a more rapid spin-down of the stars, potentially bringing the rotation in
agreement with observation. Moreover, higher mass-loss rates would enable the
stars to strip themselves off their H-rich layers and hence lead to
atmospheric H-depletion without the need of invoking very high rotation to
homogenise them. Hence, boosting the mass-loss rate would have a ”double”
effect in this regard.

To explore this further, we compared our results with MESA evolution code
(Paxton et al., 2011) calculated by Gräfener (2021), which include enhanced
mass-loss rates. The tracks were integrated into the BONNSAI tool, allowing
for a derivation of the most probable initial masses, rotation velocities, and
ages. The final results are presented in Table 3, and representative tracks
are shown in Fig. 19. Evidently, the tracks reproduce the observed HRD
positions simultaneously to the observed H-depletion without the need to
invoke high rotation. Due to the larger mass-loss rates, the initial masses
are larger than those inferred by BONNSAI, but the ages and currernt masses
are comparable to those derived with BONNSAI. Since the MESA tracks better
represent the properties of the system, and given the evidence for increased
mass-loss, we adopt the evolutionary parameters from the MESA tracks and
provide them in Table 3.

Figure 19: As Fig. 18, but showing MESA tracks with parameters corresponding
to those given in Table 3. These tracks are calculated using the setup
described by Gräfener (2021), which includes boosted mass-loss rates.

There are additional mechanisms that could be relevant in the context of
rotation and H-depletion. For example, it is also possible that additional
mixing mechanisms operates at these very high masses (e.g., Higgins & Vink,
2019; Jermyn et al., 2018; Pedersen et al., 2018; Schootemeijer et al., 2019;
Bowman et al., 2019; Gilkis et al., 2021). Moreover, tidal interaction at
periastron may have also operated in R 144, slowing the components down
towards synchronicity. At periastron, the ratio between the separation and
diameter of the stars is $a_{1}\,(1-e)\,{R_{1}}^{-1}\approx
a_{2}\,(1-e)\,{R_{2}}^{-1}\approx 4$. At such ratios, tidal synchronization of
the rotational velocities is expected to become important on a nuclear
timescale (see table 2 in Zahn, 1977). While these mechanisms do not seem to
be necessary to explain the observed H-depletion and low rotational
velocities, it is possible that they have operated in R 144.

####  4.1.2 Mass discrepancy

The derived evolutionary masses ($M_{\rm ev,cur}\approx
110\,M_{\odot},100\,M_{\odot}$) disagree with our dynamical masses ($M_{\rm
dyn}=74\pm 4\,M_{\odot}$, $69\pm 4\,M_{\odot}$), which at least formally are
determined to a high degree of accuracy. Interestingly, a similar discrepancy
was found by Shenar et al. (2017a) for the binary R 145, which shares many
similarities with R 144. The combination of high stellar luminosities and
relatively low masses imply that the binary components are situated close to
the Eddington limit, with Eddington factors of $\Gamma_{\rm e}=0.78\pm 0.10$
for both components. Considering the additional pressure from line-driving, it
is hard to imagine that the stars should not be strongly inflated at such
Eddington factors (Petrovic et al., 2006; Gräfener et al., 2011, 2012;
Grassitelli et al., 2016; Sanyal et al., 2017). However, the radii and
temperatures derived here are consistent with expectation for non-inflated
main sequence stars (Gräfener et al., 2011).

This discrepancy has two possible resolutions. On the one hand, it is possible
that the errors of one or more of the orbital parameters are underestimated.
Of all orbital parameters, the orbital inclination is probably the most
critical one here. An inclination of $50^{\circ}$ – a mere $10^{\circ}$
decrease – would yield masses that are comparable to the evolutionary masses.
However, the formal error on the inclination is only $1.5^{\circ}$, as the
eclipse light-curve model is extremely sensitive to $i$ (see Fig. 15).
Therefore, if the error is underestimated, it must be systematic in nature. We
explored various assumptions in the light-curve model, and these attempts are
described in detail in Appendix D. We could not isolate a potential source of
error to justify a $10^{\circ}$ difference. Moreover, in the framework of our
light-curve model, reducing the inclination down to $50^{\circ}$ also requires
very large mass-loss rates of the order of
$\log\dot{M}=-3.9\,$[$M_{\odot}\,{\rm yr}^{-1}$] to reconcile for the decrease
of $i$, which is inconsistent with our spectral analysis.

Taken at face value, the relatively low dynamical masses may suggest an
entirely different interpretation for the evolutionary status of the system:
the luminosities and masses would be consistent with fundamental structure
models when assuming that both stars comprise primarily He in their cores,
surrounded by relatively thin layers containing residual hydrogen (e.g.
Schootemeijer & Langer, 2018). In other words, the two components of R 144 may
be He-burning, classical WR stars rather than main-sequence WR stars. It is
interesting to note that the morphology of the spectra resembles those of
classical WR stars more than those of main-sequence WR stars (see examples in
Crowther & Smith, 1997; Crowther & Dessart, 1998).

In this context, it is also interesting to consider the less evolved ”twin” of
R 144, Mk 34 (Tehrani et al., 2019). Mk 34 was reported to host two stars of
$\approx 140+130\,M_{\odot}$ (noting that these values rely on calibration to
evolution models) that have roughly the same luminosities and temperatures as
the components of R 144. However, they exhibit significantly lower mass-loss
rates (smaller by a factor of three) and terminal velocities that are roughly
twice as large compared with the components of R 144, which, empirically, are
more typical for main sequence stars. If R 144 indeed comprises two classical
WR stars, the components of R 144 may well be highly inflated, since they
appear much larger than what one would expect ($R_{*}\approx 25\,R_{\odot}$
vs. $R_{*}\approx 1\,R_{\odot}$), as also observed for a multitude of Galactic
classical WR stars (Petrovic et al., 2006). Hence, R 144 could be a rare
binary of two very massive, inflated classical WR stars, potentially the most
massive classical WR stars ever weighed. Models by Gräfener (2021) suggest
that such stars would have initial masses comparable to those obtained when
assuming that the stars are homogeneous ($\approx 150\,M_{\odot}$ and
$130\,M_{\odot}$), but their age would be roughly 3 Myr instead of the 2 Myr
derived here.

###  4.2 Future evolution

The similarity of both binary components and the fact that the system is still
eccentric suggest that the stars have not interacted via mass-transfer in the
past. On the other hand, given the evolutionary status of the system, the
binary components of R 144 are also unlikely to interact in the future. In
fact, the final orbital period $P_{\rm f}$ between the two would only grow
with time due to wind mass-loss as $P_{\rm f}/P_{\rm cur}\propto\left(M_{\rm
tot,cur}/M_{\rm tot,f}\right)^{2}$, where $M_{\rm tot,cur}$ and $M_{\rm
tot,f}$ are the total current and final mass of the system. With with an
initial mass in excess of $140\,M_{\odot}$, the primary could be massive
enough to undergo a pair-instability supernova, leaving no black hole behind
(e.g., Langer et al., 2007). Assuming that the fate of both components is to
end their lives as $\approx 20-30$ $M_{\odot}$ BHs (Belczynski et al., 2010),
we find that $P_{\rm f}\approx 1-5\,$yr (depending on the current masses).
That is, the system would become a BH+BH binary with an orbital period of a
few years.

At such large periods, interactions with tertiary companions become much more
likely (e.g. Toonen et al., 2020). With an integrated absolute V-band
magnitude of $M_{V}\approx-8\,$mag, detecting any companions with current
masses $\lesssim 50$ $M_{\odot}$ (corresponding to light contributions smaller
than $\approx 10$%) via spectroscopy would be difficult, and it is therefore
well possible that companions with masses as large as $\approx 50$ $M_{\odot}$
reside in the vicinity of R 144. A tertiary component could act to reduce the
period of the binary through the Kozai-Lidov mechanism (Lidov, 1962; Kozai,
1962), eventually forcing the BHs to merge in a gravitational-wave event
(e.g., Vigna-Gómez et al., 2021). Alternatively, interactions with a third
star could lead to instabilities that either eject the BHs at runaway
velocities, or bind one of them to the tertiary star, forming a star+BH
binary. Such objects could be the progenitors of high-mass X-ray binaries.

###  4.3 The upper mass limit

Regardless of the evolutionary status of the system, our results imply initial
masses of the order of $\approx 130-150\,M_{\odot}$ for the system (see Sect.
4.1). The immediate conclusion is that stars as massive as 150 $M_{\odot}$
form in the Local Group beyond reasonable doubt. A similar conclusion was
reached by Tehrani et al. (2019), who reported initial masses of $\approx
150\,M_{\odot}$ for the system Mk 34 – the younger twin of R 144. The
$150\,M_{\odot}$ limit is well within the $200-300\,M_{\odot}$ limit suggested
by the bulk of the R136 population (Schneider et al., 2018b; Gräfener, 2021).

However, as discussed in Sect. 1, the derivation of masses in the range
$200-300\,M_{\odot}$ relies on the assumption that the most massive stars in
the core of the Tarantula region – R 136a1, R 136a2, R 136b, and R 136c – are
single (Crowther et al., 2010; Bestenlehner et al., 2020). Crowther et al.
(2010) argued that lack of substantial X-ray emission in these objects
indicates that they are single. However, as discussed in Sect. 3.7, R 144
provides a clear counter example to this. Future effort should be dedicated
into further investigating the multiplicity of the most massive stars, which
have so-far been largely sensitive only to short-period ($P\lesssim 50\,$d)
binaries (e.g, Schnurr et al., 2009).

###  4.4 Runaway status

Figure 20: 2MASS JHK image (Skrutskie et al., 2006), highlighting the relative
separation of R 144 from R 136 of about 4.1’ ($\approx 60\,$pc). The median-
subtracted proper motion arrow from Gaia eDR3 is also marked.

R 144 is situated in relative isolation with respect to the neighboring very
massive and young cluster R 136, roughly 4.1’ apart (i.e. $\approx 60$ pc,
Fig. 20). As discussed by Sana et al. (2013b), this could either mean that R
144 was formed in situ in relative isolation, or that it is a runaway ejected
from a nearby cluster, most likely R 136. The systemic velocity of
$V_{0}=210\pm 20$km ${\rm s}^{-1}$ measured in our study puts it roughly
$60\,$km ${\rm s}^{-1}$ below the bulk of the Tarantula clusters NGC 2070 and
R 136 (cf. Almeida et al., 2017) of $\approx 270\pm 10\,$km ${\rm s}^{-1}$.
While the measurement of absolute velocities of WR stars is model-dependent
(Sect. 3.1), the observed offset is larger than $3\sigma$ when considering our
conservative error estimate of $20\,$km ${\rm s}^{-1}$.

To further investigate the runaway status of R 144, we retrieved the proper
motion (PM) components of R 144 and all stars within 4’ from the Gaia eDR3
catalogue (Gaia Collaboration et al., 2020). We computed the mean, median, and
standard deviations of the PM components for all objects with a PM accuracy
better than 0.1 ${\rm mas}\,{\rm yr}^{-1}$ in both PM components (464 objects
in total). The standard deviation is $\overrightarrow{\sigma}_{\rm
PM,4^{\prime}}=(0.16,0.17)$ ${\rm mas}\,{\rm yr}^{-1}$ for the RA and DEC
components. The difference between the median and average is negligible. We
then subtracted the PM of R 144 from the median PM vector and obtain
$\overrightarrow{\Delta{\rm PM}}_{\rm R144}=(0.12\pm 0.02,0.28\pm 0.03)$ ${\rm
mas}\,{\rm yr}^{-1}$, where the errors are the Gaia eDR3 measurement errors of
R 144. Evidently, the residual vector has a positive component pointing to the
North-East (away from R 136, Fig. 20). Within errors, the direction of the
residual PM vector supports the idea that R 144 was ejected from R 136 at a
significance larger than 3$\sigma$, with a transverse velocity component of
$\approx 70\,$km ${\rm s}^{-1}$. If R 144 was ejected from R 136 in the past
$\approx 2\,$Myr, a distance of 4.1’ would require a runaway velocity
$\varv\gtrsim 30\,$km ${\rm s}^{-1}$. For a transverse velocity component of
$\approx 70\,$km ${\rm s}^{-1}$, the ejection needs to have happened roughly 1
Myr ago.

Assuming R 144 was indeed ejected from R 136, the ejection was most likely
initiated through dynamical interactions between R 144 and the massive R 136
cluster (Fujii & Portegies Zwart, 2011). Unlike supernovae kicks, which tend
to result in relatively modest ejection velocities (e.g., Renzo et al., 2019),
dynamical interactions of very massive binaries with massive stars in the
cluster can easily result in ejection velocities well above the escape
velocity of the cluster. While single stars are expected to be ejected first,
repeated dynamical interactions over time continuously harden the ”bully
binary” and induce a net momentum to eject it from the cluster. However, it is
worth noting that the dynamical age of the ejection is roughly half our
derived age (2-3 Myr), which is not easily explained in light of the lack of
very massive stars of comparable ages in R 136. We encourage future
investigations of this problem.

In conclusion, the systemic velocity of R 144 and its PM vector in combination
with its relative isolation supports the fact that it was ejected through
dynamical interactions from the neighboring R 136 cluster, rather than having
formed in situ.

##  5 Summary

This study provides a comprehensive spectroscopic, photometric, and
polarimetric analysis of the visually-brightest Wolf-Rayet star in the LMC,
the binary R 144. We relied on multi-epoch ESO spectra acquired with X-SHOOTER
and (FLAMES-)UVES, combined with FUSE and IUE spectra in the UV. We performed
an orbital analysis and spectral disentangling of the system, and used the
non-LTE Potsdam Wolf-Rayet (PoWR) code to derive the physical parameters of
both components. Below, we summarise out conclusions:

  * •

R 144 is a binary comprising two WR stars (WN5/6h + WN6/7h) of nearly equal
masses ($M_{1}\sin^{3}i=48.3\pm 1.8$ and $M_{2}\sin^{3}i=45.5\pm
1.8\,M_{\odot}$) orbiting their center of mass in an eccentric ($e=0.51$)
orbit with a period of $P=74.2$ d. The primary is slightly hotter and more
luminous than the secondary, and both stars exhibit substantial residual
hydrogen ($X_{\rm H}\approx 0.4$) and nitrogen enrichment. Assuming the stars
are still on the main sequence, their inferred age is $\approx 2\,$Myr.

  * •

The systemic velocity, proper-motion vector, and relative isolation of R 144
imply that it was ejected from the nearby R 136 cluster, presumably through
dynamical interactions in the cluster.

  * •

Contemporary evolution models that include enhanced mass-loss rates (Gräfener,
2021) are successful in reproducing the low rotation and significant hydrogen
depletion observed in both stars.

  * •

R 144 shows a periodically modulated light-curve that is well explained with a
hybrid model invoking excess emission from wind-wind collisions and wind
eclipses. The model enables us to derive the orbital inclination with a small
formal error, $i=60.4\pm 1.5^{\circ}$. In turn, this results in accurate
dynamical masses of $74\pm 4\,M_{\odot}$ and $69\pm 4\,M_{\odot}$, making R
144 one of the only very massive binaries in the LMC for which an independent
mass measurement is available. However, these masses are difficult to
reconcile with the luminosities of $\log L=6.44,6.39\,[L_{\odot}]$ for the
primary and secondary, which result in evolutionary masses of the order of
$100-110\,M_{\odot}$. A good agreement between evolutionary and dynamical
masses would be achieved for $i=50^{\circ}$, but this inclination does not
seem to be consistent with the light-curve (however, see discussion in
Appendix D).

  * •

Taken at face value, the derived dynamical masses and luminosities imply that
both components have classical Eddington factors of of $\Gamma_{\rm e}=0.78\pm
0.10$, and are thus expected to be inflated due to their proximity to the
Eddington limit. If the stars are on the main sequence, their derived radii
imply that they are only slightly inflated. Alternatively it is possible that
the components of R 144 are not on the main sequence, but are rather inflated
classical (i.e., core He-burning) WR stars. If so, this would mean that R 144
potentially comprises the most massive classical (i.e., core He-burning) WR
stars to have been weighed so far – the evolved counterpart of the most
massive binary thus weighed, Mk 34 (Tehrani et al., 2019).

R 144 is a one-of-a-kind laboratory to study some of the most urgent questions
at the upper-mass end, and poses several challenges to our theories of stellar
evolution. To advance, we advocate for the acquisition of high-resolution
spectra close to periastron passage, long-term high-precision photometry, and
phase-dependent polarimetry of R 144 and other very massive binaries in the
Local Group.

## Acknowledgments

Dedicated in loving memory of Dr. Simon Clark, who passed away during
preparation of this manuscript. We thank C. Evans for helpful discussions.
This research has made use of the VizieR catalogue access tool, CDS,
Strasbourg, France (DOI : 10.26093/cds/vizier). The original description of
the VizieR service was published in 2000, A&AS 143, 23. TS acknowledges
funding from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme (grant agreement numbers
772225: MULTIPLES). TVR, DMB, and PM gratefully acknowledge support from the
Research Foundation Flanders (FWO) by means of Junior and Senior Postdoctoral
Fellowships, under contract No. 12ZB620N, No. 1286521N, and No. 12ZY520N,
respectively. AFJM is grateful to NSERC (Canada) for financial aid. FRNS has
received funding from the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation programme (Grant agreement No.
945806). LM thanks the European Space Agency (ESA) and the Belgian Federal
Science Policy Office (BELSPO) for their support in the framework of the
PRODEX Programme. S.d.M. was funded in part by the European Union’s Horizon
2020 research and innovation program from the European Research Council (ERC,
Grant agreement No. 715063), and by the Netherlands Organization for
Scientific Research (NWO) as part of the Vidi research program BinWaves with
project number 639.042.728.

##  Appendix A Observation log and RV measurements

Table 2: Observation log and measured radial velocities MJD | phase | Instrument | S/N | $R$ |  ${\rm RV}_{1}$ (N v $\lambda 4945$) [km/s] |  ${\rm RV}_{2}$ (N iii $\lambda\lambda 4634,4641$) [km/s]
---|---|---|---|---|---|---
|  |  |  |  |  km ${\rm s}^{-1}$ |  km ${\rm s}^{-1}$
52956.80 | 0.41 | UVES | 100 | 55000 |  196 $\pm$ 5 |  207 $\pm$ 6
52958.82 | 0.44 | UVES | 150 | 55000 |  179 $\pm$ 5 |  212 $\pm$ 6
55308.98 | 0.11 | X-SHOOTER | 230 | 9700 |  331 $\pm$ 21 | -
55452.39 | 0.04 | X-SHOOTER | 130 | 9700 |  356 $\pm$ 19 |  -3 $\pm$ 22
55452.40 | 0.04 | X-SHOOTER | 210 | 9700 |  373 $\pm$ 21 |  18 $\pm$ 24
55580.08 | 0.76 | X-SHOOTER | 480 | 9700 |  87 $\pm$ 17 |  266 $\pm$ 21
55580.09 | 0.77 | X-SHOOTER | 160 | 9700 |  81 $\pm$ 19 |  255 $\pm$ 21
55585.15 | 0.83 | X-SHOOTER | 480 | 9700 |  96 $\pm$ 19 |  274 $\pm$ 24
55589.18 | 0.89 | X-SHOOTER | 350 | 9700 |  103 $\pm$ 20 |  302 $\pm$ 25
55604.03 | 0.09 | X-SHOOTER | 190 | 9700 |  357 $\pm$ 22 |  21 $\pm$ 24
55672.98 | 0.02 | X-SHOOTER | 170 | 9700 |  334 $\pm$ 19 |  42 $\pm$ 22
55674.48 | 0.04 | X-SHOOTER | 180 | 9700 |  327 $\pm$ 23 |  69 $\pm$ 25
55697.95 | 0.35 | FEROS | 60 | 48000 |  209 $\pm$ 10 |  163 $\pm$ 9
56210.88 | 0.27 | FLAMES-UVES | 110 | 47000 |  224 $\pm$ 6 |  146 $\pm$ 5
56217.85 | 0.36 | FLAMES-UVES | 70 | 47000 |  207 $\pm$ 6 |  181 $\pm$ 6
56243.86 | 0.71 | FLAMES-UVES | 100 | 47000 |  105 $\pm$ 6 |  296 $\pm$ 6
56256.78 | 0.88 | FLAMES-UVES | 90 | 47000 |  111 $\pm$ 7 |  273 $\pm$ 6
56257.65 | 0.90 | FLAMES-UVES | 90 | 47000 |  117 $\pm$ 7 |  291 $\pm$ 6
56277.83 | 0.17 | FLAMES-UVES | 130 | 47000 |  294 $\pm$ 6 |  84 $\pm$ 5
56283.57 | 0.25 | FLAMES-UVES | 70 | 47000 |  238 $\pm$ 7 |  139 $\pm$ 6
56294.72 | 0.40 | FLAMES-UVES | 90 | 47000 |  193 $\pm$ 6 |  204 $\pm$ 6
56295.70 | 0.41 | FLAMES-UVES | 110 | 47000 |  193 $\pm$ 5 |  210 $\pm$ 6
56305.75 | 0.54 | FLAMES-UVES | 100 | 47000 |  158 $\pm$ 6 |  243 $\pm$ 6
56306.74 | 0.56 | FLAMES-UVES | 120 | 47000 |  155 $\pm$ 17 |  234 $\pm$ 6
56316.73 | 0.69 | FLAMES-UVES | 130 | 47000 |  112 $\pm$ 7 |  272 $\pm$ 5
56338.02 | 0.98 | X-SHOOTER | 160 | 9700 |  224 $\pm$ 19 |  126 $\pm$ 22
56349.54 | 0.13 | FLAMES-UVES | 100 | 47000 |  308 $\pm$ 6 |  64 $\pm$ 6
56352.55 | 0.17 | FLAMES-UVES | 80 | 47000 |  289 $\pm$ 7 |  104 $\pm$ 6
56356.53 | 0.23 | FLAMES-UVES | 110 | 47000 |  258 $\pm$ 7 |  123 $\pm$ 6
56571.87 | 0.13 | FLAMES-UVES | 40 | 47000 |  343 $\pm$ 13 |  72 $\pm$ 22
56582.86 | 0.28 | FLAMES-UVES | 10 | 47000 |  212 $\pm$ 22 |  143 $\pm$ 24
56586.77 | 0.33 | FLAMES-UVES | 70 | 47000 |  200 $\pm$ 17 |  190 $\pm$ 6
56597.75 | 0.48 | FLAMES-UVES | 30 | 47000 |  143 $\pm$ 12 |  239 $\pm$ 7
56620.78 | 0.79 | FLAMES-UVES | 80 | 47000 |  105 $\pm$ 12 |  290 $\pm$ 6
56627.69 | 0.88 | FLAMES-UVES | 40 | 47000 |  100 $\pm$ 17 |  248 $\pm$ 15
56645.57 | 0.12 | FLAMES-UVES | 80 | 47000 |  325 $\pm$ 17 |  38 $\pm$ 20
56653.80 | 0.23 | FLAMES-UVES | 80 | 47000 |  231 $\pm$ 42 |  131 $\pm$ 6
56693.63 | 0.77 | FLAMES-UVES | 60 | 47000 |  95 $\pm$ 10 |  306 $\pm$ 6
56697.69 | 0.83 | FLAMES-UVES | 30 | 47000 |  81 $\pm$ 19 |  266 $\pm$ 10
56703.66 | 0.91 | FLAMES-UVES | 70 | 47000 |  111 $\pm$ 12 |  273 $\pm$ 6
56705.04 | 0.92 | X-SHOOTER | 520 | 9700 |  141 $\pm$ 17 |  231 $\pm$ 21
56705.05 | 0.92 | X-SHOOTER | 220 | 9700 |  147 $\pm$ 18 |  232 $\pm$ 21
56706.05 | 0.94 | X-SHOOTER | 310 | 9700 |  166 $\pm$ 18 |  209 $\pm$ 22
56707.06 | 0.95 | X-SHOOTER | 220 | 9700 |  173 $\pm$ 18 |  193 $\pm$ 22
56708.04 | 0.96 | X-SHOOTER | 260 | 9700 |  196 $\pm$ 19 |  188 $\pm$ 21
56709.04 | 0.98 | X-SHOOTER | 810 | 9700 |  215 $\pm$ 18 |  152 $\pm$ 22
56710.03 | 0.99 | X-SHOOTER | 120 | 9700 |  249 $\pm$ 18 |  107 $\pm$ 22
56711.02 | 0.01 | X-SHOOTER | 120 | 9700 |  282 $\pm$ 20 |  74 $\pm$ 22
56712.02 | 0.02 | X-SHOOTER | 160 | 9700 |  331 $\pm$ 17 |  41 $\pm$ 22
56713.03 | 0.03 | X-SHOOTER | 380 | 9700 |  345 $\pm$ 18 |  18 $\pm$ 22
56714.03 | 0.05 | X-SHOOTER | 420 | 9700 |  345 $\pm$ 8 |  23 $\pm$ 22
56714.55 | 0.05 | FLAMES-UVES | 50 | 47000 |  346 $\pm$ 13 |  30 $\pm$ 8
56715.03 | 0.06 | X-SHOOTER | 270 | 9700 |  349 $\pm$ 18 |  36 $\pm$ 22
56716.04 | 0.07 | X-SHOOTER | 190 | 9700 |  338 $\pm$ 18 |  40 $\pm$ 21
56719.54 | 0.12 | FLAMES-UVES | 80 | 47000 |  344 $\pm$ 9 |  79 $\pm$ 6
56723.70 | 0.18 | FLAMES-UVES | 40 | 47000 |  302 $\pm$ 20 | -

##  Appendix B TESS data

Figure 21: The short-cadence TESS light-curve (black dots), obtained for
sectors 27 to 30. Flagged (less reliable) data points are marked in light
gray, and the dashed gray vertical lines indicate the time stamps of the TESS
angular momentum dumps.

Figure 21 shows the entire TESS photometry, which comprises four different
sectors that are stitched together. The data as a whole covers roughly 1.5
orbital cycles and one periastron passage. With a pixel size of 21”,
contamination with nearby sources is possible, though comparison with the GAIA
early DR3 catalogue (Gaia Collaboration, 2018; Gaia Collaboration et al.,
2020) suggests the contamination does not exceed 1-2%. The high-precision TESS
photometry data reveal interesting substructures in the light-curve indicative
of pulsational activity, though the time coverage does not enable a reliable
period analysis of these structure to identify their possible origin.

As Fig. 21 suggests, combining the TESS data with the ASAS-SN is not trivial,
as the stitching and calibration of the individual sectors can lead to large
systematic errors. In fact, as Fig. 21 readily shows, the long-term increase
seen in the ASAS-SN data between $\phi=0.2$ and $0.8$ is seen only to the
right of periastron, but not to the left, implying a substantial systematic
error in the stitched light-curve.

Figure 22: As Fig. 15, but comparing with the phased TESS data instead of the
ASAS-SN data.

Because of these systematics, and since the TESS data only cover 1.5 orbital
cycles and are hence greatly affected by stochasticity in the system, we do
not use the TESS data in our model fitting. However, we do use it to refine
the orbital period of the system, as described in Sect. 3.5. Nevertheless, we
present in Fig. 22 a comparison between the phased TESS light-curve and our
hybrid light-curve model shown in Sect. 3.5. There are some apparent
deviations, but overall the model qualitatively well reproduces the light-
curve. The TESS data seem to suggest a deeper secondary minimum at $\phi=0.28$
(which implies either even higher inclinations or an increased mass-loss rate
for the secondary), but with only 1.5 cycles, we avoid an over-interpretation
of this mismatch.

##  Appendix C PHOEBE modelling of the light-curve

Figure 23: Plotted is our PHOEBE model of R 144, using the parameters given in
Table 3 as input. Due to the different scale of amplitude of the data and
model, the data are not plotted for clarity (cf. Fig. 13).

We attempted to model the ASAS-SN light-curve using the Physics of Eclipsing
Binaries (PHOEBE) light-curve modelling tool (Prša et al., 2016). Starting
from the parameters derived using the radial velocity curves and spectroscopic
fitting we endeavoured to derive a self consistent model. However, our efforts
led to the model seen in Fig 23. While the ”heartbeat” shape is recovered, the
amplitude is nearly two orders of magnitude lower than that seen in the actual
light-curve. Moreover, the minimum is obtained at periastron, while the
observed minimum is at $\phi=-0.026$ (inferior conjunction), suggesting that
tidal distortion is not the dominant mechanism governing the shape of the
light-curve.

##  Appendix D Systematic uncertainties in the light-curve model

Table 3: Inferred parameters from our hybrid light-curve model (2829 fitting
points, combined RV + light-curve).

model |  $\chi^{2}$ a | $A_{\rm WWC}$ | $\gamma_{\rm WWC}$ |  $\dot{M}_{1}$ [$M_{\odot}\,{\rm yr}^{-1}$] |  $\dot{M}_{2}$ [$M_{\odot}\,{\rm yr}^{-1}$] |  $i$ [deg] | $f_{1}/f_{2}(V)$ | $\beta$
---|---|---|---|---|---|---|---|---
adopted | 1.0000 (2827) | $0.12\pm 0.01$ | $1.5\pm 0.1$ | $-4.29\pm 0.05$ | $-4.42\pm 0.05$ | $60.4\pm 1.5$ | 0.79 (fixed) | 1 (fixed)
A2 | 1.0000 (2827) | $0.12\pm 0.01$ | $1.5\pm 0.1$ | $-4.13\pm 0.05$ | $-4.56\pm 0.05$ | $60.4\pm 1.5$ | 1.58 (fixed) | 1 (fixed)
A3 | 1.0000 (2827) | $0.12\pm 0.01$ | $1.5\pm 0.1$ | $-4.39\pm 0.05$ | $-4.22\pm 0.05$ | $60.4\pm 1.5$ | 0.40 (fixed) | 1 (fixed)
A4 | 1.0000 (2827) | $0.12\pm 0.01$ | $1.5\pm 0.1$ | -4.38 (fixed) | -4.34 (fixed) | $61.4\pm 0.3$ | $0.53\pm 0.02$ | 1 (fixed)
A5 | 1.0112 (2860) | $0.19\pm 0.01$ | $1.3\pm 0.1$ | $-3.95\pm 0.01$ | $-4.08\pm 0.01$ | 50 (fixed) | 0.79 (fixed) | 1 (fixed)
A6 | 1.0000 (2827) | $0.11\pm 0.01$ | $1.5\pm 0.1$ | $-4.34\pm 0.05$ | $-4.44\pm 0.04$ | $59.4\pm 1.5$ | 0.75 (fixed) | 2 (fixed)

101010a reduced and non-reduced $\chi^{2}$ values with 12 degrees of freedom
and 2839 fitting points. Since the photometry errors appear to be
overestimated they were reduced by a factor of 0.82 to obtain $\chi^{2}$
values of the order of unity (see Sect. 3.5).

Modelling of the ASAS-SN data enables an accurate derivation of the orbital
inclination of $i=60.4\pm 1.5^{\circ}$. However, the error reflects the formal
fitting error, and does not account for any potential systematic errors
originating in the assumptions in the model. Our numerical attempts to vary
different assumptions are shown in Table 10. Below, we summarise our results
in the context of their impact on the derived value of $i$, and discuss
additional sources of uncertainty and their potential impact.

  1. 1.

_Altering the light ratio:_ we explored the impact of fixing the light ratio
to twice and half the adopted value. This has a substantial impact on the
mass-loss rates, but has no impact on the derived inclination. The fit quality
($\chi^{2}$) remains unchanged, which is anticipated as $\dot{M}_{1}$,
$\dot{M}_{2}$ and $f_{1}/f_{2}$ are degenerate.

  2. 2.

Altering $\dot{M_{1}},\dot{M_{2}}$: similarly, we fixed $\dot{M}$ to the
values derived spectroscopically. This only leads to a slight increase (by
$1^{\circ}$) of the orbital inclination, to compensate for the net decrease of
density in the wind of the primary. The fit quality is comparable.

  3. 3.

_Altering $\beta$_: The $\beta$-exponent in the velocity law determines the
steepness of acceleration in the wind. Larger $\beta$ values result in
enhanced wind densities, and could hence result it higher inclination angles.
Lamontagne et al. (1996) derived analytical equations for $\beta=0$ and
$\beta=1$. However, larger $\beta$ values will result in denser winds, which
in turn should result in lower inclinations for the same column density. we
derive the relevant equations for $\beta=2$ (Appendix E) to explore the impact
of increasing $\beta$. While this indeed yields a slightly lower inclination,
the decrease is small (by $1^{\circ}$) and still falls within the formal error
margin. The fit quality is comparable, suggesting that $\beta$ and $i$ are
degenerate parameters. However, $\beta$ cannot be much larger than $\beta=2$
based on our spectroscopic analysis (Sect. 3.3).

  4. 4.

_Varying $\alpha$, $R_{*}$, or $v_{\infty}$:_ We have altered the number of
electrons per baryon $\alpha$, the stellar radii $R_{*}$, and the terminal
velocities $v_{\infty}$ within formal errors. Some impact, though a small one
($\approx 0.05\,$dex) could be seen on the mass-loss rates, but no notable
change was obtained in the value of $i$.

  5. 5.

_Ionization structure in the wind:_ The model relies on a constant $\alpha$,
assuming that H and He are fully ionized in the wind. Our atmosphere models
reveal that this assumption is valid for the hotter primary in the relevant
domain, but breaks in the cooler secondary a few stellar radii above the
surface, introducing a potential depth-dependence in $\alpha$. However, as
$\alpha$ is primarily a multiplicative parameter determining the strength of
absorption, any changes in $\alpha$ will at first order impact the derived
mass-loss rates from this method, with impact on the inclination being second-
order.

  6. 6.

_Impact of WWC cone:_ At the region where the two winds collide, substantial
density enhancements are expected. Moreover, the cone slices the stellar
winds, thereby altering the optical depth across the line-of-sight in a non-
trivial manner. It is difficult to assess the impact of these deviations from
spherical symmetry without a full 3D radiative transfer model. However, the
impact of this is expected to be very different at inferior and superior
conjunction. The fact that our model shows no notable deviations during both
eclipses suggests that the impact of the WWC cone in this context cannot be
substantial.

It is beyond the scope of this paper to construct a full numerical
hydrodynamic/radiative-transfer model to further investigate additional
sources of uncertainty (e.g., depth-dependent ionization, impact of density
enhancement in WWC cone, finite sizes of stellar disks). Our tests imply that
the derived value of $i$ is robust within the provided errors, but accurate
modelling of the problem in future works is encouraged.

##  Appendix E Wind-eclipse model with $\beta=2$

We provide an extension of the wind-eclipse model of Lamontagne et al. (1996)
for the case of a $\beta$-law velocity exponent of $\beta=2$. To compute the
optical depth from the wind one needs to evaluate the integral

$\displaystyle\tau=k\int_{-\epsilon}^{\infty}\frac{d(z/\mathcal{D})}{(r/\mathcal{})^{2}(1-R_{*}/r)^{\beta}},$
(3)

where

$\displaystyle\begin{aligned} \epsilon=\sin
i\cos\varpi,\quad\left(\frac{r}{\mathcal{D}}\right)^{2}=\cos^{2}i\cos^{2}\varpi+\sin^{2}\varpi+\left(\frac{z}{\mathcal{D}}\right)^{2},\end{aligned}$
(4)

$\mathcal{D}$ is the instantaneous separation between the two components, $z$
is the Cartesian coordinate towards the observer, $r$ is the radial distance
from the center of the eclipsing star (see figure 1 in Lamontagne et al.
1996), and $k$ is the Thomson opacity.

Defining $\gamma=\sqrt{\cos^{2}i\cos^{2}\varpi+\sin^{2}\varpi}$,
$\psi=\sqrt{\gamma^{2}-\delta^{2}}$ and $\delta=R_{*}/a$, the analytical
solution to the integral for the case $\beta=2$ is given by

$\displaystyle\begin{aligned}
\frac{\tau}{k}=&\frac{\delta}{\psi^{2}}\left(1+\frac{\epsilon\sqrt{\epsilon^{2}+\gamma^{2}}+\delta^{2}\epsilon^{2}}{\psi^{2}+\epsilon^{2}}\right)\\\
&+\frac{\gamma^{2}}{\psi^{3}}\left(\arcsin\left(\frac{\delta}{\gamma}\right)+\frac{|\epsilon|}{\epsilon}\arcsin\left(\sqrt{\frac{\delta^{2}\epsilon^{2}}{\gamma^{2}(\psi^{2}+\epsilon^{2})}}\right)\right)\\\
&+\frac{\pi\gamma^{2}}{2\psi^{3}}\left(\frac{\pi}{2}+\arctan\left(\frac{\epsilon}{\psi}\right)\right).\end{aligned}$
(5)

## References

  * Abbott et al. (2020) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020, Phys. Rev. Lett., 125, 101102
  * Almeida et al. (2017) Almeida, L. A., Sana, H., Taylor, W., et al. 2017, A&A, 598, A84
  * Barbá et al. (2010) Barbá, R. H., Gamen, R., Arias, J. I., et al. 2010, in Revista Mexicana de Astronomia y Astrofisica Conference Series, Vol. 38, Revista Mexicana de Astronomia y Astrofisica Conference Series, 30–32
  * Bastian & Lardo (2018) Bastian, N. & Lardo, C. 2018, ARA&A, 56, 83
  * Belczynski et al. (2010) Belczynski, K., Bulik, T., Fryer, C. L., et al. 2010, ApJ, 714, 1217
  * Bestenlehner (2020) Bestenlehner, J. M. 2020, MNRAS, 493, 3938
  * Bestenlehner et al. (2020) Bestenlehner, J. M., Crowther, P. A., Caballero-Nieves, S. M., et al. 2020, MNRAS, 499, 1918
  * Bestenlehner et al. (2014) Bestenlehner, J. M., Gräfener, G., Vink, J. S., et al. 2014, A&A, 570, A38
  * Bonanos et al. (2009) Bonanos, A. Z., Massa, D. L., Sewilo, M., et al. 2009, AJ, 138, 1003
  * Bonanos et al. (2004) Bonanos, A. Z., Stanek, K. Z., Udalski, A., et al. 2004, ApJ, 611, L33
  * Bowman et al. (2019) Bowman, D. M., Burssens, S., Pedersen, M. G., et al. 2019, Nature Astronomy, 3, 760
  * Breysacher et al. (1999) Breysacher, J., Azzopardi, M., & Testor, G. 1999, A&AS, 137, 117
  * Broos et al. (2010) Broos, P. S., Townsley, L. K., Feigelson, E. D., et al. 2010, ApJ, 714, 1582
  * Brott et al. (2011) Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115
  * Brown et al. (1982) Brown, J. C., Aspin, C., Simmons, J. F. L., & McLean, I. S. 1982, MNRAS, 198, 787
  * Brown et al. (1978) Brown, J. C., McLean, I. S., & Emslie, A. G. 1978, A&A, 68, 415
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245
  * Castor et al. (1975) Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157
  * Chené et al. (2008) Chené, A. N., Moffat, A. F. J., & Crowther, P. A. 2008, in Clumping in Hot-Star Winds, ed. W.-R. Hamann, A. Feldmeier, & L. M. Oskinova, 163
  * Cherepashchuk (1976) Cherepashchuk, A. M. 1976, Soviet Astronomy Letters, 2, 138
  * Cleveland (1979) Cleveland, W. S. 1979, Journal of the American Statistical Association, 74, 829
  * Corcoran et al. (1996) Corcoran, M. F., Stevens, I. R., Pollock, A. M. T., et al. 1996, ApJ, 464, 434
  * Crowther & Dessart (1998) Crowther, P. A. & Dessart, L. 1998, MNRAS, 296, 622
  * Crowther et al. (2010) Crowther, P. A., Schnurr, O., Hirschi, R., et al. 2010, MNRAS, 408, 731
  * Crowther & Smith (1997) Crowther, P. A. & Smith, L. J. 1997, A&A, 320, 500
  * Crowther & Walborn (2011) Crowther, P. A. & Walborn, N. R. 2011, MNRAS, 416, 1311
  * Cutri & et al. (2014) Cutri, R. M. & et al. 2014, VizieR Online Data Catalog, II/328
  * David-Uraz et al. (2012) David-Uraz, A., Moffat, A. F. J., Chené, A.-N., et al. 2012, MNRAS, 426, 1720
  * de Koter et al. (1997) de Koter, A., Heap, S. R., & Hubeny, I. 1997, ApJ, 477, 792
  * Demers et al. (2002) Demers, H., Moffat, A. F. J., Marchenko, S. V., Gayley, K. G., & Morel, T. 2002, ApJ, 577, 409
  * Doran et al. (2013) Doran, E. I., Crowther, P. A., de Koter, A., et al. 2013, A&A, 558, A134
  * Dsilva et al. (2020) Dsilva, K., Shenar, T., Sana, H., & Marchant, P. 2020, arXiv e-prints, arXiv:2006.13957
  * Evans et al. (2011) Evans, C. J., Taylor, W. D., Hénault-Brunet, V., et al. 2011, A&A, 530, A108
  * Figer (2005) Figer, D. F. 2005, Nature, 434, 192
  * Figer et al. (2002) Figer, D. F., Najarro, F., Gilmore, D., et al. 2002, ApJ, 581, 258
  * Fryer et al. (2001) Fryer, C. L., Woosley, S. E., & Heger, A. 2001, ApJ, 550, 372
  * Fujii & Portegies Zwart (2011) Fujii, M. S. & Portegies Zwart, S. 2011, Science, 334, 1380
  * Fullerton et al. (2006) Fullerton, A. W., Massa, D. L., & Prinja, R. K. 2006, ApJ, 637, 1025
  * Gaia Collaboration (2018) Gaia Collaboration. 2018, VizieR Online Data Catalog, I/345
  * Gaia Collaboration et al. (2020) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2020, arXiv e-prints, arXiv:2012.01533
  * Gieles et al. (2018) Gieles, M., Charbonnel, C., Krause, M. G. H., et al. 2018, MNRAS, 478, 2461
  * Gilkis et al. (2021) Gilkis, A., Shenar, T., Ramachandran, V., et al. 2021, arXiv e-prints, arXiv:2102.03102
  * Giménez-García et al. (2016) Giménez-García, A., Shenar, T., Torrejón, J. M., et al. 2016, A&A, 591, A26
  * González & Levato (2006) González, J. F. & Levato, H. 2006, A&A, 448, 283
  * Gräfener (2021) Gräfener, G. 2021, arXiv e-prints, arXiv:2101.03837
  * Gräfener & Hamann (2008) Gräfener, G. & Hamann, W. R. 2008, A&A, 482, 945
  * Gräfener et al. (2002) Gräfener, G., Koesterke, L., & Hamann, W. R. 2002, A&A, 387, 244
  * Gräfener et al. (2012) Gräfener, G., Owocki, S. P., & Vink, J. S. 2012, A&A, 538, A40
  * Gräfener et al. (2011) Gräfener, G., Vink, J. S., de Koter, A., & Langer, N. 2011, A&A, 535, A56
  * Grassitelli et al. (2016) Grassitelli, L., Chené, A. N., Sanyal, D., et al. 2016, A&A, 590, A12
  * Guerrero & Chu (2008) Guerrero, M. A. & Chu, Y.-H. 2008, ApJS, 177, 216
  * Hadrava (1995) Hadrava, P. 1995, A&AS, 114, 393
  * Hainich et al. (2014) Hainich, R., Rühling, U., Todt, H., et al. 2014, A&A, 565, A27
  * 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
  * Higgins & Vink (2019) Higgins, E. R. & Vink, J. S. 2019, A&A, 622, A50
  * Hillier (1991) Hillier, D. J. 1991, A&A, 247, 455
  * Howarth (1983) Howarth, I. D. 1983, MNRAS, 203, 301
  * Jenkins et al. (2016) Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9913, Software and Cyberinfrastructure for Astronomy IV, ed. G. Chiozzi & J. C. Guzman, 99133E
  * Jermyn et al. (2018) Jermyn, A. S., Tout, C. A., & Chitre, S. M. 2018, MNRAS, 480, 5427
  * Kaufer et al. (1997) Kaufer, A., Wolf, B., Andersen, J., & Pasquini, L. 1997, The Messenger, 89, 1
  * Kochanek et al. (2017) Kochanek, C. S., Shappee, B. J., Stanek, K. Z., et al. 2017, PASP, 129, 104502
  * Köhler et al. (2015) Köhler, K., Langer, N., de Koter, A., et al. 2015, A&A, 573, A71
  * Kozai (1962) Kozai, Y. 1962, AJ, 67, 591
  * Lamontagne et al. (1996) Lamontagne, R., Moffat, A. F. J., Drissen, L., Robert, C., & Matthews, J. M. 1996, AJ, 112, 2227
  * Langer et al. (2007) Langer, N., Norman, C. A., de Koter, A., et al. 2007, A&A, 475, L19
  * Larson & Starrfield (1971) Larson, R. B. & Starrfield, S. 1971, A&A, 13, 190
  * Lépine & Moffat (1999) Lépine, S. & Moffat, A. F. J. 1999, ApJ, 514, 909
  * Lidov (1962) Lidov, M. L. 1962, Planet. Space Sci., 9, 719
  * Lightkurve Collaboration et al. (2018) Lightkurve Collaboration, Cardoso, J. V. d. M., Hedges, C., et al. 2018, Lightkurve: Kepler and TESS time series analysis in Python
  * Lohr et al. (2018) Lohr, M. E., Clark, J. S., Najarro, F., et al. 2018, A&A, 617, A66
  * Magalhaes et al. (1984) Magalhaes, A. M., Benedetti, E., & Roland, E. H. 1984, PASP, 96, 383
  * Mahy et al. (2020a) Mahy, L., Almeida, L. A., Sana, H., et al. 2020a, A&A, 634, A119
  * Mahy et al. (2020b) Mahy, L., Sana, H., Abdul-Masih, M., et al. 2020b, A&A, 634, A118
  * Maíz Apellániz et al. (2017) Maíz Apellániz, J., Sana, H., Barbá, R. H., Le Bouquin, J. B., & Gamen, R. C. 2017, MNRAS, 464, 3561
  * Maíz Apellániz et al. (2019) Maíz Apellániz, J., Trigueros Páez, E., Jiménez Martínez, I., et al. 2019, in Highlights on Spanish Astrophysics X, ed. B. Montesinos, A. Asensio Ramos, F. Buitrago, R. Schödel, E. Villaver, S. Pérez-Hoyos, & I. Ordóñez-Etxeberria, 420–420
  * Marchant et al. (2019) Marchant, P., Renzo, M., Farmer, R., et al. 2019, ApJ, 882, 36
  * Marchenko et al. (2003) Marchenko, S. V., Moffat, A. F. J., Ballereau, D., et al. 2003, ApJ, 596, 1295
  * Marchenko et al. (1998) Marchenko, S. V., Moffat, A. F. J., & Eenens, P. R. J. 1998, PASP, 110, 1416
  * Martins et al. (2008) Martins, F., Hillier, D. J., Paumard, T., et al. 2008, A&A, 478, 219
  * Moffat (1989) Moffat, A. F. J. 1989, ApJ, 347, 373
  * Moffat et al. (1988) Moffat, A. F. J., Drissen, L., Lamontagne, R., & Robert, C. 1988, ApJ, 334, 1038
  * Moffat et al. (1998) Moffat, A. F. J., Marchenko, S. V., Seggewiss, W., et al. 1998, A&A, 331, 949
  * Monet et al. (1998) Monet, D., Canzian, B., Harris, H., et al. 1998, VizieR Online Data Catalog, I/243
  * Najarro et al. (2004) Najarro, F., Figer, D. F., Hillier, D. J., & Kudritzki, R. P. 2004, ApJ, 611, L105
  * Nazé et al. (2007) Nazé, Y., Corcoran, M. F., Koenigsberger, G., & Moffat, A. F. J. 2007, ApJ, 658, L25
  * Newville et al. (2014) Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2014, LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python
  * Oey & Clarke (2005) Oey, M. S. & Clarke, C. J. 2005, ApJ, 620, L43
  * Oskinova et al. (2007) Oskinova, L. M., Hamann, W. R., & Feldmeier, A. 2007, A&A, 476, 1331
  * Parker (1992) Parker, J. W. 1992, PASP, 104, 1107
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3
  * Pedersen et al. (2018) Pedersen, M. G., Aerts, C., Pápics, P. I., & Rogers, T. M. 2018, A&A, 614, A128
  * Petrovic et al. (2006) Petrovic, J., Pols, O., & Langer, N. 2006, A&A, 450, 219
  * Pojmanski (2002) Pojmanski, G. 2002, Acta Astron., 52, 397
  * Pollock (1987) Pollock, A. M. T. 1987, ApJ, 320, 283
  * Pollock et al. (2018) Pollock, A. M. T., Crowther, P. A., Tehrani, K., Broos, P. S., & Townsley, L. K. 2018, MNRAS, 474, 3228
  * Prša et al. (2016) Prša, A., Conroy, K. E., Horvat, M., et al. 2016, ApJS, 227, 29
  * Puls et al. (2006) Puls, J., Markova, N., Scuderi, S., et al. 2006, A&A, 454, 625
  * Ramachandran et al. (2019) Ramachandran, V., Hamann, W. R., Oskinova, L. M., et al. 2019, A&A, 625, A104
  * Ramiaramanantsoa et al. (2019) Ramiaramanantsoa, T., Ignace, R., Moffat, A. F. J., et al. 2019, MNRAS, 490, 5921
  * Rauw et al. (2004) Rauw, G., De Becker, M., Nazé, Y., et al. 2004, A&A, 420, L9
  * Rauw et al. (1999) Rauw, G., Vreux, J. M., & Bohannan, B. 1999, ApJ, 517, 416
  * Renzo et al. (2019) Renzo, M., Zapartas, E., de Mink, S. E., et al. 2019, A&A, 624, A66
  * Richardson et al. (2017) Richardson, N. D., Russell, C. M. P., St-Jean, L., et al. 2017, MNRAS, 471, 2715
  * Richardson et al. (2016) Richardson, N. D., Shenar, T., Roy-Loubier, O., et al. 2016, MNRAS, 461, 4115
  * Ricker et al. (2015) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003
  * Sana et al. (2013a) Sana, H., de Koter, A., de Mink, S. E., et al. 2013a, A&A, 550, A107
  * Sana et al. (2013b) Sana, H., van Boeckel, T., Tramper, F., et al. 2013b, MNRAS, 432, L26
  * Sander et al. (2015) Sander, A., Shenar, T., Hainich, R., et al. 2015, A&A, 577, A13
  * Sanyal et al. (2017) Sanyal, D., Langer, N., Szécsi, D., -C Yoon, S., & Grassitelli, L. 2017, A&A, 597, A71
  * Schmutz et al. (1989) Schmutz, W., Hamann, W. R., & Wessolowski, U. 1989, A&A, 210, 236
  * Schneider et al. (2014) Schneider, F. R. N., Langer, N., de Koter, A., et al. 2014, A&A, 570, A66
  * Schneider et al. (2018a) Schneider, F. R. N., Ramírez-Agudelo, O. H., Tramper, F., et al. 2018a, A&A, 618, A73
  * Schneider et al. (2018b) Schneider, F. R. N., Sana, H., Evans, C. J., et al. 2018b, Science, 359, 69
  * Schnurr et al. (2008a) Schnurr, O., Casoli, J., Chené, A. N., Moffat, A. F. J., & St-Louis, N. 2008a, MNRAS, 389, L38
  * Schnurr et al. (2009) Schnurr, O., Chené, A. N., Casoli, J., Moffat, A. F. J., & St-Louis, N. 2009, MNRAS, 397, 2049
  * Schnurr et al. (2008b) Schnurr, O., Moffat, A. F. J., St-Louis, N., Morrell, N. I., & Guerrero, M. A. 2008b, MNRAS, 389, 806
  * Schootemeijer & Langer (2018) Schootemeijer, A. & Langer, N. 2018, A&A, 611, A75
  * Schootemeijer et al. (2019) Schootemeijer, A., Langer, N., Grin, N. J., & Wang, C. 2019, A&A, 625, A132
  * Schwarzenberg-Czerny (1997) Schwarzenberg-Czerny, A. 1997, ApJ, 489, 941
  * Shappee et al. (2014) Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, ApJ, 788, 48
  * Shenar et al. (2020) Shenar, T., Bodensteiner, J., Abdul-Masih, M., et al. 2020, A&A, 639, L6
  * Shenar et al. (2018) Shenar, T., Hainich, R., Todt, H., et al. 2018, A&A, 616, A103
  * Shenar et al. (2014) Shenar, T., Hamann, W. R., & Todt, H. 2014, A&A, 562, A118
  * Shenar et al. (2015) Shenar, T., Oskinova, L., Hamann, W. R., et al. 2015, ApJ, 809, 135
  * Shenar et al. (2017a) Shenar, T., Oskinova, L. M., Järvinen, S. P., et al. 2017a, A&A, 606, A91
  * Shenar et al. (2017b) Shenar, T., Richardson, N. D., Sablowski, D. P., et al. 2017b, A&A, 598, A85
  * Shenar et al. (2019) Shenar, T., Sablowski, D. P., Hainich, R., et al. 2019, A&A, 627, A151
  * Simmons & Boyle (1984) Simmons, J. F. L. & Boyle, C. B. 1984, A&A, 134, 368
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
  * Smith et al. (1996) Smith, L. F., Shara, M. M., & Moffat, A. F. J. 1996, MNRAS, 281, 163
  * St. -Louis et al. (1988) St. -Louis, N., Moffat, A. F. J., Drissen, L., Bastien, P., & Robert, C. 1988, ApJ, 330, 286
  * St-Louis et al. (2005) St-Louis, N., Moffat, A. F. J., Marchenko, S., & Pittard, J. M. 2005, ApJ, 628, 953
  * Stellingwerf (1978) Stellingwerf, R. F. 1978, ApJ, 224, 953
  * Stevens et al. (1992) Stevens, I. R., Blondin, J. M., & Pollock, A. M. T. 1992, ApJ, 386, 265
  * Sundqvist et al. (2010) Sundqvist, J. O., Puls, J., & Feldmeier, A. 2010, A&A, 510, A11
  * Tehrani (2019) Tehrani, K. 2019, PhD thesis, University of Sheffield.
  * Tehrani et al. (2019) Tehrani, K. A., Crowther, P. A., Bestenlehner, J. M., et al. 2019, MNRAS, 484, 2692
  * Thompson et al. (2012) Thompson, S. E., Everett, M., Mullally, F., et al. 2012, ApJ, 753, 86
  * Todt et al. (2015) Todt, H., Sander, A., Hainich, R., et al. 2015, A&A, 579, A75
  * Toonen et al. (2020) Toonen, S., Portegies Zwart, S., Hamers, A. S., & Bandopadhyay, D. 2020, A&A, 640, A16
  * Townsley et al. (2011) Townsley, L. K., Broos, P. S., Corcoran, M. F., et al. 2011, ApJS, 194, 1
  * Tramper et al. (2016) Tramper, F., Sana, H., Fitzsimons, N. E., et al. 2016, MNRAS, 455, 1275
  * Usov (1992) Usov, V. V. 1992, ApJ, 389, 635
  * Vigna-Gómez et al. (2021) Vigna-Gómez, A., Toonen, S., Ramirez-Ruiz, E., et al. 2021, ApJ, 907, L19
  * Vink (2018) Vink, J. S. 2018, A&A, 615, A119
  * Vink et al. (2001) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574
  * Vink & Harries (2017) Vink, J. S. & Harries, T. J. 2017, A&A, 603, A120
  * Vink et al. (2011) Vink, J. S., Muijres, L. E., Anthonisse, B., et al. 2011, A&A, 531, A132
  * Šurlan et al. (2013) Šurlan, B., Hamann, W. R., Aret, A., et al. 2013, A&A, 559, A130
  * Williams et al. (1997) Williams, P. M., Dougherty, S. M., Davis, R. J., et al. 1997, MNRAS, 289, 10
  * Willis et al. (2004) Willis, A. J., Crowther, P. A., Fullerton, A. W., et al. 2004, ApJS, 154, 651
  * Woosley et al. (2007) Woosley, S. E., Blinnikov, S., & Heger, A. 2007, Nature, 450, 390
  * Zahn (1977) Zahn, J. P. 1977, A&A, 500, 121
  * Zucker & Mazeh (1994) Zucker, S. & Mazeh, T. 1994, ApJ, 420, 806