The First Dynamical Mass Determination of a Nitrogen-rich Wolf–Rayet Star Using a Combined Visual and Spectroscopic Orbit

The most fundamental parameter for a star is its mass, which is only accurately measured through application of Kepler's Laws in binary systems. While eclipsing binaries have been the standard for these analyses, some rare types of stars are not yet known to have members in eclipsing systems. In the case of massive stars, it has been shown that most massive stars reside in multiple systems (see, e.g., Mason et al. 2009; Sana et al. 2014; Aldoretta et al. 2015), and it has become an observational fact that binary evolution can dominate the end results of these stars (Sana et al. 2012; de Mink et al. 2013). In the case of massive stars, the correct interpretation of binary stellar populations could explain the observed SEDs of galaxies throughout history, as stripped stars in binary systems change the ultraviolet flux in a galaxy. However, for this to work properly, we need to have systems for which we can pinpoint stellar parameters and constrain binary evolution (Eldridge et al. 2017).

Some of the best examples of binaries showing evidence of past mass transfer are ones with classical Wolf–Rayet (WR) stars. These stars have lost their hydrogen envelopes through the processes of stellar winds and, possibly, binary interactions.

This analysis focuses on WR 133, a relatively under-studied binary classified as WN5o + O9I by Smith et al. (1996), where the "o" suffix denotes no measurable hydrogen in the WR spectrum and the "N" denotes it as nitrogen-rich. The system is a bright member of NGC 6871 (Rosslowe & Crowther 2015), a cluster residing at a distance of 2.14 ± 0.07 kpc according to Malchenko & Tarasov (2009). The system is used in the calibration of WR star parameters by Rosslowe & Crowther (2015). The binary orbit of the system was best studied by Underhill & Hill (1994), who found a 112.4 day period, a moderate eccentricity of 0.39, and a clear double-lined status. The system shows some polarimetric variability due to changing observed wind geometry throughout the orbit (Robert et al. 1989), which was recently studied by Fullard et al. (2020), who found a lower limit to the inclination to be 1159. Systems with closer to face-on inclination angles are difficult to constrain using only polarimetric variability. The associated masses for the 1159 inclination are unrealistically low, meaning that the system would likely be closer to face-on, and perhaps could be resolved with long-baseline interferometry.

In this paper, we present the first long-baseline optical interferometry of WR 133, which spatially resolves the binary and allows us to compute a combined astrometric and spectroscopic orbit. In Section 2, we present our interferometric observations and the spectroscopy used in this analysis. In Section 3, we describe the measurements made and present the orbital elements of the system. We present a combined spectroscopic model in Section 4, and then discuss these results and conclude this study in Section 5.

2.1. Spectroscopy

2.2. CHARA Array Interferometry

Interferometric observations were obtained with the CHARA Array, which has six 1 m telescopes in a Y-shaped array with baselines ranging from 34 to 330 m (ten Brummelaar et al. 2005) and the upgraded MIRC-X combiner (Monnier et al. 2004; Kraus et al. 2018; Anugu et al. 2018, 2020) in 2019 and 2020. The observations were recorded in the PRISM50 mode, which provides a spectral resolution of R = 50 in the H-band between 1.4 and 1.7 μm. The data were reduced using the MIRC-X data reduction pipeline, version 1.2–1.3¹⁵ to produce calibrated visibilities and closure phases. During the reduction, we applied the bispectrum bias correction included in the pipeline and set the number of coherent coadds to 5. Increasing the number of coherent coadds can improve the signal-to-noise, but averaging over too long a time in comparison to the atmospheric seeing can bias the calibration. In general, the difference between the binary solutions based on the 5 and 10 coherent coadd reductions was smaller than the formal errors, and the number of coherent coadds of 5 reduction typically gave a lower ${\chi }^{2}$ value in the binary fit. We applied minimum uncertainties of 5% to the visibilities and 03 to the closure phases based on the typical calibration uncertainties expected for MIRC-X data.

The calibrators observed on each night are listed in Table 1. We adopted uniform disk diameters in the H-band for the calibrators from the JMMC stellar diameter catalog (Bourges et al. 2017). We inspected the visibilities and closure phases of the calibrators to check for binarity, and also examined the results from the automated pipeline (Anugu et al. 2020), which calibrates the calibrators against each other and uses CANDID (Gallenne et al. 2015) to determine companion detection limits. We did not find any reliable binary detections for the calibrators. Due to the proximity in the sky to WR 140, we used some of the same calibrator stars as given in the analysis of Thomas et al. (2021), with the diameter estimates taken from Bourges et al. (2017). All calibrators and data sets are tabulated in Table 1. The closure phase calibration is well tested by Anugu et al. (2020), so that we are confident that the sign of the closure phases are correct. The reduction pipeline is discussed in detail by Anugu et al. (2020), and we refer the reader to that publication for a detailed discussion of the reductions.

