EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 650 (June 2021) A&A, 650 (2021) A147 Full HTML
Free Access
Issue
A&A
Volume 650, June 2021
Article Number A147
Number of page(s) 21
Section Stellar atmospheres
DOI https://doi.org/10.1051/0004-6361/202140693
Published online 21 June 2021

© ESO 2021

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 ≈ 120 M⊙ to a few thousands of solar masses, depending on the metallicity and modelling 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). Very massive stars 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. de Koter et al. (1997) found that stars with current masses in the excess of M ≈ 100 M⊙ 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 200–300 M⊙ (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 ≈ 0.5 Z⊙) 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 K1 and K2, the minimum masses M1 sin3 i and M2 sin3 i can be derived from Newtonian mechanics; we note 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 M1 = 116 ± 31 M⊙ 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 in which 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 decreases. Second, a calibration of one of the components automatically results in the mass of the second component, which allows us 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: M1 = 139−18+21$M_1\,{=}\,139^{+21}_{-18}$ M⊙ and M2 = 127 ± 17 M⊙.

With a Smith visual magnitude of v = 11.15 mag (Breysacher et al. 1999), R 144 (aka 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-infrared (NIR, 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-linedspectroscopic 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⊙) ≈ 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 LocalGroup.

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 because this binary 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 focussed 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 (ESO) Very Large Telescope (VLT), which is described in Sect. 2. Our analysis is presented in Sect. 3; the goal of this analysis is to derive its orbit, disentangle its composite spectra, assess the dynamical masses of both binary components, and derive 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

Our main spectroscopic dataset includes three different instruments spanning roughly four years. The first set of data was obtained with X-shooter through the Dutch programme for guaranteed observing time from January 2011 to February 2013, allowing thedetection 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 October 2012 to March 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 set-up, providing a wavelength coverage from 4200 to 6800 Å at a spectral resolving power of 47 000. Each observation was composed of three consecutive 980 s 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 fibres 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 days time base. The data were obtained in nodding mode; two nodding cycles were completed to guarantee optimal sky subtraction. The nod-throw was 5" and a jitter of 1" around each nodding position was included. The NDIT by DITs of 1 × 150 s and 2 × 85 s were performed at each nodding position for, respectively, the UVB/VIS arms and NIR arm. Narrow slit widths of 0.5", 0.4", and 0.4" for the UVB, VIS, and NIR spectra were adopted, delivering spectral resolving powers of 9700, 18 400, and 11 600, respectively.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 signal-to-noise ratio (S/N), after combining all nodding exposures of a given arm, is in the range 250–350.

An additional spectrum is obtained with the Fiber-fed Extended Range Optical Spectrograph (FEROS) instrument (R = 48 000, S∕N ≈ 60, Kaufer et al. 1997) mounted on the MPG/ESO 2.20 m telescope in La Silla and is 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 November 2003 with the DIC1 390+564 instrument set-up and an entrance slit of 0.7", delivering a resolving power of ~ 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 Å 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 ≈ 20 000 and S∕N ≈ 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 Åwith R ≈ 10 000 and S∕N ≈ 15. The second was acquired on 1979 Feb 14 (PI:Savage, ID: LWR03766) and covers 2000–3000 Åwith R ≈ 10 000 and S∕N ≈ 10. All three spectra are flux-calibrated and are thus useful to fit the spectral energy distribution (SED) of R 144 and to derive the wind parameters of both components.

thumbnail Fig. 1

X-shooter spectrum of R 144 taken on 20 February 2014 (MJD = 56708.04). An overview of the optical part of the spectrum is shown.

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 because it offers a broader bandwidth and the signal in this filter is more pronounced, covering a baseline of 1609 days (roughly 22 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 and instrumental peaks that 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 MAST, which span a total of 108 days. We extracted a customized light curve of R 144 from the 2 min TESS postage stamp pixel data using a 1σ threshold binary aperture mask, as defined in the lightkurve Python package (Lightkurve Collaboration 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 VizieR1: 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).

thumbnail Fig. 2

Three X-shooter spectra (MJD = 55452.89, 55585.65, 56708.04, see legend) corresponding to phases close to the RV extremes (ϕ = 0.05, 0.84 with ephemeris derived in this study) and the systemic velocity (ϕ = 0.96) of the N IV λ4058 and N V λ4945 lines. Both binary components are seen in the N IVλ4058 line, but only the hotter primary is seen in the N Vλ4945 line. Conversion to Doppler space is with reference to rest wavelengths.

2.3 X-ray data

X-ray observations of R 144 in the range ≈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). Pollock et al. (2018) and Townsley et al. (2011) give a detailed account of the observing programme 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

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λ4058 and N V λ4945 lines revealsthe 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 λ4945 line, which forms predominantly in one of the components. Another example is shown in Fig. 3 for the N V λλ4604, 4620 doublet and the N III λλ4634, 4641 triplet. Both stars appear to contribute in both line blends, but the indicative anti-phase behaviour of these lines implies thateach of them is dominated by a different stellar component. The dynamical spectra of the N V λ4945, N III λλ4634, 4641, and N V λλ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.

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 λ4945 line (see Fig. 2). For the (cooler) secondary, we use the peak of the N III λλ4634, 4641 triplet (see Fig. 3), as it is strongly dominated by the secondary. The template of the primary, focussed on the N V λ4945 line, is calibratedto 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 ≈ 20 km s−1 on the absolute value from comparison to models. The final measurements are given in Table A.1.

Relying on the RV measurements, we derive an orbital solution, constraining the time of periastron T0, eccentricity e, argument of periastron ω, systemic velocity V0, and RV amplitudes K1 and K2. The orbital period is fixed to P = 74.2074 days 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 lmfit2 (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 1, and the best-fitting orbital solution is shown in Fig. 5. We obtain a root mean square error (rms) of 15 km s−1 for both components. In contrast, the formal measurement errors have a mean of σ ≈ 10 km s−1. It is, however, not uncommon for formal errors of RV measurements to be underestimated in the case of WR stars, primarily owing to spectral line variability that stems from either intrinsic wind variability or wind-wind collisions (WWCs) 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 more realistically reflect the true measurement errors.

The inferred orbital parameters are given in Table 1. The two stars are found to have almost identical RV amplitudes and hence masses; there is an indication that the hotter companion is 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 sin3 i = 48.3 M⊙ and 45.5 M⊙, respectively.

thumbnail Fig. 3

As Fig. 2, but showing three FLAMES-UVES spectra (MJD = 56645.57, 56645.57, 56697.69), and focussing on the N V λλ4604, 4620 and N III λλ4634, 4641 lines. Both components contribute to the N V λλ4604, 4620 lines, while the cooler secondary dominates the N IIIλλ4634, 4641 lines.
thumbnail Fig. 4

Dynamical spectra for the N Vλ4945, N III λλ4634, 4641, and N V λλ4604, 4620 lines. The red and green curves depict the orbital solution derived in Sect. 3.1 for the primary and secondary components, respectively.
thumbnail Fig. 5

Radial-velocity measurements of the N Vλ4945 line (primary) and N IIIλλ4634, 4641 triplet (secondary),along with the RV curves corresponding to the orbital parameters given in Table 1.

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. 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 the orbital parameters were derived in Sect. 3.1, we can in principle fix these parameters to separate the spectra. However, by allowing K1 and K2 to vary, we canuse additional line blends to verify the results obtained in Sect. 3.1. We do this by calculating the χ2 of the residuals obtained when shifting-and-adding the disentangled spectra and subtracting the resulting spectrum from the individual observationsiteratively (see e.g. Shenar et al. 2020).

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

The disentangled X-shooter spectra are shown in Fig. 8. The spectra can be disentangled up to a scaling factor that depends on the optical light ratio of the two stars. For reasons that we outline in Sect. 3.3, we adopt the optical light ratio of f1∕f2(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. WWCs), 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 very 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 ionization stages N. 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 v∞ 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 spectrapresented 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.

Table 1

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

thumbnail Fig. 6

Upper panels: comparison between the disentangled spectra obtained for the best-fitting K1 and K2 values, their sum, and observations at two phases close to RV extremes for the N IV λ4058. Lower panel:reduced χ2 as a functionof K1 and K2 for the N IV λ4058 line. The analysis is performed using the X-shooter data. The minimum is at K1 = 129 ± 6 km s−1, K2 = 134 ± 4 km s−1.

3.3 Spectral analysis

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). A 1D code, PoWR solves the radiative transfer problem in spherical geometry. This code 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 WWCs, mutual irradiation, or other non-spherical effects, these effects are globally negligible. They can however 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 grids3 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): XP = 2.91 × 10−6, XSi = 3.21 × 10−4, XS = 1.55 × 10−4, XFe = 7.02 × 10−4; the rest are fitting parameters. The goodness of the fit is judged through visual inspection; producing χ2 -like estimates in the framework of full non-LTE calculations, especially in a binary, is not feasible and would in any case neglect systematic errors originating in the uncertain structure of the winds, which are the dominant factor of uncertainty in this work.

