key: cord-0326683-to0ihncz
authors: Bestenlehner, Joachim M.; Crowther, Paul A.; Broos, Patrick S.; Pollock, Andrew M. T.; Townsley, Leisa K.
title: Melnick 33Na: a very massive colliding wind binary system in 30 Doradus
date: 2021-11-30
journal: nan
DOI: 10.1093/mnras/stab3521
sha: c5c0f39e1e7d41a5d175164c7b9fb155384a8506
doc_id: 326683
cord_uid: to0ihncz

We present spectroscopic analysis of the luminous X-ray source Melnick 33Na (Mk 33Na, HSH95 16) in the LMC 30 Doradus region (Tarantula Nebula), utilising new time-series VLT/UVES spectroscopy. We confirm Mk 33Na as a double-lined O-type spectroscopic binary with a mass ratio $q = 0.63 pm 0.02$, $e = 0.33 pm 0.01$ and orbital period of $18.3 pm 0.1$ days, supporting the favoured period from X-ray observations obtained via the Tarantula -- Revealed by X-rays (T-ReX) survey. Disentangled spectra of each component provide spectral types of OC2.5 If* and O4 V for the primary and secondary respectively - unusually for an O supergiant the primary exhibits strong CIV 4658 emission and weak NV 4603-20, justifying the OC classification. Spectroscopic analysis favours extreme physical properties for the primary ($T_{rm eff} = 50$ kK, $log L/L_{odot} = 6.15$) with system components of $M_{1} = 83 pm 19 M_{odot}$ and $M_{2} = 48 pm 11 M_{odot}$ obtained from evolutionary models, which can be reconciled with results from our orbital analysis (e.g. $M_{1} sin^3 i = 20.0 pm 1.2 M_{odot}$) if the system inclination is $sim 38^{circ}$ and it has an age of 0.9 to 1.6 Myr. This establishes Mk 33Na as one of the highest mass binary systems in the LMC, alongside other X-ray luminous early-type binaries Mk34 (WN5h+WN5h), R144 (WN5/6h+WN6/7h) and especially R139 (O6.5,Iafc+O6,Iaf).

The mass of a star is its fundamental defining property, but direct measurements remain elusive with the exception of visual binaries or eclipsing spectroscopic binaries (Andersen 1991; Serenelli et al. 2021) . Consequently mass estimates usually follow from spectroscopic results or comparison with evolutionary predictions. Uncertainties of evolutionary models increase with stellar mass (e.g Martins & Palacios 2013) due to uncertainties in nuclear reaction rates, stellar structure, internal mixing processes and mass-loss properties (Langer 2012) , so stellar mass estimates of high mass stars remain poorly constrained. Eclipsing double-lined spectroscopic binaries (SB2s) are the ideal gold standard calibrators to test single-star evolutionary models, because they permit accurate determinations of the dynamical mass, stellar parameters and chemical abundances for both components. Most studies of such systems have focused on stars up to 15 which only consider convective core-overshooting and rotational mixing (e.g. Southworth et al. 2004; Tkachenko et al. 2014; Pavlovski et al. 2018) . Higgins & Vink (2019) also included the effect of mass loss, which strips the stellar envelope and reveals CNO-processed material. They calibrated their grid of stellar models with the 40 + 34 eclipsing massive supergiant binary HD 166734 (O7.5If + O9I(f)) using the results for this system from Mahy et al. (2017) .

Despite the high binary frequency amongst massive stars (Sana Way there are very few massive eclipsing binaries whose component masses exceed a few tens of solar masses, including NGC 3603 A1 with system components of 116+89 (Schnurr et al. 2008) , F2 in the Arches cluster with 82+60 (Lohr et al. 2018) and WR20a with 71+69 (Bonanos et al. 2004) . Such systems are also scarce in the Magellanic Clouds, although a few massive systems have been identified via MACHO, OGLE or targeted photometric surveys (e.g. Massey et al. 2002; Ostrov & Lapasset 2003; Bonanos 2009 ). By way of example, as part of the TMBM (Almeida et al. 2017 ) follow up to the VFTS survey of OB stars in the Tarantula Nebula (Evans et al. 2011) , only 4 of ∼ 100 O-type binaries were identified as detached eclipsing binaries (Mahy et al. 2020 ).

Double-lined spectroscopic binaries (SB2) are relatively common amongst massive stars and provide accurate orbital parameters, although the inclination must be inferred in order to derive the dynamical masses of the individual components. If the orbital solution is known, SB2 spectra can be disentangled and stellar parameters and chemical compositions obtained for the primary and secondary components. In the LMC there are a few very massive SB2s, the majority of which are located in the Tarantula Nebula (Crowther 2019), including Melnick 34 whose primary component has been estimated to be 139 (Tehrani et al. 2019) , R144 with a primary mass of 74 (Shenar et al. 2021 ) and R139 whose primary exceeds 69 (Mahy et al. 2020) . In principle, inclinations and in turn dynamical masses can be obtained for binaries lacking photometric eclipses via optical (e.g. Shenar et al. 2021) or X-ray light curve modelling. Indeed, Melnick 34 (Pollock et al. 2018 ) was first established as a colliding wind binary via the Chandra X-ray imaging survey of the Marchi et al. 2011) and Chandra ACIS (red, Townsley+ in prep) indicating the target of the current study Mk 33Na and its close visual companions Mk 33Nb (O6.5 V), Mk 33Sa (O3 III(f*), Massey et al. 2005) and Mk 33Sb (WC4+OB, Smith et al. 1990 ). Mk 33Na and Mk 33Sa are bright X-ray sources, while Mk 33Nb and Mk 33Sb remain undetected in >2Ms of Chandra data. The field of view is 8×8 arcsec (2×2 parsec at the distance of the LMC), with North up and East to the left. The rich star cluster R136 lies 3 parsec to the SW.