Table 1.  Calibrator Stars Used for the CHARA Data

Night	${N}_{\mathrm{datasets}}$	Baselines	Calibrator	${\theta }_{\mathrm{UD}}$ (mas)
2019 Jul 01	1	E1-W2-W1-S2-E2	HD 178538	0.2487 ± 0.0062
			HD 191703	0.2185 ± 0.0055
			HD 197176	0.2415 ± 0.0058
			HD 201614	0.3174 ± 0.0074
2019 Jul 02	2	E1-W2-W1-S2-S1-E2	HD 178538	0.2487 ± 0.0062
			HD 191703	0.2185 ± 0.0055
			HD 197176	0.2415 ± 0.0058
2019 Sep 04	1	E1-W2-W1-S2-S1-E2	HD 170585	0.3576 ± 0.0087
			HD 170977	0.3409 ± 0.0082
			HD 191703	0.2185 ± 0.0055
			HD 197176	0.2415 ± 0.0058
			HD 201614	0.3174 ± 0.0074
2019 Sep 05	2	E1-W2-W1-S2-S1-E2	HD 170585	0.3576 ± 0.0087
			HD 170977	0.3409 ± 0.0082
			HD 178538	0.2487 ± 0.0062
			HD 191703	0.2185 ± 0.0055
			HD 197176	0.2415 ± 0.0058
			HD 201614	0.3174 ± 0.0074
2019 Sep 13	1	E1-W1-S2-S1-E2	HD 191703	0.2185 ± 0.0055
			HD 197176	0.2415 ± 0.0058
2020 Jun 23	1	E1-W2-W1-S2-S1-E2	HD 178538	0.2487 ± 0.0062
			HD 191703	0.2185 ± 0.0055
			HD 197176	0.2415 ± 0.0058
			HD 201614	0.3174 ± 0.0074

Download table as:  ASCIITypeset image

For the WR star, we measured radial velocities of the He ii λ5411 emission line through bisecting the line (see Figure 1). We measured the bisector at five different heights above the normalized continuum, with levels chosen between 3% and 5% above the continuum. For each level, the velocity on the red and blue wing was interpolated over a smoothed emission line, and then averaged between the two line wings. The average velocity was then used for the measurement. Errors were estimated and then combined in quadrature for both the standard deviation of these measurements and the differences between the measured velocities of both interstellar Na D absorption lines. We reanalyzed all the data reported by Underhill & Hill (1994) for consistency of the measurements, and because the former measurements were made with Gaussian fits of the line, which could cause discrepancies between the observations. All velocities are tabulated in an online table.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The O star velocities were measured through Gaussian fits to the He ii λ5411, O iii λ5592, and He i λ5876 absorption lines. The errors for each individual measurement were estimated through the combination in quadrature of the errors of the centroid position and the error in the wavelength calibration, estimated by the difference in velocity for the fit of the two interstellar Na D absorption lines. The final averaged error was then calculated by combining the standard deviation of these measurements and the maximum individual error of the measurements in quadrature.

The relative astrometry of the binary was measured by using the calibrated interferometric data with the same approach as Richardson et al. (2016a) and Thomas et al. (2021). This uses an adaptive grid search to find the best-fit binary position and flux ratio using software described by Schaefer et al. (2016).¹⁶ During the binary fit, we fixed the uniform disk diameters based on the Gaia distance and expected radius of the O star. For the O9I star, we adopt ${R}_{* }=22.6{R}_{\odot }$ based on the calibration by Martins et al. (2005), and adopt the same radius for the WR star. While we have better radii from our models (Section 4), the adopted angular diameters were 0.11 mas, which can be considered point sources with the CHARA Array and have little effect on the derived measurements.

The uncertainties in the binary fit were derived by adding in quadrature errors computed from four sources: the formal covariance matrix from the binary fit, the variation in parameters when changing the coherent integration time used to reduce the data (5 and 10 coherent coadds), and the variation in parameters when changing the wavelength scale by 0.5% (the formal uncertainty reported by Anugu et al. (2020), for the MIRC-X instrument) as well as the variation in binary parameters when changing the visibility calibration by 5%. In scaling the uncertainties in the position, we added these values in quadrature for the major axis of the error ellipse (${\sigma }_{\mathrm{major}}$) and scaled the minor axis (${\sigma }_{\mathrm{minor}}$) to keep the axis ratio and position angle fixed according to the values derived from the covariance matrix. The results of the astrometric measurements are given in Table 2, with significant figures dependent on the individual measurements. In addition to the previously discussed parameters, we include the position angle of the error ellipse (${\sigma }_{\mathrm{PA}}$) in Table 2.