The main parameters defining a model atmosphere are the effective temperature T*, the luminosity L, the mass-loss rate Ṁ, and the chemical abundances, which are expressed as mass fractions in this work. Because of the optically thick winds of WR stars, their photospheres (defined at a mean optical depth of τross ≈ 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 τRoss = 204. The stellar radius R* is related to T* and L via the Stefan–Boltzmann equation L/L⊙ = (R/R⊙)2 (T/T⊙)4$L/L_{\odot}\,{=}\,(R/R_{\odot})^2\,(T/T_{\odot})^4$. The wind velocity v(r) assumes the functional form of a β-law (Castor et al. 1975), characterised by a unitlessparameter β of the order of unity and the terminal velocity v∞. We exploredvarious β values in the range 1–10 and can conclude that β ≲ 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 β = 1, lacking evidence to assume otherwise. The Doppler widths of the opacity and emissivity profiles are defined as vDop = vth2+ξ2$\varv_{\textrm{Dop}}\,{=}\,\sqrt{\varv_{\textrm{th}}^2 + \xi^2}$, where vth (r) is the thermal motion of the chemicalspecies and ξ(r) is a prespecified microturbulence (Shenar et al. 2015). We assume ξ(R*) = 30 km s−1 and that ξ(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 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 ρ (P-Cygni lines, electron-scattering wings) and ρ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 in here. While we explored various clumping parameters up to D = 100, we cannot derive D unambiguously, and hence we 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 M˙D$\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 Ṁ, but also on the size of the stellar surface (∝R*2$\propto R_*^2$). The transformed radius Rt (Schmutz et al. 1989), defined as Rt = R*[v∞2500 km s−1 /M˙D10−4 M⊙ yr−1]2/3,\begin{equation*} R_{\text{t}}\,{=}\,R_* \left[ \frac{v_{\infty}}{2500\,\textrm{km}\,\textrm{s}^{-1}\,} \middle/ \frac{\dot{M} \sqrt{D}}{10^{-4}\,M_{\odot}\,\textrm{yr}^{-1}} \right]^{2/3},\end{equation*}(1)

provides a useful parameter that preserves the strength of emission recombination lines in the model. Models with the same values of T*, Rt, and chemical abundances would exhibit emission lines of comparable equivalent width, irrespective of D, Ṁ, v∞, 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 XH of the stars. Fixing XH to the values derived in our analysis (see below), we find that the luminosity of the hotter primary should be larger byabout 0.05 dex. Alternatively, if we assume that both stars are comprised 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 and56% for the cooler secondary, which we adopt in this work.

The effective temperatures are determined from the balance of lines belonging to N III, IV, V. The wind parameters Rt and v∞ are determined from the strengths and shapes of recombination lines and P-Cygni lines. The mass fractions XN, XH, and XC 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 EB−V,Gal = 0.03 mag and RV = 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: EB−V = EB−V,Gal + EB−V,LMC = 0.03 + 0.17 = 0.20 mag,

A comparison between the final models and observations is shown in Fig. 9, and our derived parameters are given in Table 1. 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 these lines consistently appears to come at a cost of reproducing diagnostic lines such as N V λλ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 λλ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 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). 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 XH 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 disentangledspectra (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 v sin i using emission lines that form close to the stellar surface, such as the N V λ4945 and N IV λ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 v sin i calculations for the primary component and observations is shown in Fig. 11, focussing on the N V λ4945 line. Evidently, v sini ≈ 60 km s−1 best reproduces the observed profile. However, there is a significant degeneracy between the adopted microturbulence ξ and v sin i, which prevents us from providing actual estimates for vsin i. Nevertheless, upper limits can be obtained. For the secondary component, the N IV λ4058 line was used. These upper limits are provided in Table 1.

The uncertainties on T*, Rt, logL, EB−V, XH, XC, XN, XO 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.

thumbnail Fig. 7

As Fig. 6, but for N V λλ4604, 4620. The analysis is performed using the FLAMES-UVES data and X-shooter data. The minimum is at K1 = 130 ± 3 km s−1, K2 = 134 ± 4 km s−1.
thumbnail Fig. 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).

3.4 Wind-wind collision

When two dense winds in a massive binary collide, they form a WWC cone; its tip is located at the region at which 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. 12the 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 λλ4604, 4620 doublet. The WWC emission may increase the line flux by up to ≈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 λ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 ≈ 10%, as we observe (Sect. 3.5).

thumbnail Fig. 9

Comparison between observed SED (blue squares and lines, upper panel) and normalized FUSE, IUE, and X-shooter (ϕ = 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 indicated in red.
thumbnail Fig. 10

Schematic of R 144 as seen during primary inferior conjunction (ϕ ≈ −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 WWC cone. The line-forming regions of several diagnostic lines for the primary star, as calculated from our tailored model atmosphere (Sect. 3.3), are denoted. The line-forming regions of the secondary star are comparable.
thumbnail Fig. 11

Comparison between an observed spectrum of the N Vλ4945 multiplet and model spectra for the primary calculated with the parameters in Table 1, 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).
thumbnail Fig. 12

Variations of EW with phase for He Iλ4472 (4465–4490 Å), He Iλ5876 (5875–5900 Å), He IIλ4686 (4670–4720 Å), Hβ (4840–4890 Å), N IIIλλ4634, 4641 (4620–4650 Å), and N V λλ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.
thumbnail Fig. 13

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

3.5 Visual light curve

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 ± 0.0043 days (Table 1). 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 1 and binned at Δϕ = 0.002. The light curve mimics that of heartbeat stars (Thompson et al. 2012) – a class of eccentric binaries that show rapid periodic variations intheir light curve, which are typically associated with tidal distortion of the star(s) during periastron passage. We attempted to model the ASAS-SN light curve with 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; this suggests 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 in which adiabatic conditions prevail, this emission is expected to grow roughly as EWWC∝1/D(ϕ),$E_{\textrm{WWC}} \propto 1/\mathcal{D}(\phi),$ reaching a maximum at periastron (ϕ = 0), where D(ϕ)$\mathcal{D}(\phi)$ is the instantaneous separation between the stars (Usov 1992). We model the excess emission stemming in WWC as EWWC(ϕ) = AWWC⋅(DminD(ϕ))γWWC = AWWC⋅(1+e1+e cos ν)γWWC,\begin{equation*} E_{\textrm{WWC}} (\phi)\,{=}\,A_{\textrm{WWC}}\cdot \left(\frac{\mathcal{D}_{\textrm{min}}}{\mathcal{D}(\phi)}\right)^{\gamma_{\textrm{WWC}}}\,{=}\,A_{\textrm{WWC}}\cdot \left(\frac{1 + e}{1 + e\,\cos\,\nu}\right)^{\gamma_{\textrm{WWC}}}, \end{equation*}(2)

where ν(ϕ) is the true anomaly, AWWC is the maximum excess at periastron, and γWWC is a number of the order of unity. The quantities AWWC and γ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 (ϕ = −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. In this work, we extend this model to elliptical orbits by substituting the angle 2πϕ with the angle5 ϖ = ν + ω + π∕2, defined to be 0 when the primary eclipses the secondary (i.e. at ν = −ω − π∕2), and substituting the constant separation a in the case of a circular orbit with D(ϖ)$\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 1. Aside from the orbital parameters P, e, ω, T0, the hybrid model contains six fitting parameters: AWWC, γWWC, Ṁ1, Ṁ2, i, and the visual light ratio f1 ∕f2(V). However, we cannot fit the light ratio simultaneously to both Ṁ1, Ṁ2 as the model becomes degenerate. Therefore, we fix the light ratio to the value estimated in Sect. 3.3 (f1 ∕f2 = 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 ≈ 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 β law. Analytical solutions are provided by Lamontagne et al. (1996) for the β = 0 (constant velocity) and β = 1, but arbitrary values of β may be assumed if the integration is performed numerically. For now, we adopt β = 1, but explore β = 2 in Appendices D and E. (v) The wind is fully ionized, which yields the number of electrons per baryon α ≃(1 + XH)∕2, resulting in α1 = 0.68 and α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 α 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.

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 for two reasons. First, the time span of the TESS light curve only covers 1.5 orbital periods, that is only one wind-eclipse is seen. Second, 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 AWWC = 0.12 ± 0.01, γWWC = 1.5 ± 0.1, log Ṁ1 = − 4.29 ± 0.05 M⊙ yr−1, log Ṁ2 = −4.42 ± 0.05 M⊙ yr−1, and i = 60.4° (given also in Table D.1).

The agreement between the ASAS-SN data and our analytical model is excellent; a formal reduced χ2 for the RV curves plus light curve combined is χ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 χ2 values of the order of unity, which ensures that the relative weighting between the photometric and RV measurements is adequate6. 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 significant frequency peaks.

While it is not immediately evident, the light curve is impacted by both eclipses. The eclipse at primary superior conjunction (ϕ = 0.28, secondary isin front) is responsible for the rapid decrease of flux after periastron and to the slight discontinuity observed at ϕ = 0.28. This eclipse is wider than the eclipse at ϕ = −0.026 as a result of the longer time spent at apastron, and it is weaker primarily because of 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 Ṁ1 = −4.29 ± 0.05 (M⊙ yr−1) and log Ṁ2 = −4.42 ± 0.05 (M⊙ yr−1), agree within1σ with those derived from spectroscopy. The WWC excess reaches a maximum of ≈ ± 10%, which is roughly consistent with the flux changes in recombination lines (Fig. 12, Sect. 3.4). The exponent of γWWC = 1.5 ± 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 ± 5° when fixing all other parameters. Our fit implies an inclination of i = 60.4 ± 1.5°. This, in turns, implies dynamical masses of 74 ± 4 M⊙ and 69 ± 4 M⊙ 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.

thumbnail Fig. 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) plus WWC excess emission (green dotted line), which are multiplied by 0.5 for clarity. Lower panel: residuals between observation and model (O–C).
thumbnail Fig. 15

Same as Fig. 14, but showing the impact of changing the inclination by ± 5° while fixing all other parameters (green and blue lines, see legend).

3.6 Polarimetric analysis

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 owing 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 wind7. 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, T0, e, ω, and i, but also on the longitude of the ascending node Ω. Moreover, the model depends on the density and structure of the scattering medium, which is typically parameterized with two free fitting parameters: an effective optical depth τ*, and an exponent related to the structure of the wind, γ. For brevity, a full account of the relevant equations is available in Shenar et al. (2017b).

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 1), leaving only Ω, τ*, γ, and two arbitrary offset constants Q0 and U0 as fitting parameters. The fitting procedure, performed using the Python lmfit package, reveals that γ and τ* cannot be constrained well with the available data, while Ω can be. Therefore, we provide our preliminary value of Ω in Table 1. In Fig. 16, we show the comparison between the Q∕U data and our best-fitting models for two typical values of γ = −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

As discussed in Sect. 3.4, with shock velocities of the order of ≈ 1000 km 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 often, although not always, observed in binaries with two massive components (e.g. Pollock 1987; Corcoran et al. 1996; Nazé et al. 2007; Guerrero & Chu 2008). Similar to the optical emission excess in recombination lines, we 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 1. Tehrani (2019) reported R 144 to be a fairly hard X-ray source (kT = 4.2 ± 0.4 keV) and reported with an X-ray luminosity of log LX = 33.4 erg 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 portraying a comparable X-ray hardness (kT = 3.16 ± 0.03 keV), exceeds thisX-ray luminosity by almost two orders of magnitudes at peak luminosity (log LX = 35.3 erg s−1, Pollock et al. 2018). Moreover, while the X-ray flux of Mk 34 increases shortly before periastron (followed by a rapid decrease), this flux 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 in this work. 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, the winds likely do not reach their terminal value at the collision region. Assuming β = 1 for the velocity law, the wind velocities are roughly 75% (i.e. ≈ 1000 km s−1) at collision during periastron, and even lower for larger β 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.

thumbnail Fig. 16

Comparison between the observed and modelled Stokes parameters Q (blue line) and U (green line) for two distinct exponents γ. For γ = 1 (upper panel), we obtain Ω = 120 ± 11°, τ* = 0.22, Q0 = − 0.24 ± 0.02, U0 = 0.16 ± 0.02. For γ = −1 (lower panel), we obtain Ω = 125 ± 11°, τ* = 0.10 ± 0.02, Q0 = − 0.23 ± 0.02, U0 = 0.15 ± 0.02. The dashed lines bounding the curves show the impact of changing the inclination by ± 10°. See text for details.
thumbnail Fig. 17

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

4 Discussion

4.1 Evolutionary status

We estimate the evolutionary masses and ages of both components via the Bayesian statistics tool BONNSAI8 (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⊙ and over a wide range of initial rotation velocities. As input parameters, we use the derived values of T*, log L, XH, and upper limits on veq. 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.

Our final results are shown in Table 1. 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 significant eccentricity. We derive current masses of Mev, cur = 110 ± 9 and 100 ± 10 M⊙ and initial masses of Mini = 130 ± 13 and 119 ± 12 M⊙ 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 veq, ini ≈ 350–400 km s−1, which tend to chemically homogenise the stars. This in turn results in current rotation velocities of veq, cur ≈ 250 km s−1, which are substantially larger than the conservative upper limits derived in this work (veq ≲ 100 km 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. The current mass-loss rates retrieved by the BONNSAI tool are of the order of log Ṁ = −4.8 (M⊙ yr−1), which are 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 to invoke very high rotation to homogenise the stars. 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 1, 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. As a consequence of the larger mass-loss rates, the initial masses are larger than those inferred by BONNSAI, but the ages and current 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 1.

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 of the separation to diameter of the stars is a1 (1−e) R1−1 ≈ a2 (1−e) R2−1 ≈ 4$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.

thumbnail Fig. 18

HRD positions of the primary and secondary components of R 144, along with evolution tracks calculated by Köhler et al. (2015) for Mini = 100, 125, and 150 M⊙ with an initial rotation of veq, ini = 200 km s−1 and 350 km s−1. The colour corresponds to the surface hydrogen mass fraction.
thumbnail Fig. 19

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

4.1.2 Mass discrepancy

The derived evolutionary masses (Mev, cur ≈ 110 M⊙, 100 M⊙) disagree with our dynamical masses (Mdyn = 74 ± 4 M⊙, 69 ± 4 M⊙), 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; the Eddington factors are Γe = 0.78 ± 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 in this case. An inclination of 50° – a mere 10° decrease – would yield masses that are comparable to the evolutionary masses. However, the formal error on the inclination is only 1.5°, 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° difference. Moreover, in the framework of our light-curve model, reducing the inclination down to 50° also requires very large mass-loss rates of the order of log Ṁ = −3.9 (M⊙ 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 ≈140 + 130 M⊙ (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, these two stars exhibit significantly lower mass-loss rates by a factor of 3 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 comprises two classical WR stars, the components of R 144 may well be highly inflated because they appear much larger than what we would expect (R* ≈ 25 R⊙ vs. R* ≈ 1 R⊙); this was 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 (≈ 150 M⊙ and 130 M⊙), 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. The final orbital period Pf between the two would only grow with time due to wind mass loss as Pf/Pcur∝(Mtot,cur/Mtot,f)2$P_{\textrm{f}} / P_{\textrm{cur}} \propto \left(M_{\textrm{tot, cur}}/M_{\textrm{tot, f}}\right)^2$, where Mtot, cur and Mtot, f are the total current and final mass of the system. With an initial mass in excess of 140 M⊙, the primarycould be massive enough to undergo a pair-instability supernova, leaving no BH behind (e.g. Langer et al. 2007). Assuming that the fate of both components is to end their lives as ≈ 20–30 M⊙ BHs (Belczynski et al. 2010), we find that Pf ≈ 1–5 yr (depending on the current masses). That is, the system would become a BH plus 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 MV ≈ −8 mag, detecting any companions with current masses ≲ 50 M⊙ (corresponding to light contributions smaller than ≈10%) via spectroscopy would be difficult, and it is therefore very possible that companions with masses as large as ≈ 50 M⊙ 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, thereby forming a star plus BH binary. Such objects could be the progenitors of high-mass X-ray binaries.

thumbnail Fig. 20

2MASS JHK image (Skrutskie et al. 2006), highlighting the relative separation of R 144 from R 136 of about 4.1′ (≈ 60 pc). The median-subtracted PM arrow from Gaia eDR3 is also denoted.

4.3 Upper mass limit

Regardless of the evolutionary status of the system, our results imply initial masses of the order of ≈ 130–150 M⊙ for the system (see Sect. 4.1). The immediate conclusion is that stars as massive as 150 M⊙ form in the Local Group beyond reasonable doubt. A similar conclusion was reached by Tehrani et al. (2019), who reported initial masses of ≈ 150 M⊙ for the system Mk 34 – the younger twin of R 144. The 150 M⊙ limit is well within the 200–300 M⊙ 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⊙ 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 efforts 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 ≲ 50 days) binaries (e.g., Schnurr et al. 2009).

4.4 Runaway status

R 144 is situated in relative isolation with respect to the neighbouring very massive and young cluster R 136, roughly 4.1′ apart (i.e. ≈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 V0 = 210 ± 20 km s−1 measured in our study puts it roughly 60 km s−1 below the bulk of the Tarantula clusters NGC 2070 and R 136 (cf. Almeida et al. 2017) of ≈ 270 ± 10 km s−1. While the measurement of absolute velocities of WR stars is model-dependent (Sect. 3.1), the observed offset is larger than 3σ when considering our conservative error estimate of 20 km 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 2021). We computed the mean, median, and standard deviations of the PM components for all objects with a PM accuracy better than 0.1 mas yr−1 in both PM components (464 objects in total). The standard deviation is σ→PM,4′ = (0.16,0.17)$\overrightarrow{\sigma}_{\textrm{PM, 4'}}\,{=}\,(0.16, 0.17)$mas 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 ΔPM→R144 = (0.12 ± 0.02,0.28 ± 0.03)$\overrightarrow{ \Delta {\textrm{PM}}}_{\textrm{R144}}\,{=}\,(0.12 \,{\pm}\, 0.02, 0.28\,{\pm}\,0.03)$mas 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 higher than 3σ and has a transverse velocity component of ≈ 70 km s−1. If R 144 wasejected from R 136 in the past ≈ 2 Myr, a distance of 4.1′ would require a runaway velocity v ≳ 30 km s−1. For a transverse velocity component of ≈ 70 km s−1, the ejection needs to have happened roughly 1 Myr ago.

Assuming R 144 was 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 neighbouring 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 WR 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 PoWR code to derive the physical parameters of both components. Below, we summarize out conclusions:

  • R 144 is a binary comprising two WR stars (WN5/6h + WN6/7h) of nearly equal masses (M1 sin3 i = 48.3 ± 1.8 and M2 sin3 i = 45.5 ± 1.8 M⊙) orbiting their centre of mass in an eccentric (e = 0.51) orbit with a period of P = 74.2 days. The primary is slightly hotter and more luminous than the secondary, and both stars exhibit substantial residual hydrogen (XH ≈ 0.4) and nitrogen enrichment. Assuming the stars are still on the main sequence, their inferred age is ≈2 Myr;

  • The systemic velocity, PM 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 WWCs and wind eclipses. The model enables us to derive the orbital inclination with a small formal error, i = 60.4 ± 1.5°. In turn, this results in accurate dynamical masses of 74 ± 4 M⊙ and 69 ± 4 M⊙, makingR 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⊙] for the primary and secondary, which result in evolutionary masses of the order of 100–110 M⊙. A good agreement between evolutionary and dynamical masses would be achieved for i = 50°, but this inclination does not seem to be consistent with the light curve; however, see the discussion in Appendix D;

  • Taken at face value, the derived dynamical masses and luminosities imply that both components have classical Eddington factors of Γe = 0.78 ± 0.10 and are thus expected to be inflated as a result of 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.

Acknowledgements

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 Ochsenbein et al. (2000). The research leading to these results has received 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 programme 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 programme BinWaves with project number 639.042.728.

Appendix A Observation log and RV measurements

Table A.1

Observation log and measured RVs.

Appendix B TESS data

Figure B.1 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, although a comparison with the Gaia early DR3 catalogue (Gaia Collaboration 2018, 2021) suggests the contamination does not exceed 1–2%. The high-precision TESS photometry data reveal interesting substructures in the light curve that are indicative of pulsational activity, although the time coverage does not enable a reliable period analysis of these structure to identify their possible origin.

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

thumbnail Fig. B.1

Short-cadence TESS light curve (black dots), obtained for sectors 27–30. Flagged (less reliable) data points are denoted in light grey, and the dashed grey vertical lines indicate the time stamps of the TESS angular momentum dumps.

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 these data to refine the orbital period of the system, as described in Sect. 3.5. Nevertheless, we present in Fig. B.2 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 reproduces the light curve well. The TESS data seem to suggest a deeper secondary minimum at ϕ = 0.28, which implies either even higher inclinations or an increased mass-loss rate for the secondary. But because this is the case only for 1.5 cycles, we avoid an over-interpretation of this mismatch.

thumbnail Fig. B.2

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

Appendix C PHOEBE modelling of the light curve

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 RV curves and spectroscopic fitting, we endeavored to derive a self-consistent model. However, our efforts led to the model seen in Fig. C.1. 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 ϕ = −0.026 (inferior conjunction), suggesting that tidal distortion is not the dominant mechanism governing the shape of the light curve.

thumbnail Fig. C.1

Our PHOEBE model of R 144 plot, using the parameters given in Table 1 as input. Owing to the differentscale of amplitude of the data and model, the data are not plotted for clarity (cf. Fig. 13).

Appendix D Systematic uncertainties in the light-curve model

Modelling of the ASAS-SN data enables an accurate derivation of the orbital inclination of i = 60.4 ± 1.5°. 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 D.1. Below we summarize 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.

    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 (χ2) remains unchanged, which is anticipated as Ṁ1, Ṁ2 and f1∕f2 are degenerate.

  • 2.

    Altering M1.,M2.$\dot{M_1}, \dot{M_2}$: Similarly, we fixed Ṁ to the values derived spectroscopically. This only leads to a slight increase (by 1°) of the orbital inclination to compensate for the net decrease of density in the wind of the primary. The fit quality is comparable.

  • 3.

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

  • 4.

    Varying α, R*, or v∞: We altered the number of electrons per baryon α, stellar radii R*, and terminal velocities v∞ within formal errors. Some impact, although a small one (≈0.05 dex) could be seen in the mass-loss rates, but no notable change was obtained in the value of i.

  • 5.

    Ionization structure in the wind: The model relies on a constant α, 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 α. However, as α is primarily a multiplicative parameter determining the strength of absorption, any changes in α at first order impact the derived mass-loss rates from this method; the impact on the inclination is second order.

  • 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.

Table D.1

Inferred parameters from our hybrid light-curve model (2829 fitting points, combined RV plus light curve).

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, and finite sizes of stellar discs). 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 β = 2

We providean extension of the wind-eclipse model of Lamontagne et al. (1996) for the case of a β-law velocity exponent of β = 2. To compute the optical depth from the wind we need to evaluate the integral τ = k∫−ϵ∞d(z/D)(r/)2(1−R*/r)β ,\begin{eqnarray*} \tau\,{=}\,k\int_{-\epsilon}^{\infty}\frac{d(z/\mathcal{D})}{(r/\mathcal{})^2(1-R_*/r)^{\beta}}, \end{eqnarray*}(E.1)

where ϵ = sinicosϖ, (rD)2 = cos2icos2ϖ+sin2ϖ+(zD)2,\begin{eqnarray*} \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{eqnarray*}(E.2)

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

Definingγ = cos2icos2ϖ+sin2ϖ$\gamma\,{=}\,\sqrt{\cos^2 i \cos^2\varpi + \sin^2 \varpi}$, ψ = γ2−δ2$\psi\,{=}\,\sqrt{\gamma^2-\delta^2}$ and δ = R*∕a, the analytical solution to the integral for the case β = 2 is given by τk = δψ2(1+ϵϵ2+γ2+δ2ϵ2ψ2+ϵ2)+γ2ψ3(arcsin(δγ)+|ϵ|ϵarcsin(δ2ϵ2γ2(ψ2+ϵ2)))+πγ22ψ3(π2+arctan(ϵψ)). \begin{eqnarray*} \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)\nonumber \\ &&+\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).\nonumber \end{eqnarray*}(E.3)

References

  1. Abbott, R., Abbott, T. D., Abraham, S., et al. 2020, Phys. Rev. Lett., 125, 101102 [CrossRef] [Google Scholar]
  2. Almeida, L. A., Sana, H., Taylor, W., et al. 2017, A&A, 598, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Barbá, R. H., Gamen, R., Arias, J. I., et al. 2010, Rev. Mex. Astron. Astrofis. Conf. Ser., 38, 30 [Google Scholar]
  4. Bastian, N., & Lardo, C. 2018, ARA&A, 56, 83 [Google Scholar]
  5. Belczynski, K., Bulik, T., Fryer, C. L., et al. 2010, ApJ, 714, 1217 [Google Scholar]
  6. Bestenlehner, J. M. 2020, MNRAS, 493, 3938 [Google Scholar]
  7. Bestenlehner, J. M., Gräfener, G., Vink, J. S., et al. 2014, A&A, 570, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Bestenlehner, J. M., Crowther, P. A., Caballero-Nieves, S. M., et al. 2020, MNRAS, 499, 1918 [Google Scholar]
  9. Bonanos, A. Z., Stanek, K. Z., Udalski, A., et al. 2004, ApJ, 611, L33 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bonanos, A. Z., Massa, D. L., Sewilo, M., et al. 2009, AJ, 138, 1003 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bowman, D. M., Burssens, S., Pedersen, M. G., et al. 2019, Nat. Astron., 3, 760 [NASA ADS] [CrossRef] [Google Scholar]
  12. Breysacher, J., Azzopardi, M., & Testor, G. 1999, A&AS, 137, 117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Broos, P. S., Townsley, L. K., Feigelson, E. D., et al. 2010, ApJ, 714, 1582 [NASA ADS] [CrossRef] [Google Scholar]
  14. Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Brown, J. C., McLean, I. S., & Emslie, A. G. 1978, A&A, 68, 415 [NASA ADS] [Google Scholar]
  16. Brown, J. C., Aspin, C., Simmons, J. F. L., & McLean, I. S. 1982, MNRAS, 198, 787 [NASA ADS] [CrossRef] [Google Scholar]
  17. Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245 [NASA ADS] [CrossRef] [Google Scholar]
  18. Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [NASA ADS] [CrossRef] [Google Scholar]
  19. 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 [Google Scholar]
  20. Cherepashchuk, A. M. 1976, Sov. Astron. Lett., 2, 138 [NASA ADS] [Google Scholar]
  21. Cleveland, W. S. 1979, J. Am. Stat. Assoc., 74, 829 [Google Scholar]
  22. Corcoran, M. F., Stevens, I. R., Pollock, A. M. T., et al. 1996, ApJ, 464, 434 [NASA ADS] [CrossRef] [Google Scholar]
  23. Crowther, P. A., & Dessart, L. 1998, MNRAS, 296, 622 [NASA ADS] [CrossRef] [Google Scholar]
  24. Crowther, P. A., & Smith, L. J. 1997, A&A, 320, 500 [NASA ADS] [Google Scholar]
  25. Crowther, P. A., & Walborn, N. R. 2011, MNRAS, 416, 1311 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  26. Crowther, P. A., Schnurr, O., Hirschi, R., et al. 2010, MNRAS, 408, 731 [NASA ADS] [CrossRef] [Google Scholar]
  27. Cutri, R. M., Wright, E. L., Conrow, T., et al. 2014, VizieR Online Data Catalog: II/328 [Google Scholar]
  28. David-Uraz, A., Moffat, A. F. J., Chené, A.-N., et al. 2012, MNRAS, 426, 1720 [NASA ADS] [CrossRef] [Google Scholar]
  29. de Koter, A., Heap, S. R., & Hubeny, I. 1997, ApJ, 477, 792 [Google Scholar]
  30. Demers, H., Moffat, A. F. J., Marchenko, S. V., Gayley, K. G., & Morel, T. 2002, ApJ, 577, 409 [NASA ADS] [CrossRef] [Google Scholar]
  31. Doran, E. I., Crowther, P. A., de Koter, A., et al. 2013, A&A, 558, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Dsilva, K., Shenar, T., Sana, H., & Marchant, P. 2020, A&A, 641, A26 [CrossRef] [EDP Sciences] [Google Scholar]
  33. Evans, C. J., Taylor, W. D., Hénault-Brunet, V., et al. 2011, A&A, 530, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Figer, D. F. 2005, Nature, 434, 192 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  35. Figer, D. F., Najarro, F., Gilmore, D., et al. 2002, ApJ, 581, 258 [NASA ADS] [CrossRef] [Google Scholar]
  36. Fryer, C. L., Woosley, S. E., & Heger, A. 2001, ApJ, 550, 372 [NASA ADS] [CrossRef] [Google Scholar]
  37. Fujii, M. S., & Portegies Zwart, S. 2011, Science, 334, 1380 [NASA ADS] [CrossRef] [Google Scholar]
  38. Fullerton, A. W., Massa, D. L., & Prinja, R. K. 2006, ApJ, 637, 1025 [NASA ADS] [CrossRef] [Google Scholar]
  39. Gaia Collaboration 2018, VizieR Online Data Catalog, I/345 [Google Scholar]
  40. Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Gieles, M., Charbonnel, C., Krause, M. G. H., et al. 2018, MNRAS, 478, 2461 [NASA ADS] [CrossRef] [Google Scholar]
  42. Gilkis, A., Shenar, T., Ramachandran, V., et al. 2021, MNRAS, 503, 1884 [CrossRef] [Google Scholar]
  43. Giménez-García, A., Shenar, T., Torrejón, J. M., et al. 2016, A&A, 591, A26 [EDP Sciences] [Google Scholar]
  44. González, J. F., & Levato, H. 2006, A&A, 448, 283 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Gräfener, G. 2021, A&A, 647, A13 [CrossRef] [EDP Sciences] [Google Scholar]
  46. Gräfener, G., & Hamann, W. R. 2008, A&A, 482, 945 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Gräfener, G., Koesterke, L., & Hamann, W. R. 2002, A&A, 387, 244 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Gräfener, G., Vink, J. S., de Koter, A., & Langer, N. 2011, A&A, 535, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Gräfener, G., Owocki, S. P., & Vink, J. S. 2012, A&A, 538, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Grassitelli, L., Chené, A. N., Sanyal, D., et al. 2016, A&A, 590, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Guerrero, M. A., & Chu, Y.-H. 2008, ApJS, 177, 216 [NASA ADS] [CrossRef] [Google Scholar]
  52. Hadrava, P. 1995, A&AS, 114, 393 [NASA ADS] [Google Scholar]
  53. Hainich, R., Rühling, U., Todt, H., et al. 2014, A&A, 565, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Hamann, W. R., & Gräfener, G. 2003, A&A, 410, 993 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Hamann, W. R., & Gräfener, G. 2004, A&A, 427, 697 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Hawcroft, C., Sana, H., Mahy, L., et al. 2021, A&A, submitted [Google Scholar]
  57. Higgins, E. R., & Vink, J. S. 2019, A&A, 622, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Hillier, D. J. 1991, A&A, 247, 455 [NASA ADS] [Google Scholar]
  59. Howarth, I. D. 1983, MNRAS, 203, 301 [NASA ADS] [CrossRef] [Google Scholar]
  60. Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, in SPIE Conf. Ser., 9913, Software and Cyberinfrastructure for Astronomy IV, ed. G. Chiozzi & J. C. Guzman, 99133E [Google Scholar]
  61. Jermyn, A. S., Tout, C. A., & Chitre, S. M. 2018, MNRAS, 480, 5427 [NASA ADS] [CrossRef] [Google Scholar]
  62. Kaufer, A., Wolf, B., Andersen, J., & Pasquini, L. 1997, The Messenger, 89, 1 [NASA ADS] [Google Scholar]
  63. Kochanek, C. S., Shappee, B. J., Stanek, K. Z., et al. 2017, PASP, 129, 104502 [Google Scholar]
  64. Köhler, K., Langer, N., de Koter, A., et al. 2015, A&A, 573, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Kozai, Y. 1962, AJ, 67, 591 [Google Scholar]
  66. Lamontagne, R., Moffat, A. F. J., Drissen, L., Robert, C., & Matthews, J. M. 1996, AJ, 112, 2227 [NASA ADS] [CrossRef] [Google Scholar]
  67. Langer, N., Norman, C. A., de Koter, A., et al. 2007, A&A, 475, L19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Larson, R. B., & Starrfield, S. 1971, A&A, 13, 190 [NASA ADS] [Google Scholar]
  69. Lépine, S., & Moffat, A. F. J. 1999, ApJ, 514, 909 [NASA ADS] [CrossRef] [Google Scholar]
  70. Lidov, M. L. 1962, Planet. Space Sci., 9, 719 [Google Scholar]
  71. Lightkurve Collaboration (Cardoso, J. V. d. M., et al.) 2018, Lightkurve: Kepler and TESS time series analysis in Python [Google Scholar]
  72. Lohr, M. E., Clark, J. S., Najarro, F., et al. 2018, A&A, 617, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Magalhaes, A. M., Benedetti, E., & Roland, E. H. 1984, PASP, 96, 383 [NASA ADS] [CrossRef] [Google Scholar]
  74. Mahy, L., Almeida, L. A., Sana, H., et al. 2020a, A&A, 634, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Mahy, L., Sana, H., Abdul-Masih, M., et al. 2020b, A&A, 634, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Maíz Apellániz, J., Sana, H., Barbá, R. H., Le Bouquin, J. B., & Gamen, R. C. 2017, MNRAS, 464, 3561 [NASA ADS] [CrossRef] [Google Scholar]
  77. Maíz Apellániz, J., Trigueros Páez, E., Jiménez Martínez, I., et al. 2019, in Highlights on Spanish Astrophysics X, eds. B. Montesinos, A. Asensio Ramos, F. Buitrago, R. Schödel, E. Villaver, S. Pérez-Hoyos, & I. Ordóñez-Etxeberria, 420 [Google Scholar]
  78. Marchant, P., Renzo, M., Farmer, R., et al. 2019, ApJ, 882, 36 [Google Scholar]
  79. Marchenko, S. V., Moffat, A. F. J., & Eenens, P. R. J. 1998, PASP, 110, 1416 [NASA ADS] [CrossRef] [Google Scholar]
  80. Marchenko, S. V., Moffat, A. F. J., Ballereau, D., et al. 2003, ApJ, 596, 1295 [NASA ADS] [CrossRef] [Google Scholar]
  81. Martins, F., Hillier, D. J., Paumard, T., et al. 2008, A&A, 478, 219 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  82. Moffat, A. F. J. 1989, ApJ, 347, 373 [NASA ADS] [CrossRef] [Google Scholar]
  83. Moffat, A. F. J., Drissen, L., Lamontagne, R., & Robert, C. 1988, ApJ, 334, 1038 [NASA ADS] [CrossRef] [Google Scholar]
  84. Moffat, A. F. J., Marchenko, S. V., Seggewiss, W., et al. 1998, A&A, 331, 949 [NASA ADS] [Google Scholar]
  85. Monet, D., Canzian, B., Harris, H., et al. 1998, VizieR Online Data Catalog: I/243 [Google Scholar]
  86. Najarro, F., Figer, D. F., Hillier, D. J., & Kudritzki, R. P. 2004, ApJ, 611, L105 [NASA ADS] [CrossRef] [Google Scholar]
  87. Nazé, Y., Corcoran, M. F., Koenigsberger, G., & Moffat, A. F. J. 2007, ApJ, 658, L25 [NASA ADS] [CrossRef] [Google Scholar]
  88. Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2014, LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python [Google Scholar]
  89. Ochsenbein, F., Bauer, P., & Marcout, J. 2000, A&AS, 143, 23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Oey, M. S., & Clarke, C. J. 2005, ApJ, 620, L43 [NASA ADS] [CrossRef] [Google Scholar]
  91. Oskinova, L. M., Hamann, W. R., & Feldmeier, A. 2007, A&A, 476, 1331 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Parker, J. W. 1992, PASP, 104, 1107 [NASA ADS] [CrossRef] [Google Scholar]
  93. Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
  94. Pedersen, M. G., Aerts, C., Pápics, P. I., & Rogers, T. M. 2018, A&A, 614, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Petrovic, J., Pols, O., & Langer, N. 2006, A&A, 450, 219 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Pojmanski, G. 2002, Acta Astron., 52, 397 [NASA ADS] [Google Scholar]
  97. Pollock, A. M. T. 1987, ApJ, 320, 283 [NASA ADS] [CrossRef] [Google Scholar]
  98. Pollock, A. M. T., Crowther, P. A., Tehrani, K., Broos, P. S., & Townsley, L. K. 2018, MNRAS, 474, 3228 [NASA ADS] [CrossRef] [Google Scholar]
  99. Prša, A., Conroy, K. E., Horvat, M., et al. 2016, ApJS, 227, 29 [NASA ADS] [CrossRef] [Google Scholar]
  100. Puls, J., Markova, N., Scuderi, S., et al. 2006, A&A, 454, 625 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  101. Ramachandran, V., Hamann, W. R., Oskinova, L. M., et al. 2019, A&A, 625, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  102. Ramiaramanantsoa, T., Ignace, R., Moffat, A. F. J., et al. 2019, MNRAS, 490, 5921 [CrossRef] [Google Scholar]
  103. Rauw, G., Vreux, J. M., & Bohannan, B. 1999, ApJ, 517, 416 [NASA ADS] [CrossRef] [Google Scholar]
  104. Rauw, G., De Becker, M., Nazé, Y., et al. 2004, A&A, 420, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  105. Renzo, M., Zapartas, E., de Mink, S. E., et al. 2019, A&A, 624, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. Richardson, N. D., Shenar, T., Roy-Loubier, O., et al. 2016, MNRAS, 461, 4115 [CrossRef] [Google Scholar]
  107. Richardson, N. D., Russell, C. M. P., St-Jean, L., et al. 2017, MNRAS, 471, 2715 [NASA ADS] [CrossRef] [Google Scholar]
  108. Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, J. Astron. Telescopes Instrum. Syst., 1, 014003 [Google Scholar]
  109. Sana, H., de Koter, A., de Mink, S. E., et al. 2013a, A&A, 550, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  110. Sana, H., van Boeckel, T., Tramper, F., et al. 2013b, MNRAS, 432, L26 [NASA ADS] [CrossRef] [Google Scholar]
  111. Sander, A., Shenar, T., Hainich, R., et al. 2015, A&A, 577, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  112. Sanyal, D., Langer, N., Szécsi, D.,-C Yoon, S., & Grassitelli, L. 2017, A&A, 597, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  113. Schmutz, W., Hamann, W. R., & Wessolowski, U. 1989, A&A, 210, 236 [Google Scholar]
  114. Schneider, F. R. N., Langer, N., de Koter, A., et al. 2014, A&A, 570, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  115. Schneider, F. R. N., Ramírez-Agudelo, O. H., Tramper, F., et al. 2018a, A&A, 618, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  116. Schneider, F. R. N., Sana, H., Evans, C. J., et al. 2018b, Science, 359, 69 [Google Scholar]
  117. Schnurr, O., Casoli, J., Chené, A. N., Moffat, A. F. J., & St-Louis, N. 2008a, MNRAS, 389, L38 [NASA ADS] [CrossRef] [Google Scholar]
  118. Schnurr, O., Moffat, A. F. J., St-Louis, N., Morrell, N. I., & Guerrero, M. A. 2008b, MNRAS, 389, 806 [NASA ADS] [CrossRef] [Google Scholar]
  119. Schnurr, O., Chené, A. N., Casoli, J., Moffat, A. F. J., & St-Louis, N. 2009, MNRAS, 397, 2049 [NASA ADS] [CrossRef] [Google Scholar]
  120. Schootemeijer, A., & Langer, N. 2018, A&A, 611, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  121. Schootemeijer, A., Langer, N., Grin, N. J., & Wang, C. 2019, A&A, 625, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  122. Schwarzenberg-Czerny, A. 1997, ApJ, 489, 941 [NASA ADS] [CrossRef] [Google Scholar]
  123. Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, ApJ, 788, 48 [NASA ADS] [CrossRef] [Google Scholar]
  124. Shenar, T., Hamann, W. R., & Todt, H. 2014, A&A, 562, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  125. Shenar, T., Oskinova, L., Hamann, W. R., et al. 2015, ApJ, 809, 135 [Google Scholar]
  126. Shenar, T., Oskinova, L. M., Järvinen, S. P., et al. 2017a, A&A, 606, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  127. Shenar, T., Richardson, N. D., Sablowski, D. P., et al. 2017b, A&A, 598, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  128. Shenar, T., Hainich, R., Todt, H., et al. 2018, A&A, 616, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  129. Shenar, T., Sablowski, D. P., Hainich, R., et al. 2019, A&A, 627, A151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  130. Shenar, T., Bodensteiner, J., Abdul-Masih, M., et al. 2020, A&A, 639, A6 [CrossRef] [EDP Sciences] [Google Scholar]
  131. Simmons, J. F. L., & Boyle, C. B. 1984, A&A, 134, 368 [NASA ADS] [Google Scholar]
  132. Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163 [NASA ADS] [CrossRef] [Google Scholar]
  133. Smith, L. F., Shara, M. M., & Moffat, A. F. J. 1996, MNRAS, 281, 163 [NASA ADS] [CrossRef] [Google Scholar]
  134. St.-Louis, N., Moffat, A. F. J., Drissen, L., Bastien, P., & Robert, C. 1988, ApJ, 330, 286 [NASA ADS] [CrossRef] [Google Scholar]
  135. St-Louis, N., Moffat, A. F. J., Marchenko, S., & Pittard, J. M. 2005, ApJ, 628, 953 [NASA ADS] [CrossRef] [Google Scholar]
  136. Stellingwerf, R. F. 1978, ApJ, 224, 953 [Google Scholar]
  137. Stevens, I. R., Blondin, J. M., & Pollock, A. M. T. 1992, ApJ, 386, 265 [NASA ADS] [CrossRef] [Google Scholar]
  138. Sundqvist, J. O., Puls, J., & Feldmeier, A. 2010, A&A, 510, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  139. Šurlan, B., Hamann, W. R., Aret, A., et al. 2013, A&A, 559, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  140. Tehrani, K. 2019, PhD thesis, University of Sheffield, UK [Google Scholar]
  141. Tehrani, K. A., Crowther, P. A., Bestenlehner, J. M., et al. 2019, MNRAS, 484, 2692 [NASA ADS] [CrossRef] [Google Scholar]
  142. Thompson, S. E., Everett, M., Mullally, F., et al. 2012, ApJ, 753, 86 [Google Scholar]
  143. Todt, H., Sander, A., Hainich, R., et al. 2015, A&A, 579, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  144. Toonen, S., Portegies Zwart, S., Hamers, A. S., & Bandopadhyay, D. 2020, A&A, 640, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  145. Townsley, L. K., Broos, P. S., Corcoran, M. F., et al. 2011, ApJS, 194, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  146. Tramper, F., Sana, H., Fitzsimons, N. E., et al. 2016, MNRAS, 455, 1275 [NASA ADS] [CrossRef] [Google Scholar]
  147. Usov, V. V. 1992, ApJ, 389, 635 [NASA ADS] [CrossRef] [Google Scholar]
  148. Vigna-Gómez, A., Toonen, S., Ramirez-Ruiz, E., et al. 2021, ApJ, 907, L19 [Google Scholar]
  149. Vink, J. S. 2018, A&A, 615, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  150. Vink, J. S., & Harries, T. J. 2017, A&A, 603, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  151. Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  152. Vink, J. S., Muijres, L. E., Anthonisse, B., et al. 2011, A&A, 531, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  153. Williams, P. M., Dougherty, S. M., Davis, R. J., et al. 1997, MNRAS, 289, 10 [NASA ADS] [Google Scholar]
  154. Willis, A. J., Crowther, P. A., Fullerton, A. W., et al. 2004, ApJS, 154, 651 [NASA ADS] [CrossRef] [Google Scholar]
  155. Woosley, S. E., Blinnikov, S., & Heger, A. 2007, Nature, 450, 390 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  156. Zahn, J. P. 1977, A&A, 500, 121 [NASA ADS] [Google Scholar]
  157. Zucker, S., & Mazeh, T. 1994, ApJ, 420, 806 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
4

τ = 20 is a standard convention; any value above ≈10 would be appropriate owing to the exponential increase in τ(r).

5

Sometimes also denoted as λ.

6

The 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.

7

The equations were later corrected by Simmons & Boyle (1984).

8

The BONNSAI web-service is available at https://www.astro.uni-bonn.de/stars/bonnsai

All Tables

Table 1

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

In the text
Table A.1

Observation log and measured RVs.

In the text
Table D.1

Inferred parameters from our hybrid light-curve model (2829 fitting points, combined RV plus light curve).

In the text

All Figures

thumbnail Fig. 1

X-shooter spectrum of R 144 taken on 20 February 2014 (MJD = 56708.04). An overview of the optical part of the spectrum is shown.
In the text
thumbnail Fig. 2

Three X-shooter spectra (MJD = 55452.89, 55585.65, 56708.04, see legend) corresponding to phases close to the RV extremes (ϕ = 0.05, 0.84 with ephemeris derived in this study) and the systemic velocity (ϕ = 0.96) of the N IV λ4058 and N V λ4945 lines. Both binary components are seen in the N IVλ4058 line, but only the hotter primary is seen in the N Vλ4945 line. Conversion to Doppler space is with reference to rest wavelengths.
In the text
thumbnail Fig. 3

As Fig. 2, but showing three FLAMES-UVES spectra (MJD = 56645.57, 56645.57, 56697.69), and focussing on the N V λλ4604, 4620 and N III λλ4634, 4641 lines. Both components contribute to the N V λλ4604, 4620 lines, while the cooler secondary dominates the N IIIλλ4634, 4641 lines.
In the text
thumbnail Fig. 4

Dynamical spectra for the N Vλ4945, N III λλ4634, 4641, and N V λλ4604, 4620 lines. The red and green curves depict the orbital solution derived in Sect. 3.1 for the primary and secondary components, respectively.
In the text
thumbnail Fig. 5

Radial-velocity measurements of the N Vλ4945 line (primary) and N IIIλλ4634, 4641 triplet (secondary),along with the RV curves corresponding to the orbital parameters given in Table 1.
In the text
thumbnail Fig. 6

Upper panels: comparison between the disentangled spectra obtained for the best-fitting K1 and K2 values, their sum, and observations at two phases close to RV extremes for the N IV λ4058. Lower panel:reduced χ2 as a functionof K1 and K2 for the N IV λ4058 line. The analysis is performed using the X-shooter data. The minimum is at K1 = 129 ± 6 km s−1, K2 = 134 ± 4 km s−1.
In the text
thumbnail Fig. 7

As Fig. 6, but for N V λλ4604, 4620. The analysis is performed using the FLAMES-UVES data and X-shooter data. The minimum is at K1 = 130 ± 3 km s−1, K2 = 134 ± 4 km s−1.
In the text
thumbnail Fig. 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).
In the text
thumbnail Fig. 9

Comparison between observed SED (blue squares and lines, upper panel) and normalized FUSE, IUE, and X-shooter (ϕ = 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 indicated in red.
In the text
thumbnail Fig. 10

Schematic of R 144 as seen during primary inferior conjunction (ϕ ≈ −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 WWC cone. The line-forming regions of several diagnostic lines for the primary star, as calculated from our tailored model atmosphere (Sect. 3.3), are denoted. The line-forming regions of the secondary star are comparable.
In the text
thumbnail Fig. 11

Comparison between an observed spectrum of the N Vλ4945 multiplet and model spectra for the primary calculated with the parameters in Table 1, 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).
In the text
thumbnail Fig. 12

Variations of EW with phase for He Iλ4472 (4465–4490 Å), He Iλ5876 (5875–5900 Å), He IIλ4686 (4670–4720 Å), Hβ (4840–4890 Å), N IIIλλ4634, 4641 (4620–4650 Å), and N V λλ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.
In the text
thumbnail Fig. 13

ASAS-SN light curve of R 144 phased with the ephemeris given in Table 1. The phases of inferior and superior conjunction (hotter primary in front and behind the cooler secondary, respectively) and periastron passage are indicated in the inset.
In the text
thumbnail Fig. 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) plus WWC excess emission (green dotted line), which are multiplied by 0.5 for clarity. Lower panel: residuals between observation and model (O–C).
In the text
thumbnail Fig. 15

Same as Fig. 14, but showing the impact of changing the inclination by ± 5° while fixing all other parameters (green and blue lines, see legend).
In the text
thumbnail Fig. 16

Comparison between the observed and modelled Stokes parameters Q (blue line) and U (green line) for two distinct exponents γ. For γ = 1 (upper panel), we obtain Ω = 120 ± 11°, τ* = 0.22, Q0 = − 0.24 ± 0.02, U0 = 0.16 ± 0.02. For γ = −1 (lower panel), we obtain Ω = 125 ± 11°, τ* = 0.10 ± 0.02, Q0 = − 0.23 ± 0.02, U0 = 0.15 ± 0.02. The dashed lines bounding the curves show the impact of changing the inclination by ± 10°. See text for details.
In the text
thumbnail Fig. 17

X-ray light curve of R 144, also shown aggregated into adaptive wider phase intervals. See text for details.
In the text
thumbnail Fig. 18

HRD positions of the primary and secondary components of R 144, along with evolution tracks calculated by Köhler et al. (2015) for Mini = 100, 125, and 150 M⊙ with an initial rotation of veq, ini = 200 km s−1 and 350 km s−1. The colour corresponds to the surface hydrogen mass fraction.
In the text
thumbnail Fig. 19

As Fig. 18, but showing MESA tracks with parameters corresponding to those given in Table 1. These tracks are calculated using the set-up described by Gräfener (2021), which includes boosted mass-loss rates.
In the text
thumbnail Fig. 20

2MASS JHK image (Skrutskie et al. 2006), highlighting the relative separation of R 144 from R 136 of about 4.1′ (≈ 60 pc). The median-subtracted PM arrow from Gaia eDR3 is also denoted.
In the text
thumbnail Fig. B.1

Short-cadence TESS light curve (black dots), obtained for sectors 27–30. Flagged (less reliable) data points are denoted in light grey, and the dashed grey vertical lines indicate the time stamps of the TESS angular momentum dumps.
In the text
thumbnail Fig. B.2

As Fig. 15, but comparing with the phased TESS data instead of the ASAS-SN data.
In the text
thumbnail Fig. C.1

Our PHOEBE model of R 144 plot, using the parameters given in Table 1 as input. Owing to the differentscale of amplitude of the data and model, the data are not plotted for clarity (cf. Fig. 13).
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.