Tarantula Nebula, T-ReX (Townsley et al., Broos & Townsley, in prep.) . Mk 33Na is located at a projected distance of ∼3.4 pc from the rich star cluster R136 in the Tarantula Nebula (Massey & Hunter 1998; Crowther & Dessart 1998) . Melnick (1985) first identified "Mk 33N" as an early O-type star (O4), with "Mk 33S" ∼2 arcsec to the south west. Massey & Hunter (1998) refined its spectral type to O3 If* from HST/FOS spectroscopy (alias # 16 from Hunter et al. 1995) with "a" added by Crowther & Dessart (1998) to distinguish it from "b", an O6.5 V star (alias # 32 from Hunter et al. 1995) 1.2 arcsec away as shown in Fig. 1 . Other aliases for Mk 33Na include # 1140 in Parker (1993) , # 33 in Selman et al. (1999 ) and # 1943 in Castro et al. (2018 .

Shallow Chandra X-ray observations revealed Mk33 Na as a bright X-ray source (CX9 in Portegies Zwart et al. 2002; # 133 in Townsley et al. 2006) , with a similar X-ray luminosity ( X ) to R139, suggesting the presence of a colliding-wind binary (CWB). Deep X-ray imaging via the T-ReX survey (Townsley et al. in prep.) confirmed Mk33Na as a bright X-ray source ( Fig. 1) corresponding to an attenuated X-ray luminosity of 2.8×10 33 erg s −1 at the distance of the LMC (Crowther et al. in prep.) . On the basis of its Kuiper statistic, a variant of the Kolmogorov-Smirnov (KS) test (Paltani 2004) shown in Figure 2 , Broos & Townsley (in prep.) identified significant X-ray variability suggesting a potential period of 18.3 days for Mk 33Na most likely due to orbital modulation in a binary system. Figure 1 shows that Mk33Sa is similarly bright in X-rays suggestive perhaps of another binary, although it is less variable than its neighbour.

Here we undertake an analysis of VLT/UVES time-series spectroscopy of Mk 33Na motivated by the T-ReX observations, confirm that it is a double-lined spectroscopic binary with a period of ∼18.3 days, estimate its physical and wind properties, and reveal that it comprises a colliding wind binary with a high X / bol = 1.2 × 10 −6 (Crowther et al. in prep.) . In Section 2 we present our orbital analysis of new VLT/UVES spectroscopy and X-ray variability of Mk 33Na, Section 3 presents the spectral types of each disentangled component plus a spectroscopic analysis from which physical, wind and chemical properties are obtained. Stellar masses, ages and orbital inclination are presented in Section 4. We discuss our results in the context of other massive binaries in the LMC in Section 5 and brief conclusions are drawn in Section 6.

Using optical multi-epoch observations (Sect. 2.1) we determine the radial velocity variations of both components of the Mk 33Na system (Sect. 2.2). The orbital parameters are calculated in Sect. 2.3 and the X-ray observations are presented in Sect. 2.4.

To determine the orbital parameters we spectroscopically followed up Mk 33Na (RA 05:38:44.34, DEC -69:05:54.6, J2015.5, GaiaDR2 Gaia Collaboration et al. 2018 with the Ultraviolet and Visual Echelle Spectrograph (UVES, Dekker et al. 2000) mounted at the Very Large Telescope (VLT) in Paranal, Chile. 15 epochs of VLT/UVES spectra were obtained in service mode over a period of 17 days in October 2020 (Program ID 106.21MH.001 aka 5106.D-0835(A), PI Bestenlehner) . The observations were initially planned to sample the orbital period over 2.5 months with preferences on dates where service mode was available. Each epoch or orbital phase could be potentially observed on 2 or 3 dates. According to the Phase 2 guidelines 1 to 2 epochs per orbital phase were scheduled on days when service mode was not available. In the unlikely event that all epochs are obtained in one go including dates where service mode was not available, possible coverage periods were between 17 and 19 days.

Due to Covid-19, regular observations with the VLT paused earlier in 2020. VLT/UVES was one of the first instruments to be brought back in operation for service mode only, from the beginning of October 2020. All 15 epochs were obtained on the nights of 8th through 25th October or in Modified Julian Date (MJD) from 59131.25 to 59148.18 d, even though some nights were originally flagged as no service mode available, as indicated in the observations log (Table 1) . The slit was placed north to south with a slight angle avoiding stars directly in the slit (Fig. 1 ). Observations were taken in Dichroic 2 mode with central wavelengths of 437 and 760nm, providing spectral ranges of 3730 to 4990 Å (blue-arm, EEV 2k×4k detector), 5670 to 7570 Å (lower red-arm, EEV 2k×4k detector) and 7670-9450 Å (upper red-arm, MIT/LL 2k×4k detector) with an integration time of 593 seconds. A 0.8 arcsec slit was used for all epochs with a spectral resolution between 50,000 and 60,000. A signal to noise ratio (S/N) between 30 to 50 was achieved for the blue and lower red-arm.

For the data reduction we used the ESO Reflex pipeline (Freudling et al. 2013) with standard calibration data for service mode. Mk 33Na is located in a crowded region within the Tarantula nebula (Fig. 1) . The default pipeline setting for flux calibrated spectra took the entire slit into account. Even though we placed the slit in a way that minimised contamination from nearby sources, some of the spectra were still contaminated owing to high air mass (up to 2.0) and variable seeing conditions (DIMM 1 0.5 to 1.3 arcsec). The WC star Mk33Sb was the most prominent source which affected He 4686 and the C 5801/5812 spectral lines (Fig. 1) . Therefore, we chose the 2D slit extraction and co-added the regions which were the least contaminated by nearby sources. However, this led to a reduction in the S/N that was achieved (35 and 25 in the blue and lower red-arm, respectively).

In the final step we applied to each epoch the barycentric correction to the wavelength array. To bring the wavelength array of the individual epoch to the same scale we used the interstellar band Ca 3934 for the blue-arm and Na 5890 for the lower red-arm. As a reference value we chose the average radial velocity of the interstellar lines.

Visual inspection showed that N 4058 and C 5801 hardly show any signatures of the secondary star. With the stellar atmosphere and radiative transfer code (Hillier & Miller 1998) we computed a synthetic reference spectrum based on estimated stellar parameters from our spectral classification (Sect. 3.2). Using the synthetic spectrum we performed 2 -minimisation analyses to derive the radial velocities (RV) for N 4058 and C 5801.

The observation taken on 25 October (MJD = 59147.3122d) revealed the largest separation of the He 4542 lines of the primary and secondary (Fig. C2 ). This epoch was used to estimate the optical flux ratio of ∼ 0.25 for the individual components by scaling the templates of the primary and secondary by a factor of 0.8 and 0.2 to match the He line strength. According to the spectral classification we chose a template for the secondary with an effective temperature of 45 000 K. As an initial RV guess for the primary we used the average of the N 4058 and C 5801 RV measurements for each epoch. Subsequently, we performed an iterative 2 -minimisation analysis of He 4542 and H by alternately co-adding and shifting the template spectra of the primary and secondary until convergence. The nebular emission was clipped in the case of H . Lower excitation levels of the Balmer series showed too strong nebular lines and it was not possible to derive RV for all observations. The results are given in Table 1 .

In Fig. 3 we visualise our RV measurements. It clearly shows the anti-phase between the primary and secondary of the Mk 33Na system. The uncertainties of the RV measurements are larger for the secondary and in particular for H . 

To determine the Keplerian orbital parameters we used the Markov Chain Monte Carlo (MCMC) fitting routine developed, utilised and described in Tehrani et al. (2019) . As input prior for the period we used a uniform distribution of 18 ± 2 days according to the Xray period from Fig. 2 . A flat, uninformative prior was adopted for the other parameters. The MCMC fitting was performed using the uncertainties weighted averaged RV measurements from Table 3 for the primary and secondary. The routine delivers a full set of Keplerian orbital parameters: semi-amplitudes of the velocities of the primary ( 1 ) and secondary star ( 2 ), system velocity , eccentricity , longitude of periastron , orbital period , and time of periastron 0 .

In Fig. 4 we show the best-fit RV curve in time and phase space using the most probable parameters from MCMC fitting. The corresponding posterior probability corner plot is given in Fig. 5 and shows that each orbital parameter is well constrained. Within the uncertainties the orbital period agrees well with the derived from the X-ray light-curve (Sect. 2.4). Based on the orbital parameters we are able to infer minimum masses through the standard formula 1,2 sin 3 = /(2 )(1 − 2 ) 1.5 ( 1 + 2 ) 2 2,1 .

From Eq. 1 we are able to derive the mass ratio of the system, = 1 / 2 = 2 / 1 = 0.63 ± 0.02. The eccentricity is low enough that Mk 33Na can be assumed to be a detached binary system. The system velocity of 266.5±1.0 km/s is consistent with the mean value of 267.7 km/s within 5 pc around the stellar cluster R136 (Hénault-Brunet et al. 2012) .

All orbital parameters and minimum masses including uncertainties are given in Table 2 .

The 2 Ms Chandra Visionary Program T-ReX was executed over 630 days between 2014 May 3 and 2016 Jan 22 using the ACIS instrument (16.9 × 16.9 field of view) centred on R136a, the central star at the heart of the Tarantula Nebula cluster (Townsley et al.; in prep) . Here we confirm the 18.3 day orbital period of Mk 33Na from UVES, supporting the X-ray results from the Kuiper statistic (Fig. 2 ) from the T-ReX dataset (Broos & Townsley, in prep) . We present its X-ray light curve folded on the best-fit optical radial-velocity period in Fig. 6 and 7. The T-ReX campaign itself covered 34 orbital cycles suggesting a high degree of reproducibility for an X-ray source of mean count rate close to 1 count per kilosecond. This conclusion is reinforced by the 3 relatively bright measurements taken close to an implied periastron about a decade and 165 orbital cycles earlier in 2006 January. The X-ray sourcelist and source properties from the 2006 Chandra data on 30 Doradus were published in Townsley et al. (2014) . At first sight, Fig. 6 and 7 appear consistent with the luminosity variations expected from a colliding-wind binary system in which the X-ray cooling time exceeds the flow time out of the system leading to a luminosity inversely proportional to the binary separation (e.g. Stevens et al. 1992 ). However, the shape of the observed X-ray spectrum varies with phase as shown in Fig. 8 . The rise in flux in the approach to periastron is mostly due to soft X-rays near 1 keV. Spectral analysis performed in four broad phase bins reveals that the data are consistent with a model in which a separation-dependent luminosity is further modulated by a phase-dependent absorbing column density, where H /10 22 = 1.5, 1.9, 2.1, 1.7 cm −2 in phase bins -0.1:0.1, 0.1:0.3, 0.3:0.7, 0.7:0.9 (Fig. 8) .

Although other models can fit the data, this variable-absorber model offers a simple interpretation. The lowest absorbing column (maximum count rate) occurs when the shocked X-ray emission oc-curring in the space between the stars is viewed through the weaker wind of the secondary O4V star which passes in front across the lineof-sight 1.4 days before periastron, according to the orbital elements in Table 2 . High reddening shows that the interstellar medium accounts for much of the observed X-ray absorption. The time-averaged X-ray luminosity of Mk 33Na is log X = 33.95 ± 0.05 erg s −1 (Crowther et al. in prep.) , a value more typical of binary systems than the lower intrinsic luminosities of single stars.

Before we determine the spectral types (Sect. 3.2) and perform the spectroscopic analysis (Sect. 3.5) we disentangled the spectra into primary and secondary components (Sect. 3.1). Line broadening and line-profile variations are described in Sect. 3.3 and 3.4). The section concludes with Sect. 3.6 where we obtain reddening parameters and luminosities for primary and secondary star.

For an accurate spectral classification and reliable spectroscopic analysis, we disentangled the spectra into primary and secondary components. Individual components exhibit large RV variations which are enhanced by their anti-phase ( Fig. 3 ) making disentanglement relatively straightforward.

Based on the orbital solution from Sect. 2.3, we calculate more precise radial velocities for each epoch of the primary and secondary than measured in Sect. 2.2. In the first step we shifted the epochs to the rest-frame of the primary according to its RV. Co-adding and averaging the spectra would have included signatures of the secondary and for example cosmics which are not removed by 2D slit extraction (Sect. 2.1). Therefore, we calculated a median spectrum which should largely remove the secondary component. The difference between the both methods is shown in Fig. B1 . To obtained the spectrum for the secondary we applied the RV-shift for each epoch to the primary median spectrum from the first step and divided the epoch by the primary spectrum (Fig. B2 ). Then we shifted the epochs to the rest-frame of the secondary and calculated the median spectrum of the secondary.

For an improved disentangling we applied the same methodology Table 1 for the primary and secondary, respectively.

for the primary by dividing all epochs by the secondary median spectrum before calculating the median spectrum of the primary. Those steps were repeated until no changes to the previous iteration were observed. The spectra of Mk 33Na 1 and Mk 33Na 2 are rescaled according to their flux ratio (Sect. 3.6).

A similar methodology was applied to remove the broad C 5801-12 emission caused by the WC4 star Mk33Sb ( Fig. 1 and B3 ). The disentangled spectra are shown in Fig. 11 and 12. A detailed description of the spectra disentangling methodology is given in Appendix B.

Spectral types for the unresolved source Mk 33N from the literature include O4 from Melnick (1985) , with spatial crowding hindering further progress until Massey & Hunter (1998) Fig. 9 ). VFTS 016 and Mk33Na 1 have similar C 4658 line strength, but the nitrogen lines are stronger in VFTS 016 while weaker in the Mk 33Na 1 spectrum (C 4658 > N 4603-20) . The carbon abundances ( C = log(C/H) + 12 = 7.82) of VFTS 016 and Mk33Na 1 are comparable, but the nitrogen abundance ( N = log(C/H)+12 = 7.76) of VFTS 016 is higher by ∼0.5 dex , Table 3 ). Therefore, we classified Mk 33Na 1 as OC2.5 If*.