Table 2.  Astrometric Measurements Derived from CHARA Data

UT	MJD	θ	ρ	Error Ellipse (${\sigma }_{\mathrm{PA}}$)			Flux Ratio
Date		(°)	(mas)	major (mas)	minor (mas)	PA (mas)	${f}_{\mathrm{WR}}/{f}_{{\rm{O}}}$
2019 Jul 01	58665.310	272.069	0.9152	0.0048	0.0032	155.3	0.2850 ± 0.0027
2019 Jul 02	58666.222	270.154	0.9184	0.0070	0.0030	95.4	0.2664 ± 0.0088
2019 Jul 02	58666.417	270.013	0.9050	0.0046	0.0039	21.4	0.2859 ± 0.0022
2019 Sep 04	58730.273	356.868	0.8832	0.0052	0.0035	99.4	0.2774 ± 0.0046
2019 Sep 05	58731.255	354.493	0.8986	0.0050	0.0042	69.2	0.2823 ± 0.0038
2019 Sep 05	58731.320	354.341	0.8979	0.0049	0.0025	99.8	0.2788 ± 0.0011
2019 Sep 13	58739.226	339.177	0.9797	0.0058	0.0017	161.5	0.2835 ± 0.0181
2020 Jun 23	59023.262	217.971	0.6811	0.0050	0.0043	62.2	0.2831 ± 0.0059

Download table as:  ASCIITypeset image

Underhill & Hill (1994) found that the systemic velocities for the orbit of the Wolf–Rayet and O star components differed by about 100 km s⁻¹. In order to minimize the errors present in our fit, we began by fitting the orbit of each component star with the method of Morbey & Brosterhus (1974). This allowed us to then remove the systemic velocity from the radial velocity measurements we made in order to combine all astrometric and radial velocity measurements into a combined orbital fit. The fitting routine we use was also described by Schaefer et al. (2016). We found that fitting the orbit resulted in a large reduced ${\chi }^{2}$-statistic. We were able to get this to a reasonable value by increasing the errors of the WR radial velocity measurements by a factor of 2. The errors of the O-star velocities and the astrometry yield orbital fits with a reduced ${\chi }^{2}$-statistic of less than unity when considered as a fit with a visual orbit and an SB1, so our errors in the orbit seem to be dominated by the errors in the WR radial velocities. We also note that Underhill & Hill (1994) found that the time of periastron passage was different for the different component stars when fitting the orbit. This has been observed in other WR binaries and comes from colliding wind excess emission altering the emission line profiles from the WR star, and could explain our large reduced ${\chi }^{2}$-statistic from the combined fit.

The orbital elements from the simultaneous fit to the visual orbit and spectroscopic radial velocities are given in Table 3 in the column labeled VB+SB2. These parameters are consistent with the orbital elements of Underhill & Hill (1994), which are included in the last column of Table 3 for reference, and our fit is shown in Figure 2. With a visual orbit, we can combine the spectroscopic parameters of a_WR sin i and a_O sin i to compare with the angular semimajor axis from interferometry. This yields an orbital parallax yielding a distance of 1.73 ± 0.17 kpc, which is in agreement with the Gaia DR2 distance of ${1.85}_{-0.14}^{+0.16}$ kpc (Rate & Crowther 2020) and the Gaia early DR3 results that imply the same distance with a smaller uncertainty of 1.86 ± 0.08 kpc (Gaia Collaboration et al. 2020). The dynamical masses computed from M_WR sin³ i, M_O sin³ i, and the inclination are M_O = 18.1 ± 6.5 M_⊙ and M_WR = 8.2 ± 2.8 M_⊙. These masses have large uncertainties due to the error on the inclination and the nearly face-on system geometry and the uncertainty in the semiamplitudes. We can reduce the uncertainties in the masses by using the Gaia EDR3 distance to compute the total mass through Kepler's third law and combine with the spectroscopic mass ratio (K_O/K_WR) to get slightly larger masses of M_O = 22.5 ± 2.5 M_⊙ and M_WR = 10.2 ± 1.3 M_⊙ (last two rows of Table 3). Lastly, in this table, we also include the Roche radii of the component stars from the approximation of Eggleton (1983), showing that these stars are not filling their Roche lobes.

In Table 3, we also present a solution where we constrain the orbital parallax to agree with the Gaia distance of 1.85 kpc. We varied the distance within the 1σ uncertainties, recomputed the orbit fits, and added the difference in the parameters in quadrature with the formal errors from the fitting routine. The orbital parameters at the Gaia distance are consistent with the fit where the distance is not constrained. We adopt this solution, combined with the uncertainties from the Gaia ERD3 parallax, to produce final masses of M_O = 22.6 ± 3.2 M_⊙ and M_WR = 9.3 ± 1.6 M_⊙.

Table 3.  Orbital Elements

Orbital Element	VB+SB2	At Gaia Distance	Underhill & Hill (1994)
Elements of the System
Period (day)	112.780 ± 0.036	112.736 ± 0.073	112.4 ± 0.2
T (periastron, MJD)	58701.58 ± 0.38	58701.58 ± 0.56	47420.5 ± 3.6 (O-star only)
Eccentricity	0.3558 ± 0.0050	0.3646 ± 0.0103	0.39 ± 0.07
${\omega }_{{\rm{O}}}$ (°)	225.3 ± 6.1	238.9 ± 14.9	$198.9\pm 10.$
Elements of the Visual Orbit
a (visual, mas)	0.7863 ± 0.0060	0.7791 ± 0.0024	⋯
Inclination (°)	160.44 ± 1.86	162.08 ± 1.74	⋯
Ω (°)	171.5 ± 6.5	186.1 ± 15.6	⋯
Elements of the Spectroscopic Orbit
${K}_{{\rm{O}}}$ (km s−1)	14.63 ± 1.51	13.41 ± 2.34	16.9 ± 2.1
${K}_{\mathrm{WR}}$ (km s−1)	32.30 ± 3.02	32.56 ± 3.98	34.4 ± 7.4
${\gamma }_{{\rm{O}}}$ (km s−1)	−15.09 ± 0.48	−15.09 ± 0.48	−20.9 ± 0.7
${\gamma }_{\mathrm{WR}}$ (km s−1)	78.1 ± 3.0	78.1 ± 3.0	70.2 ± 4.6
${\chi }_{\mathrm{red},\mathrm{fit}}^{2}$	1.37	1.07	⋯
Derived Properties
${a}_{\mathrm{WR}}$ (au)	0.934 ± 0.122	1.021 ± 0.157	⋯
${a}_{{\rm{O}}}$ (au)	0.423 ± 0.058	0.421 ± 0.083	⋯
${a}_{\mathrm{WR}}$ (${R}_{\odot }$)	200.7 ± 26.2	219.5 ± 33.7	⋯
${a}_{{\rm{O}}}$ (${R}_{\odot }$)	90.9 ± 12.5	90.5 ± 17.8	⋯
Distance (kpc)	1.73 ± 0.17	1.85 (fixed)	⋯
${M}_{{\rm{O}}}$ (${M}_{\odot }$)	18.1 ± 6.5	22.3 ± 9.4	⋯
${M}_{\mathrm{WR}}$ (${M}_{\odot }$)	8.2 ± 2.8	9.2 ± 4.0	⋯
${R}_{\mathrm{roche},\mathrm{WR}}$ (${R}_{\odot }$)	58.8	60.1	⋯
${R}_{\mathrm{roche},{\rm{O}}}$ (${R}_{\odot }$)	84.3	90.1	⋯
Adopting Gaia EDR3 distance of 1.86 ± 0.08 kpc
MO (M⊙)	22.5 ± 2.5	22.6 ± 3.2	⋯
MWR (M⊙)	10.2 ± 1.3	9.3 ± 1.6	⋯

Download table as:  ASCIITypeset image

To derive the physical parameters of both stellar components of WR 133, we utilize the Potsdam Wolf–Rayet (PoWR) model atmosphere code (Hamann & Gräfener 2003; Gräfener et al. 2002; Sander et al. 2015) to analyze the available FUSE, IUE, and optical spectra. PoWR is a 1D code that solves the radiative transfer problem in spherical geometry, and is applicable to any hot star with an expanding atmosphere (e.g., Sander et al. 2014; Shenar et al. 2015; Ramachandran et al. 2019). The analysis takes advantage of precalculated grids for Galactic WN stars (Todt et al. 2015) and OB-type stars (Hainich et al. 2019), which enable an estimation of the errors and provide good starting models from which tailored models can be constructed. The analysis of a binary is possible by combining two PoWR models such that the sum of their spectra reproduce the observations. In the following, we briefly describe the concepts of our PoWR binary analysis. A more detailed description of the methodology is given in Shenar et al. (2016, 2019).