The spectral type of the secondary (Mk 33Na 2 ) is more difficult to determine. The absence of strong signatures of carbon and nitrogen lines is indicative for dwarfs with a relatively weak stellar wind. Based on the spectral types by Walborn et al. (2014) we noticed that VFTS 797 (O3.5 V((n))((fc))) shows similar spectral properties such as weak C, N and He 4471 lines and relatively strong He 4542 absorption line. N 4058 is not or hardly visible in the spectrum and there might be a hint of N 4634-41 and C 4650. We concluded that the spectral type of Mk 33Na 2 is approximately O4 V.

If the binary system is tidally locked, the mean rotational velocity is 2 eff / based on stellar radius and period of the system (Table 2 and 3). In an eccentric orbit, the projected rotational velocity ( sin ) should vary throughout the orbit and lower or equal to the mean rotational velocity, e.g. 50 km/s for the Mk33Na 1 . We used IACOB-BROAD (Simón-Díaz & Herrero 2014) to derive sin and macro-turbulent velocity ( mac ). IACOB-BROAD is an interactive analysis tool that combines Fourier transformation (FT) and goodness-of-fit methods. The code is publicly available and is written in the interactive data language (IDL).

We were only able to derive the line broadening for the primary Table 3 . Stellar parameters obtained from the spectroscopic analysis. Masses under the assumption of chemically homogeneous evolution ( hom ) are derived with the mass-luminosity relation by Gräfener et al. (2011) . Most probable stellar parameters derived with BONNSAI (Schneider et al. 2014 ) based on stellar evolutionary models by Brott et al. (2011); Köhler et al. (2015) . a log could not be derived, because the wings of the Balmer lines are in emission. b Approximate upper limit for the mass-loss rate. c The probability distribution function shows a wide spread of possible composition, but the most probable value is the initial chemical composition for the LMC (Brott et al. 2011) . based on the C 5801 line. Unfortunately, we were unable to employ metal lines in absorption to determine the line broadening which would provide more accurate measurements. We inferred sin using the first minimum of the FT. The macro-turbulent velocity was determined based on the goodness-of-fit while sin was set to the FT value. The measurements are listed in the Appendix (Table A1 ). Figure 10 shows sin and mac as a function of phase. For reference the RV variations for the primary are included as well. Shortly after periastron both broadening velocities drop below 60 km/s. With the exception of a couple of outliers the average for sin and mac lies around 95 km/s. Averaged sin = 75.5 ± 0.6 and mac = 110.3 ± 4.1 based on the measurements of N and C 5801-12 (Table A1 ). The stellar spectra of the primary show large line profile variations (LPV) including the C 5801-12 doublet (Sect. 3.4). The outliers could be due to the dynamics within the stellar atmosphere and wind. However, we can conclude that the primary is not tidally locked.