A PoWR model is primarily defined by its temperature ${T}_{* }$, gravity g, luminosity L, mass-loss rate $\dot{M}$, terminal velocity ${v}_{\infty }$, and chemical abundances. The temperature ${T}_{* }$ refers to the effective temperature at the innermost layer of the model, which is defined at a mean Rosseland optical depth of ${\tau }_{\mathrm{Ross}}=20$. For OB-type stars, it is almost identical to the effective temperature ${T}_{2/3}$, defined at ${\tau }_{\mathrm{Ross}}=2/3$. For WR stars, the photosphere can be significantly above the hydrostatic layers, which are concealed by the stellar wind. The stellar radius ${R}_{* }$ relates to ${T}_{* }$ and L via the Stephan–Boltzmann equation $L\propto {R}_{* }^{2}{T}_{* }^{4}$.

A helpful parameter used to described PoWR models is the so-called transformed radius ${R}_{{\rm{t}}}$, which characterizes the strength of emission recombination lines:

$$\begin{eqnarray}&&{R}_{{\rm{t}}}={R}_{* }{\left[\displaystyle \frac{{v}_{\infty }}{2500\mathrm{km}{{\rm{s}}}^{-1}}/\displaystyle \frac{\dot{M}\sqrt{D}}{{10}^{-4}{M}_{\odot }{\mathrm{yr}}^{-1}}\right]}^{2/3}.\end{eqnarray} \tag{ 1 }$$

Models with given ${T}_{* }$ and ${R}_{{\rm{t}}}$ values will tend to exhibit a similar emission spectrum irrespective of L, ${v}_{\infty }$, and D.

The velocity field $v(r)$ in the model smoothly connects the subsonic regime, where the density is determined from quasi-hydrostatic equilibrium (Sander et al. 2015), and the supersonic regime, where the density follows from the wind velocity described by a $\beta$-law (Castor et al. 1975). The winds are assumed to be clumped. Clumping is treated in the microclumping approximation (Hamann & Koesterke 1998) and is described by the clumping factor D, which gives the density ratio between a clumped wind and the equivalent smooth wind. The radial clumping stratification is as described by Shenar et al. (2015) for the O-type component, and is assumed to be constant for the WR component.

The line profiles are characterized by the Doppler velocity ${v}_{\mathrm{Dop}}(r)=\sqrt{{v}_{\mathrm{th}}^{2}+{\xi }^{2}}$, where ξ is the microturbulent velocity and ${v}_{\mathrm{th}}$ is the thermal velocity (Shenar et al. 2015). The microturbulence is depth-dependent, starting from photospheric values of 100 and 14 km s⁻¹ for the WR and O-type components, respectively (Todt et al. 2015; Hainich et al. 2019), and scaling with the wind velocity as $\xi (r)=0.1\cdot v(r)$ beyond.

We include the chemical species H, He, C, N, O, and the iron group elements for both components, and add Mg, Si, P, and S for the O-type secondary. The abundances of heavy elements are fixed to solar (Asplund et al. 2009), but the H, C, N, and O abundances are examined as free parameters in the analysis.

An important quantity in a binary analysis is the flux ratio of both components in a certain band. In principle, it can be estimated from the relative strengths of the O-star features in the spectrum, which would become more or less diluted depending on the adopted light ratio. However, the interferometric measurements described in Section 3 provide us with a direct measure of the light ratio in the H-band of ${f}_{\mathrm{WR}}/{f}_{{\rm{O}}}(H)=0.28$. We fix the light ratio to this measured value.

The effective temperatures ${T}_{* }$ are derived by reproducing the ionization balance of different ionization stages of a given species, primarily N iv and N v for the WR primary and He i and He ii for the O-type secondary. With the fixed light ratio, the wind parameters of the WR primary are determined from the line strength and shape of P-Cygni lines in the UV (most prominently N v $\lambda \lambda 1239,1243$, and C iv $\lambda \lambda 1548,1551$) and the optical recombination lines. The wind parameters of the O-type secondary can also be derived accurately, owing to the presence of significant Hα emission, along with the resonance P-Cygni lines N v $\lambda \lambda 1239,1243$, Si iv $\lambda \lambda 1394,1403$, and C iv $\lambda \lambda 1548,1551$. The gravity of the O-type secondary is determined from its Balmer absorption wings. While they are entangled with the WR emission, the gravity can be estimated due to the now fixed WR wind parameters and is found to be typical for O-type supergiants or bright giants (luminosity classes I or II, e.g., Martins et al. 2005). Lacking photospheric features, the gravity of the WR component cannot be derived. Finally, the abundances are determined from the overall strength of corresponding spectral lines.

The projected rotational velocity $v\sin i$ and macroturbulent velocity ${v}_{\mathrm{mac}}$ of the O-type secondary are determined by utilizing the iacob-broad tool (Simón-Díaz & Herrero 2007, 2014), which relies on the Fourier technique and goodness-of-fit analysis. We use the isolated O iii $\lambda 5592$ line for this purpose. The models are convolved with a rotation profile and a radial-tangential profile of the derived Doppler widths (Gray 1973). The values agree with a qualitative inspection between the synthetic and observed line profiles. For the WR primary, $v\sin i$ and ${v}_{\mathrm{mac}}$ cannot be derived.

For the WR component, we find no clear indications for β values that significantly differ from unity, although we note that there is a degeneracy between β and the other parameters. We therefore adopt the typical $\beta =1$ (Todt et al. 2015). For the O-type component, we find indications for β values on the order of 3. Lower β values overestimate the width of the Hα line, while larger β values underestimate it. The terminal velocity is fixed from the P-Cygni lines, and hence β is altered to provide a qualitative agreement. The clumping contrasts can only be estimated via electron scattering wings for the WR primary and the relative strength of the P-Cygni and Hα line for the O-type secondary. They are found to be typical for the respective spectral types (e.g., Hamann & Koesterke 1998; Bouret et al. 2012).

We find very clear indications that X-ray radiation impacts the UV spectrum of the O-type secondary (Figure 3). Through K-Shell transitions (Auger ionization), X-rays can alter the ionization balance in the wind, most prominently populating levels associated with high-ionization lines such as N v $\lambda \lambda 1239,1243$ (Cassinelli & Olson 1979; Oskinova et al. 2011). X-rays are modeled as described in detail in Baum et al. (1992) and Shenar et al. (2015). In this approach, the X-ray emission is assumed to originate in optically thin filaments of shocked plasma embedded in the wind, characterized by three parameters: the X-ray temperature ${T}_{{\rm{X}}}$, the filling factor ${X}_{\mathrm{fill}}$, and the onset radius R₀. We set ${T}_{{\rm{X}}}$ to 1 MK and the onset radius to ${R}_{0}=2\,{R}_{* }$. The filling factor is adjusted such that the observed X-ray luminosity roughly agrees with the observed value of ${L}_{{\rm{X}}}\approx {10}^{33}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ (Pollock 1995). In Figure 3, we show the impact of including X-rays in the PoWR model of the O-type secondary. While the treatment of X-rays is approximate and does not account for nonspherical X-ray irradiation from wind–wind collisions, it illustrates the necessity of X-rays in reproducing the observed spectrum.

Zoom In Zoom Out Reset image size

Finally, the reddening and total luminosity of the system are derived by comparing the observed spectral energy distribution (SED) to the sum of the SEDs of both models. The reddening is modeled following Cardelli et al. (1989). From the total luminosity, the fixed light ratio and the other determined parameters, we then derive the individual luminosities.

The final spectral fit is shown in Figure 4. The inferred parameters are given in Table 4. The spectroscopic mass for the O-type secondary is derived from ${M}_{\mathrm{spec}}=g\,{R}_{* }^{2}/G$, where G is the gravitational constant. The spectroscopic mass of the WR primary stems from mass–luminosity relations calculated for chemically homogeneous WR stars by Gräfener et al. (2011), which for a H-free star depends only on L. The error on the spectroscopic mass of the WR star is propagated from the error in L, and is possibly underestimated in its lower limit given that the WR primary may be chemically inhomogeneous, which would lower M for a given L.

Zoom In Zoom Out Reset image size

Table 4.  Inferred Physical Parameters of WR 133 from the PoWR Spectral Analysis