For the secondary we were only able to measure sin ≈ 125 km/s and mac ≈ 10 km/s at the orbital phase ∼ 0.1 using He 4542. The He line is also broadened by the Stark effect and could lead to an overestimation of sin and mac . The sin ≈ 125 km/s can be considered as an upper limit for the secondary. sin and mac based on the median spectrum is around 78 and 56 km/s for the C 5812 line (Table A1) .

Line profiles of some lines varied significantly throughout the orbit, e.g. He 4686, N 4058 and wings of Balmer lines. Shortly after periastron He 4686 is weakest and N 4058 has nearly disappeared as well. At phase ∼ 0.1 the line strength suddenly increases and settles down again after 1 day. This could be an outburst of material, but N 4058 is hardly affected. The emission-line strength increases until apastron at orbital phase 0.5 and decreases again with some variation towards periastron (Fig. A1 ).

There are noticeable variations in the wings of the absorption lines H , H , He 4200 and 4542 with filling in of an emission line component. As seen in Fig. 12 those features seem to belong to the secondary. Either the mass-loss rate varies significantly or those LPVs are caused by the stellar wind and ionising flux of the primary star. However, the C 5801-12 doublet is very sensitive to the massloss rate and a small increase would turn the lines from absorption into emission.

A future study will try to understand the nature and physical processes of those variations including the X-ray variability (Sect. 2.4).