Parameter	Primary	Secondary
Spectral type	WN5	O9 I
distance (kpc)	1.856 (fixed)
flux $f/{f}_{\mathrm{tot}}(H)$	0.22 (fixed)	0.78 (measured)
${T}_{* }$ (kK)	71 ± 5	30.0 ± 1.0
${T}_{2/3}$ (kK)	67 ± 5	29.5 ± 1.0
$\mathrm{log}g$ (cgs)	−	3.30 ± 0.15
$\mathrm{log}{R}_{{\rm{t}}}$ $({R}_{\odot })$	0.90 ± 0.05a	−
$\mathrm{log}\,L\,({L}_{\odot })$	5.42 ± 0.10	5.29 ± 0.10
${R}_{* }\,({R}_{\odot })$	3.4 ± 0.5	16.6 ± 1.0
D	$4\pm 0.3$ dex	$20\pm 0.3$ dex
β	1 (fixed)	3 ± 1
$\mathrm{log}\dot{M}$ (M⊙ yr−1)	−5.05 ± 0.15	−6.40 ± 0.15
${v}_{\infty }$ (km s−1)	1600 ± 100	1800 ± 200
${X}_{{\rm{H}}}$	${0}_{-0}^{+0.05}$	0.74 ± 0.01
${X}_{{\rm{C}}}/{10}^{-3}$	0.1 ± 0.05	2 ± 1
${X}_{{\rm{N}}}/{10}^{-3}$	15 ± 5	0.7 ± 0.3
${X}_{{\rm{O}}}/{10}^{-3}$	−	6 ± 3
$v\sin i$ (km s−1)	−	96 ± 10
${v}_{\mathrm{mac}}$ (km s−1)	−	40 ± 20
${M}_{\mathrm{spec}}$	14 ± 2	${20}_{-5}^{+8}$
${E}_{B-V}$ (mag)	0.43 ± 0.02

Note.

^aThe transformed radius should not be confused with the physical radius of the star; see the text for details.

Download table as:  ASCIITypeset image

Sana et al. (2012) found that about 75% of massive OB stars should have their evolution impacted by a companion. The general assumption used in the analysis of binary populations is that two stars will interact if the system is born with a semimajor axis of less than ∼10 au and the star evolves into a red supergiant (e.g., Sana 2017). During a binary interaction, the system initially becomes more compact until the mass ratio is unity and then increases until mass transfer eventually ceases. This has been seen in an active phase for binaries that could become similar to WR 133 with systems such as the LBV-like binaries MWC 314 (Richardson et al. 2016b) and HDE 326823 (Richardson et al. 2011). The supergiant companion to the WR star in WR 133 could have a larger rotation rate than an average field star. The inferred value of $v\sin i$ of 96 km s⁻¹ is similar to the average value found for O supergiants by Simón-Díaz & Herrero (2014). However, if there were indeed past binary interactions, the O9I star should be rotating in the orbital plane, and thus would likely be a rapid rotator. In fact, if we use our derived value for $v\sin i$ and the inclination of the system imply a modestly fast rotation rate of 280 km s⁻¹, which is a reasonably fast rotation for a supergiant O star. However, if the O star had a magnetic field from the interaction, then the rotation rate could be drastically less than critical in a short time (Potter et al. 2012).

Eldridge (2009) used the stellar mass and parameters for the O star in the ${\gamma }^{2}$ Velorum system to constrain the age of the ${\gamma }^{2}$ Vel stars using both binary as well as single-star evolution models. As the O-star can be rejuvenated during binary interactions, a comparison with single-star models only provides a lower limit for the age of the system. In the case of ${\gamma }^{2}\,$ Vel, Eldridge (2009) found the system to be about 2 Myr older when taking binary evolution into account. Consequently, we can infer a minimum age of our system when we use our measured O star parameters from Table 3 as an input for the BONNSAI interface (Schneider et al. 2014)¹⁷ to explore any evolutionary constraints that might be possible. Given the input parameters from the orbit, along with the model parameters from Table 4 (e.g., $\mathrm{log}g$), the star was able to be replicated for the single-star BONNSAI models. With that, we find an initial mass of ${27.6}_{-1.6}^{+1.9}{M}_{\odot }$, and an age of ${5.08}_{-0.40}^{+0.47}$ Myr. The age for the ${\gamma }^{2}$ Velorum system was calculated by Eldridge (2009) to be $5.5\pm 1\,\mathrm{Myr}$, in line with the age of the Vela OB association. This means that the WR 133 system is about as old or even older than ${\gamma }^{2}\,$ Vel. Moreover, the mass exchange and envelope stripping of the WR component in WR 133 must have been considerably weaker than in ${\gamma }^{2}\,$ Vel.

It seems that WR 133 is similar to the ${\gamma }^{2}$ Velorum system in that it has had past interactions, and the orbit now has a small eccentricity and a period of ∼100 days. The main difference in these systems is likely the progenitor masses. In the case of WR 133, we see a less massive system than in the case of ${\gamma }^{2}$ Vel, as well as a nitrogen-rich WR star instead of a C-rich WR star. This means that either WR 133 was not massive enough to undergo the same evolution as ${\gamma }^{2}$ Vel, or that it is not as far along in its evolutionary path. In the case of ${\gamma }^{2}$ Vel, we may expect that the larger initial masses would lead to a larger radiative force to drive the stellar winds and evolve to the carbon-rich WR we observe today.