The spectra were analysed with the spherical, non-LTE stellar atmosphere and radiative transfer code CMFGEN (Hillier & Miller 1998) . We used the grid of stellar atmospheres from Bestenlehner et al. (2014) to roughly estimate the stellar parameters. Based on the initial guesses we computed a grid with an extended atomic model We were only able to confidently derive stellar parameters for the primary. The effective temperature is based on the presence of C 4658 and the weakness of N 4604-20 relative to N 4058. The surface gravity is derived from wings of H − . The mass-loss rate is determined from the emission line strength of He 4686 which varies throughout the orbit. H was strongly contaminated by nebular lines and was not suitable. Due to CN cycle, the C abundance decreases and N increases during the main-sequence lifetime. However, the primary is still carbon-rich close to the LMC baseline of ( − ) LMC = 0.85 dex (Brott et al. 2011) , with − = ( − ) LMC −0.35±0.4 dex, suggesting a young age for the binary system of less than ∼ 1 Myr.

The determination of stellar parameters for the secondary is less straightforward because of the emission line components in the wings of the Balmer and helium lines. As the log is derived from the broadening of the wings of the Balmer lines, we were unable to determine the surface gravity for the secondary. The absence or the difficult identification of N 4634-41, N 4058, N 4604-20, C 4650, the absorption of C 5801-12 and the presence of He 4471 suggest a eff = 45 000 ± 2500 K for the secondary. C 5801-12 turns from absorption into emission if the mass-loss rate is increased. Other lines are hardly effected and we are only able to determine an upper mass-loss rate limit of log [ / ] < −6.5.

Physical, wind and chemical properties of Mk 33Na are provided in Table 3 . The only previous spectroscopic analysis of the system was undertaken by Castro et al. (2021) 

Bolometric luminosities ( ) of both components are derived by matching optical and near-infrared (near-IR) photometric data with the combined theoretical spectral energy distribution (SED). We Python 2 . We used the nomenclature according to Maíz Apellániz et al. (2014) instead of the more common and ( − ). The combined model SED is inferred using the mass ratio of the Mk 33Na system (Sect. 2.3) and the luminosity -mass relation (LMR) by Gräfener et al. (2011, Eq. 9 with coefficients from row 1 of Table A .1) for chemical-homogeneously core-hydrogen burning stars to calculate the flux ratio of the binary system for a given primary stellar mass (Fig. 13) . The LMR requires a hydrogen surface abundances as input. We assumed the baseline abundances for hydrogen in mass-fraction according to Brott et al. (2011, = 0 .7391), which 2 https://lmfit.github.io/lmfit-py/index.html is supported by the unprocessed C and N abundances (Section. 3.5 and 4).

Under these assumptions the luminosity ratio is 2 / 1 = 0.43 with a flux ratio in the optical of 0.62. The flux ratio is larger than the value of 0.25 assumed for the RV measurements in Sect. 2.2. Looking at Fig. 12 and the spectroscopic fit we can see that the Balmer and helium absorption lines of the secondary are filled in by an emission line component.