The component masses we derive for WR 133 are lower than those inferred with the PoWR model described in Section 4. We note that the large errors in the masses are currently dominated by the seemingly small error of the orbital inclination of 185, but due to the nearly face-on system geometry, it propagates to large errors, especially when we also consider the errors in the spectroscopic semiamplitudes. While our orbital mass of 22.3 ${M}_{\odot }$ for the O star agrees with the PoWR model results (20 ${M}_{\odot }$), it is noteworthy that Martins et al. (2005) predicts a higher mass of 32 ${M}_{\odot }$ for the O9I spectral type. If the spectral type was O9III instead of O9I, the expected mass from Martins et al. (2005) would be $21{M}_{\odot }$, which is more in line with our findings. We note that many of the spectral classification lines are likely contaminated by the presence of the WR star, potentially blurring the stellar classification of the component stars in the WR 133 system.

Fullard et al. (2020) recently published models of the polarimetric variability of WR 133. With our better-constrained orbital parameters, we can then use the same models for the system to constrain the mass-loss rate of the WR star. Using the parameters from our orbit, such as Ω, the flux ratio, ${a}_{\mathrm{system}}$, and the inclination, the polarimetric model is best fit with a mass-loss rate of the WR star of $(4.0\pm 0.6)\times {10}^{-6}{M}_{\odot }{\mathrm{yr}}^{-1}$, which is about a factor of 2 lower than the PoWR-derived mass-loss rate.

The orbit of WR 133 offers us the chance to better understand the nature of these interactions. The system resides in the Galactic open cluster NGC 6871 (Rosslowe & Crowther 2015). This cluster has a line-of-sight velocity of −11 km s⁻¹ (Kharchenko et al. 2005), comparable to the value we derive for γ_O = −15.3 km s⁻¹. The cluster has an age of ≈10 Myr (Kharchenko et al. 2005), meaning that the system had to form and interact in that relatively short time frame.

Finally, we note that there are three well-established masses for WR stars, where the mass comes from a combination of a visual (interferometric) orbit and a double-lined spectroscopic orbit. The two additional systems are ${\gamma }^{2}$ Velorum (Lamberts et al. 2017), which contains a WC8 star with a mass of 8.9 ${M}_{\odot }$, and WR 140 (Thomas et al. 2021), which contains a WC7 star with a mass of 10.3 ${M}_{\odot }$. It is interesting to note that the two WC stars appear to have similar masses to the WN star in WR 133. However, we continue to be in the range of very small number statistics, and any differences in these masses could be reflecting different conditions for binary interactions in the past.

WR 133 provides a rare chance to constrain stellar wind structure models, in particular for the spectroscopically unconstrained masses of WN stars, which is a key ingredient for our perception of WR-type mass loss (Sander et al. 2020). For these purposes, an improvement of our orbital solution will be vital. This can be done with improved spectroscopic measurements to both constrain the spectral types and better measure the radial velocities. We will need more measurements with the CHARA Array to improve our visual orbit and better constrain the inclination. The two proposed orbits in this paper differ by about 0.02 mas near periastron, which should be possible to measure with the MIRC-X combiner in the future. It also allows us to use polarimetry to provide independent constraints on the mass-loss rate of the WR star, in turn providing for more calibrated measurements for many WR+O binaries in the future.

We thank our anonymous referee for an insightful report that led to better clarity in this paper. We also wish to thank Douglas Gies for clarification about Roche radii estimates. This work is based in part upon observations obtained with the Georgia State University Center for High Angular Resolution Astronomy Array at Mount Wilson Observatory. The CHARA Array is supported by the National Science Foundation under Grant No. AST-1636624 and AST-1715788. Institutional support has been provided from the GSU College of Arts and Sciences and the GSU Office of the Vice President for Research and Economic Development. Some of the time at the CHARA Array was granted through the NOAO community access program (NOAO PropID: 17B-0088; PI: Richardson). MIRC-X received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant No. 639889). This work is also based in part on spectral data retrieved from the ELODIE archive at Observatoire de Haute-Provence (OHP). This research has made use of the Jean-Marie Mariotti Center searchcal and aspro2 service.

Some of L.L.'s participation was funded by the Embry-Riddle Aeronautical University's Undergraduate Research Institute. T.S. acknowledges funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement number 772225). A.A.C.S. is supported by STFC funding under grant No. ST/R000565/1. J.D.M. acknowledges support from NASA-XRP NNX16AD43G and NSF-AST 1909165. S.K. acknowledges support from ERC Starting Grant "ImagePlanetFormDiscs" (grant Agreement Number 639889).

Facilities: CHARA - Center for High Angular Resolution Astronomy interferometer, DAO:1.83m - .