Reddening parameters and total bolometric luminosity of the Mk 33Na system are given in Table 4 while individual for the primary and secondary are listed in Table 3 . The error provided for the combined only represents the uncertainties of the reddening parameters. 5495 and (4405 − 5495) are typical values for the stellar cluster R136 and its immediate surroundings (e.g. Doran et al. 2013; . Armed with log X = 34.01 ± 0.05 erg s −1 from Sect. 2.4, we confirm that Mk 33Na has an usually high X / Bol = 10 −5.9 consistent with excess X-ray emission from wind-wind collisions.

Our spectroscopic results for Mk 33Na favour masses of 106 and 66

for the primary and secondary based on the LMR for chemical homogeneously evolving hydrogen burning stars of Gräfener et al. (2011) . Therefore, those masses can be seen as an upper mass-limit for the Mk 33Na system. This can be reconciled with results from our orbital solution, 1,2 sin 3 = 20.0 and 12.6 if the orbital inclination of Mk 33Na is ∼ 35 • .

For LMC metallicity and He baseline chemical abundances the star approaches the Eddington limit at an Eddington parameter Γ e ≈ 0.7 considering only the electron opacity e . With the assumption that a star cannot exist above the Eddington limit we derive lower mass limits of ∼ 39 and ∼ 17 for Mk 33Na 1,2 . Taking the dynamical mass ratio into account the minimum for Mk 33Na 2 would correspond to ∼ 25 . This results in a maximum orbital inclination of < 53 • which is a true upper limit as stars in proximity to Eddington limit experience enhanced mass loss with strong emission lines in their spectra and would be of spectral type WN rather than Of* (e.g. Bestenlehner et al. 2014; . We conclude that the orbital inclination is in the range between 35 and 53 • .

We have also estimated the physical properties of the Mk 33Na components by utilising BONNSAI 3 Schneider et al. (2014) which applies a Bayesian approach to the determination of stellar masses and ages based on evolutionary models, for which we employ LMC metallicity models from Brott et al. (2011) and Köhler et al. (2015) . BONNSAI results are presented in Table 3 and reveal evolutionary masses of 83 and 48 for the components of Mk 33Na plus individual ages of 0.9±0.6 and 1.6 +0.7 −1.3 Myr, such that the favoured age of the system is ∼1.1 Myr. Evolutionary masses favour an orbital inclination of ∼ 38 • . Alternate evolutionary masses and ages based on the optical flux ratio solution are also presented in Table 3 , 3 https://www.astro.uni-bonn.de/stars/bonnsai/ from which a similar age is obtained, albeit with 91 and 42 for systemic components, contrary to the dynamical mass ratio.

Near Mk 33Na is the star cluster R136. By analysing the stellar population of R136 Bestenlehner et al. (2020) discovered a positive mass-discrepancy, where the spectroscopic masses are systematically larger than evolutionary masses for stars more massive than ∼ 40 . Based on a sample of eclipsing binaries with known inclination Mahy et al. (2020) showed that spectroscopic masses are in better agreement with dynamical than evolutionary masses. Therefore, the actual masses of both components might be larger than the values obtained above with BONNSAI.

Unusually for an O supergiant, the primary in Mk 33Na has an OC spectral type since C 4658 > N 4603-20. The LMC baseline abundance difference between carbon and nitrogen is ( − ) LMC = 0.85 dex, with − = ( − ) LMC −0.35±0.4 dex obtained from our spectroscopic analysis. The favoured BONN-SAI solutions involve no chemical enrichment for Mk 33Na, as summarised in Table 3 .

Historically HSH95 38 (O3 V+O6 V) was the only massive eclipsing binary system in the LMC with a primary whose mass exceeded ∼ 50 (Massey et al. 2002) . A number of other eclipsing binary systems have been discovered via OGLE, MACHO or other photometric surveys, including Sk -67 • 105 (Ostrov & Lapasset 2003) and W61 28-22 (Bonanos 2009 ). The advent of multi-epoch optical spectroscopic surveys involving hundreds of OB stars, such as VFTS (Evans et al. 2011) , and subsequent follow up studies (TMBM, Almeida et al. 2017) , has led to the discovery of several non-eclipsing systems involving O supergiants or H-rich WN stars such as R139 (Taylor et al. 2011; Mahy et al. 2020 ) and R144 (Shenar et al. 2021) . A compilation of massive binary systems is shown in Fig. 14 and listed in Table 5 .

Optical spectroscopic studies are ideal methods of identifying relatively short period binaries with high inclinations. X-ray monitoring, such as T-ReX (Broos & Townsley, in prep.) , provides an alternative approach to identifying binaries via the detection of excess X-ray emission, arising from wind-wind collisions (Pollock 1987) . Crowther et al. (2021) have shown that SB2 systems have higher / Bol ratios than single or SB1 systems detected in T-ReX. T-ReX has previously identified Mk34 as an exceptional colliding Table 5 . Physical properties of double-lined spectroscopy binaries in the LMC whose minimum primary mass exceeds ∼ 20 . These have been identified from optical photometric eclipses (Phot: e.g. MACHO, OGLE), optical spectroscopic radial velocity variability (Spec: e.g. VFTS, Evans et al. 2011; TMBM: Almeida et al. 2017) or X-ray variability (X-ray: Broos & Townsley, in prep . Here, we have confirmed the suspicion from T-ReX that Melnick 33Na is another massive short-period (18.3 day) binary comprising an OC2.5 If* primary with an evolutionary mass of 83±19 and an O4 V secondary with a mass of 48 ± 11

. Table 5 provides a summary of the most massive double-lined spectroscopic binaries in the LMC, for which the primary of Mk 33Na has the (joint) third highest mass according to spectroscopic/evolutionary models, after Mk34 and R144, and comparable to R139. Figure 14 shows a Hertzsprung-Russell diagram (HRD) of binary systems in the LMC, where stellar parameters are available for both components. Mk33 Na 1 has the earliest O star spectral type known in a binary system. In the HRD the star is located in a region which is populated by Of/WN and WNh stars (Shenar et al. 2017; Tehrani et al. 2019; Shenar et al. 2019 Shenar et al. , 2021 . The only O stars with a similar luminosity is the O supergiant binary system R139 (Taylor et al. 2011; Mahy et al. 2020 ), but of later spectral type. Mk33 Na 1 still has a high carbon abundance at the surface while stars of similar luminosity show chemical enrichment due to the CNO-cycle, N enriched and C deficient. Some of those luminous stars are even He enriched or potentially He core-burning. This supports the young age of ∼1 Myr for the Mk 33Na system.

Looking at the top and high luminosity end of the HRD (Fig. 14 ) and the stellar masses and mass ratio of those objects (Table 5) it is apparent that most of those binary systems have close to unity and consist of stars with similar luminosities which is likely an observational bias. Binaries can be easier identified if either both components have similar luminosities or the spectrum has high SNR and high spectral resolution (SB2). Photometrically eclipsing binaries with high inclinations or a large monitoring program such as TMBM are able to identify and characterise binary systems over a wide range of mass ratios with between 0.35 and 1 (e.g. Massey et al. 2002; Bonanos 2009 ; Mahy et al. 2020, Table 5 ). Those circumstances have left Mk 33Na undetected as a binary system for several decades. Only now has the X-ray variability due to colliding winds finally revealed the multiple nature of Mk 33Na.

We establish Melnick 33Na (Mk 33Na) in the Tarantula Nebula as a double-lined spectroscopic binary, from an analysis of VLT/UVES time-series spectroscopy, with the following properties • An orbital period of 18.3±0.1 days, supporting the period estimate from analysis of X-ray observations obtained with the T-ReX survey, eccentricity of ∼0.33 with a mass ratio of = 2 / 1 = 0.63.

• A primary component with spectral type OC2.5 If* with a minimum dynamical mass of 20 , although a much higher mass obtained from spectroscopic or evolutionary models (83±19 and an age of 0.9±0.6 Myr via BONNSAI) since our favoured spectroscopic luminosity is log / = 6.15 ± 0.18. Unusually, − = ( − ) LMC − 0.35 ± 0.4 dex, also supporting a young age for the system of ≤1 Myr.

• A secondary component with spectral type O4 V with a minimum dynamical mass of 12.6

, and again a substantially higher mass from spectroscopic or evolutionary models (48±11 and an age of 1.6 +0.7 −1.3 Myr via BONNSAI) since our favoured spectroscopic luminosity is log / = 5.78 ± 0.18. • Evolutionary mass estimates of the primary of Mk 33Na rank it (joint) third in terms of the most massive SB2 systems in the LMC after Mk 34 (WN5h+WN5h) and R144 (WN5-6+WN6-7h) and comparable to R139 (O6.5 Iafc+O6 Iaf). Evolutionary and dynamical solutions can be reconciled if ∼ 38 • . Mk 33Na is confirmed as a colliding wind binary with excess X-ray emission from T-ReX and X / Bol = 10 −5.9 . In the following disentangling work-flow we ignore any crosscontamination by nearby sources and only consider the fluxes of the primary and secondary. The SED of each components can be written as the product of relative flux and continuum SED (SED 1,2 = 1,2 1,2 ). The combined flux of the secondary and primary are recorded with VLT/UVES and after normalisation we ob- tained the following spectrum for each epoch i:

In the first step we shifted the epochs to the rest-frame of the primary and calculate the median spectrum:

In the case of shift and add (e.g. González & Levato 2006) the median is replaced with mean. A comparison is shown in Fig. B1 .

To obtain the spectrum of secondary we divide each epoch by the median primary spectrum (Med 1 ). We could choose a subtraction instead of division, but this could add numerical noise to the spectrum because computers can only represent numbers with limited precision and 2 nearly equal subtracted numbers become indistinguishable. A division leads to the following equation and secondary median spectrum (Med 2 ):

which has the same form as the results of the primary in the first step (Eq. B2). An example of an epoch spectrum after division by the primary median spectrum is given in Fig. B2 . The primary and secondary median spectra were iteratively used to divide each epoch before the median was calculated. The iteration was stopped, if no changes were observed to the previous iteration. The median spectra of both components are rescaled according to their flux ratio (Sect. 3.6). We divided the median primary spectrum by relative Flux Figure B1 . Blue solid line is the primary median spectrum while red dashed line is the averaged spectrum over all epochs (shift and add). Both spectra are similar, but the averaged spectrum (red dashed line) contains more artefacts. or 2 /( 1 + 2 ) for the secondary and than subtract 2 / 1 (primary) or 1 / 2 (secondary):

The broad C 5801-12 emission caused by the WC4 star Mk33Sb was removed in a similar way by just calculating the median without shifting the epochs (Fig. B3) . We do not know the actual relative flux of Mk33Sb to the combined Mk33Na system as the light of Mk33Sb leaked into the slit depending on seeing and airmass. In addition, we tried to minimise its contamination during the reduction process by using the 2D slit extraction (Sect. 2.1). A rescaling might not make a noticeable difference.

The disentangled spectra for Mk 33Na 1 and Mk 33Na 2 are shown in Fig. 11 and 12 . 

This paper has been typeset from a T E X/L A T E X file prepared by the author. 

Optical and IR Telescope Instrumentation and Detectors

We thank the referee, Paco Najarro, for detailed and helpful comments which improved the clarity and content of the manuscript. JMB would like to thank STFC for financial support through STFC consolidated grant ST/V000853/1. PSB and LKT are supported by the Penn State ACIS Instrument Team Contract SV4-74018, issued by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060.

The data underlying this article are available in the article and in its supplementary material. Spectroscopic data will be available via the ESO archive facility while synthetic spectra can be requested from the lead author. X-ray data are available in the Chandra Data Archive and the Chandra Source Catalogue and will be published in Townsley et al., Broos et al., in preparation.