© V. Ramachandran et al. 2022

1. Introduction

The majority of massive stars are commonly found in binary or higher-order multiple systems (e.g., Sana et al. 2012, 2014). The evolution of these massive stars is largely governed by the loss of mass via stellar winds, and by the exchange of mass among binary components (e.g., Langer 2012; Postnov & Yungelson 2006). When the initially more massive (i.e. primary) component in a binary system collapses into either a neutron star (NS) or a black hole (BH) at the end of its evolution, some systems remain bound after the supernova explosion. The accretion of mass transferred from the former secondary star onto its compact companion can power strong X-ray emission (Shklovsky 1967). Such systems are called high-mass X-ray binaries (HMXBs). The majority of HMXBs contain NSs, while systems containing BHs are quite rare. In this paper, we focus on M33 X-7, an HMXB with a BH companion.

Black hole X-ray binaries (BH XRBs) are essential natural laboratories for studying stellar-mass BHs. These systems provide an indispensable observational test for BHs that are formed by core collapse, therefore providing observational constraints for binary evolution and compact object formation. The majority of the detected BH XRBs have low-mass donor stars and are X-ray transients (Remillard & McClintock 2006; Özel et al. 2010; Corral-Santana et al. 2016). BH XRBs with massive donor stars are persistently strong X-ray sources that are powered by massive stellar winds and/or Roche lobe overflow. Cyg X-1 was the first discovered BH HMXB system (BH + O9I companion) and it remains the only confirmed one in the Milky Way (e.g., Paczynski 1974; Gies & Bolton 1986; Miller-Jones et al. 2021). Other BH HMXB candidates in the Galaxy are SS433 (Cherepashchuk et al. 2020), Cyg X-3 (Zdziarski et al. 2013), and MWC656 (Casares et al. 2014). In the extra-galactic population of BH HMXBs, dynamical evidence has been presented for LMC X-1 (BH + O7 III companion; Orosz et al. 2009), LMC X-3 (BH + B3–5, V star Orosz et al. 2014), and M33 X-7 (BH with an O7-8 III companion; Orosz et al. 2007). There is indirect evidence for the presence of a BH with a Wolf-Rayet companion in NGC300 X-1 (Binder et al. 2021; Crowther et al. 2010) and IC10 X-1 (Laycock et al. 2015; Silverman & Filippenko 2008) systems.

If the donor star in a BH HMXB is massive enough, its core collapse can result in the formation of a relativistic binary comprising two BHs. In the case of a tight system, the binary components may merge within a Hubble time. Such a merger may be observed as a gravitational wave (GW) event (e.g., Abbott et al. 2016a,b, 2017). Thus, HMXBs with BHs represent a key transitional stage toward a binary BH merger. It is crucial to understand the physics of these systems and their dependence on metallicity in order to establish whether the BHs in present-day massive binaries and those detected by GWs during a merger are part of the same stellar progenitor pool.

In HMXBs, accretion onto the compact object is typically powered by the stellar wind of the donor star. For Bondi-Hoyle-Lyttleton accretion (Bondi & Hoyle 1944; Davidson & Ostriker 1973; Edgar 2004), the X-ray luminosity crucially depends on the wind parameters of the donor star. In particular, Bondi-Hoyle mass accretion is highly dependent on the velocity at which the donor wind is captured by the gravitational field of the compact object (v⁻³, from Bondi & Hoyle 1944). The X-ray luminosity is directly dependent on the accretion rate and accretion mode and can vary rapidly due to the complexity of the stellar wind and the possible clumping of material. Hence it is of the utmost importance to understand the wind properties of the massive donor. Sophisticated spectral analyses of the donor stars using state-of-the-art model atmospheres to determine the stellar and wind parameters of the donor stars in HMXBs are therefore a pivotal instrument for a better understanding of these systems.

The winds of OB stars are driven by the scattering of UV photons in metal lines (Castor et al. 1975). At lower metallicity, the winds of the OB stars are much weaker (Bouret et al. 2015; Ramachandran et al. 2019) and hence remove less mass and angular momentum during stellar evolution. Thus, a metal-poor star is expected to undergo its core collapse at a higher mass and at faster rotation than a star with initially the same mass, but higher metallicity. As of yet, virtually nothing is known about the winds of metal-poor donor stars in HMXBs with BH companions. Our current understanding of the wind of massive stars in BH binaries is mostly based on the Cygnus X-1 system in our Milky Way. Yet, this prototypical system has not been analyzed with the current generation of expanding stellar atmosphere codes. Spectroscopic analyses of the BH HMXBs were limited to the Sobolev with Exact Integration (SEI) method (e.g., Cyg X-1 in Gies et al. 2008; Vrtilek et al. 2008) or using plane-parallel models (e.g., for M33 X-7 by Orosz et al. 2007 and for Cyg X-1 by Caballero-Nieves et al. 2009) which are optimized for hot stars with no significant wind. More advanced analyses of donors done in the recent past are limited to NS HMXB systems. For example, a sample of HMXBs with NS companions was recently studied by Hainich et al. (2020) with detailed expanding stellar atmosphere models for radiation-driven winds. In this study, we carry out a detailed atmospheric analysis of multiwavelength spectra of the metal-poor donor star in M33 X-7, taking into account line-driven winds, wind clumping, complex effects of millions of spectral lines, and X-ray photoionization.

The M33 X-7 system. M33 X-7 is reported to contain one of the most massive known stellar mass BHs with a very massive O7 III star in a tight binary orbit (Orosz et al. 2007). The system is located in the metal-poor galaxy of Messier 33 (M33). M33 is known to have a steep metallicity gradient (Cioni 2009; Magrini et al. 2007). The M33 X-7 system is located in an OB association (HS 13) approximately 2 kpc from the center of M 33. Based on H II region abundance measurements by Magrini et al. (2007), the metallicity at this galactic radius is roughly half solar.

Importantly, M33 X-7 is the only known eclipsing BH HMXB (Pietsch et al. 2004, 2006). Pietsch et al. (2006) derived a period of 3.45 days and an X-ray luminosity of 5 − 11 × 10³⁷ erg s⁻¹ out of eclipse. The study conducted by Orosz et al. (2007) provided a dynamical estimate of the mass of the system, reporting a ≈70 M_⊙ O supergiant donor, and an ≈16 M_⊙ BH. The mass of the donor is very high for an O7 III spectral type (Markova et al. 2018; Ramachandran et al. 2018) as well as compared to O stars in other BH HMXBs. Such high masses for the donor and compact companion are exceptional among HMXBs, making M33X-7 a unique laboratory for studying the potential progenitors of massive BH mergers detected by LIGO.

Despite a rather tight orbit, it was reported that the O-type donor does not fill its Roche lobe radius (R_* ≈ 0.77R_rl, Orosz et al. 2007), and therefore the X-ray luminosity is expected to be powered by accretion of the donor’s stellar wind alone. Orosz et al. (2007) estimated the donor mass-loss rate from the X-ray luminosity by adopting an ad hoc accretion efficiency. The estimated mass-loss rate turned out to be five times higher than theoretically expected. Furthermore, an ad hoc wind velocity law (and consequently the radial density structure of the wind) was adopted for the donor star. On the basis of these assumptions, the eclipse width was measured from the light curve with a suggested accuracy of 1 degree in terms of orbital phase angle. As a result, it was claimed that the component masses are estimated with 10% accuracy.

M33 X-7 has been subject to evolutionary modeling in past studies (Abubekerov et al. 2009; Valsecchi et al. 2010; de Mink et al. 2010; Bogomazov 2014). A large theoretical effort is ongoing to explain the large masses of the M33X-7 components, its tight orbit, and the large BH spin. Measuring the stellar wind parameters of M33 X-7 will significantly constrain evolutionary models and allow unique tests of the evolution of BH binaries at low metallicity.

For this study, we have secured quasi-simultaneous multiwavelength spectroscopy of M33 X-7 at three key binary orbital phases (see Fig. 1). Using the combined analysis of the multiwavelength data (UV, optical, and X-rays) we measured the properties of the donor star wind and the mass accretion rate onto the BH, and we study the effect of X-rays on the donor star wind. These wind parameter measurements are required to understand how much wind is quenched by the X-rays, so as to test stellar wind theory, refine the measurements of the BH mass, and determine the accretion efficiency onto the BH. Subsequently, our new realistic wind parameters serve as input for sophisticated massive binary evolution models. The basic parameters of the system which are used in this study from the abovementioned literature are listed in Table 1. Throughout the paper, we assume a distance of 840 kpc to M33 (Galleti et al. 2004; Gieren et al. 2013).

Fig. 1. Sketch of the binary orbit in M33 X-7. The observer is located along the horizontal axis on the right. The thick lines span the range of orbital phases covered by the HST observations. The orbital separation is taken as ≈1.75 stellar radius as obtained in this work (see Table 4).

This paper is organized as follows: We discuss the new and archival multiwavelength spectroscopic observations of M33 X-7 in Sect. 2. Based on the phase-resolved UV spectroscopy, we report orbital variations in the wind line profiles in Sect. 3. In Sect. 4, we describe the detailed analysis of the spectra. In Sect. 5 we present the updated stellar parameters of the system. We address the detailed modeling of the donor star wind in Sect. 6 and estimate the accretion rate in Sect. 7. Section 8 discusses the evolution and formation of the system using binary evolution models. Section 9 contains our final conclusions about this BH HMXB system.

2. Observations

We carried out simultaneous HST and XMM-Newton observations of M33 X-7. These includes observations at two different phases of the orbit: quadrature (ϕ ≈ 0.33) and inferior conjunction (ϕ ≈ 0.53). These two observations were taken four days apart in July 2019. In addition we acquired HST spectra during X-ray eclipse (ϕ ≈ 0.95) of the system in September 2020. From the latter observations we want to obtain stellar and wind parameters of the O-star largely undisturbed by the BH. A schematic illustration of orbit with the observed phase ranges is shown in Fig. 1. The details of the spectroscopic observations are summarized in Table 2.

The planning of the joint HST and XMM-Newton campaign relied on the X-ray light-curve and ephemeris which were based on observations with the Chandra X-ray telescope in 2005 (Pietsch et al. 2006). According to Pietsch et al. (2006), the mid-eclipse time is occurring at HJD(2453639.119 ± 0.005)±N(3.453014 ± 0.000020). However, the BH binaries are dynamic systems, therefore the ephemeris might not be stable and the BH might change states.

To check the validity of ephemeris we obtained monitoring observations of M33 X-7 with the Neil Gehrels Swift Observatory. Swift visited M33 X-7 seven times between 2019-07-02 and 2019-07-08. During five visits, the count rates were consistent with that of the Chandra’s ACIS-I detections in 2005. On the other hand, on HJD 2458666.828 and HJD 2458670.147 the system was not detected in X-rays which we interpret as X-ray eclipses. Using the ephemeris in Pietsch et al. (2006), these correspond to N = 1456.035 and N = 1456.996 cycles respectively, that is to say they are in agreement within uncertainties (the eclipse duration is ≈0.2 in phase). Thus we verified that the orbital parameters of the system are still valid and hence were applied for scheduling our observations at three key orbital phases.

2.1. UV spectroscopy using HST/COS

We secured a total of 12 HST orbits in cycle 26 under the GO program 15636 (PI: Oskinova). Each of the three selected phases (see Fig. 1) was observed during four orbits. We selected the COS FUV G140L grating with a central wavelength setting of 1280 to cover the entire 912–2150 Å wavelength range (Green et al. 2012). This configuration includes two spectral segments A and B with a gap in the range of 1192–1265 Å, omitting the Lyα. We disregard wavelengths λ < 1100 Å and > 1800 Å due to low signal-to-noise ratio. The resulting UV spectra have a spectral resolving power of 1500–2000. The observations were carried out in TIME-TAG mode.

Crowding in the COS aperture. Since M33 X-7 resides in an OB association (Ma et al. 2001; Pietsch et al. 2006), separating the contamination of the COS UV spectrum by nearby sources is challenging. The COS circular aperture is 1.25 arcsec in radius, but the aberrated beam entering the aperture allows objects up to 2 arcsec from the center to contribute to the recorded target spectrum. The most noticeable contribution comes from a bright neighboring cluster of stars, which is separated by only 0.9″ to the southwest (see Fig. 2). The position angle of the instrument during observation further reduced the distance along the cross-dispersion axis to 0.8″.

Fig. 2. HST/COS aperture around M33 X-7. The background image is an RGB composite of HST images in F336W, F475W, and F814W bands. A nearby cluster falls within the FOV of the COS spectrograph.

We extracted photometry of M33 X-7 and the cluster from high-resolution HST images in F275W, F336W, F475W, and F814W filters. The fluxes were measured using 0.2″ aperture photometry and applied to the appropriate zero point. At 814 nm, the brightness of the nearby cluster is just 0.2 mag lower than M33,X-7, however at 275 nm, the cluster becomes 0.6 mag fainter than our target. In the FUV, only GALEX photometry of the region is available, with a spatial resolution of ≈6″. Even though M33 X-7 dominates in the FUV, we expect a non-negligible contribution from the nearby cluster. This contamination could be larger in the optical and thus affect the parameters obtained by Orosz et al. (2007).

In order to separate contributions from nearby sources in the COS aperture, we optimally extracted the one-dimensional target spectrum from the two-dimensional flat-fielded detector images. For this purpose we used observations carried out at ϕ = 0.33 and 0.53. Observation at these two phases were carried out with a position angle of ≈116°, implying a separation of 0.8″ between our target and the nearby cluster along the cross-dispersion axis. At ϕ = 0.95 the position angle is ≈96°, so the corresponding separation is only ≲0.6″. In the two-dimensional spectra we could see a double peak profile on the cross-dispersion axis indicating contributions from both M33 X-7 and neighboring cluster. Hence we simultaneously fit for both our target and the cluster by multiplying them with the (wavelength-dependent) cross-dispersion instrumental profile and subsequently subtract from the combined observed spectra. This has been done along the cross-dispersion axis for all wavelengths. Finally, we extracted two spectral components such that the residual is minimum. We found that the disentangled spectra of M33 X-7 have the same spectral features as the pipeline extracted spectra. Based on all exposures taken at ϕ = 0.33 and 0.53, we found that on average M33 X-7 accounts only for ≈75% of total flux, whereas the nearby cluster contributes the remaining 25%.

2.2. X-ray spectroscopy using XMM-Newton

We carried out dedicated X-ray observations of M33 X-7 with the XMM-Newton telescope with a goal of obtaining X-ray spectroscopy as close in time as possible to the HST observations in quadrature and inferior conjunction. The log ofthe XMM-Newton observations is given in Table 2.

The XMM-Newton data were analyzed using the Science Analysis System (SAS)¹. The standard analysis steps were performed to extract the spectra measured by the European Photon Imaging Camera (EPIC). EPIC consists of a PN (Strüder et al. 2001) and two MOS (Turner et al. 2001) CCD detectors in the focal plane, which operate in the 0.2–12 keV energy range. During our observations, all cameras were operating in the standard, full-frame mode with a thin filter. To analyze the EPIC spectra we used the standard spectral fitting software XSPEC (Arnaud 1996).

The extracted spectra and the best fit spectral models are shown in Figs. 3 and 4, while the best fit model parameters are listed in Table 3. It is interesting to note that unlike HMXBs with neutron star companions accreting from the stellar wind of OB-type supergiants (Giménez-García et al. 2015; Sanjurjo-Ferrín et al. 2021), no Fe Kα line is seen in the XMM-Newton EPIC spectra of M33 X-7. After testing various models often used for describing the spectra of accreting BHs, the simplest model which well describes the spectra is a combination of thermal plasma apec and a multiple blackbody disk model diskpbb. The latter model is characterized by the inner disk temperature T_in and the scaling of the local disk temperature depends on the radius, T ∝ r^−p, with p being a free parameter (Kubota et al. 2005).

Fig. 3. XMM-Newton EPIC PN, MOS1, and MOS2 spectra of M33 X-7 (black, red, and green curves, respectively) obtained at phase ϕ ≈ 0.47. The error bars correspond to 3σ; the best-fit combined thermal plasma (apec) and diskpbb model is shown by solid lines. The model parameters are shown in Table 3. Lower panel: residuals between the data and the best-fit model.

Fig. 4. Same as in Fig. 3, but for phase ϕ ≈ 0.33.

In the observation at phase ϕ ≈ 0.47, i.e. close to the inferior conjunction when the observer has the direct view on the accreting BH, the emission measure of the collisionally ionized plasma (described by the apec model) appears to be higher than at ϕ ≈ 0.33, i.e. close to quadrature. We interpret this apec plasma component as originating in the accretion disk wind. Furthermore, close to the inferior conjunction, the parameter norm which scales with the accretion disk inclination as cos θ (with θ = 0 being face-on) is twice as high as quadrature. This might imply that we see a larger fraction of the disk at inferior conjunction. Interestingly, the X-ray absorption, as characterized by the parameter N_H, does not change between the two orbital phases. This implies an insignificant contribution to the X-ray absorption from the general wind of the O-star beyond the BH orbit which could be due to the wind quenching by strong X-ray irradiation (see Sect. 6.1).

2.3. Archival optical spectra

We made use of the optical spectra of M33 X-7 published by Orosz et al. (2007) obtained from August 18 to November 16, 2006, with the 8.2-m Gemini North Telescope at wavelengths λ = 4000 − 5000 Å. The final spectra is an average of spectra taken over 22 nights, and used these data for our spectral analysis. Although the optical data is not phase-resolved like UV or X-ray, it is sufficient for our spectral analysis since we do not expect any notable variation in the absorption lines at the blue optical range.

3. Stellar atmosphere modeling

We carefully analyzed the UV spectra obtained at three different phases (eclipse, quadrature, and inferior conjunction). As a first approximation, we used the stellar parameters obtained by Orosz et al. (2007), where they used plane-parallel models from the TLUSTY grid to analyze the spectra. This approximation is only suitable for stars without winds where the spectrum is purely formed below the sonic point in a quasi-hydrostatic photosphere. However, since the donor in M33 X-7 is an O giant or supergiant with a dense wind, the spectrum is altered in the layers beyond. This can cause an impact on the derived stellar parameters such as log g. So-called unified model atmospheres are necessary to consistently describe the outermost layers of the star and their winds (e.g., Gabler et al. 1989; Sander et al. 2017). Finally, X-rays need to be accounted for the spectral analysis of the HMXB donor since they can have a noticeable effect on the ionization structure of the wind. In turn this affects the derived parameters of the donor and its wind.

To analyze the UV and optical spectra, we use the Potsdam Wolf–Rayet (PoWR) model², which is a state-of-the-art non-local thermodynamic equilibrium (NLTE) code. The PoWR code assumes a spherically symmetric outflow and accounts for iron-line blanketing, wind inhomogeneities, a consistent stratification in the quasi-hydrostatic part, and irradiation by X-rays, all of which are necessary to measure stellar and wind parameters. To achieve a consistent solution, the equations of statistical equilibrium and radiative transfer are iteratively solved to yield the population numbers, while accounting for energy conservation. The radiative transfer is solved in the comoving frame, which avoids simplifications such as the Sobolev approximation. Once an atmosphere model is converged, the synthetic spectrum is calculated via a formal integration along emerging rays. Details of the PoWR code are described in Gräfener et al. (2002), Hamann & Gräfener (2003), Todt et al. (2015), and Sander et al. (2015). However, it should be noted that the standard version of the PoWR code does not account for deviations from spherical geometry.

A PoWR model is specified by the stellar temperature T_*, the bolometric luminosity L, the surface gravity g_*, the mass-loss rate Ṁ, the wind terminal velocity v_∞, the velocity law, and the chemical abundances. The stellar radius R_* (inner boundary) is defined at a Rosseland continuum optical depth of τ_Ross = 20 from T_* and L. Following the Stefan-Boltzmann law $L=4\pi {\sigma }_{\text{SB}}{R}_{\ast }^2{T}_{\ast }^4$ the stellar temperature T_* is the effective temperature corresponding to the stellar radius R_*. The outer boundary in our models is set to R_max = 100R_*.

In the main iteration, thermal broadening and turbulence are approximately accounted for by assuming Gaussian line profiles with a Doppler width of 40 km s⁻¹. In the formal integral for the calculation of the emergent spectrum, the Doppler velocity is split into the depth-dependent thermal velocity and a microturbulence velocity ξ(r). A microturbulent velocity of 10 km s⁻¹ was adopted in the photosphere, growing proportional with the wind velocity.

The velocity field in PoWR models consists of two parts. In the inner part of the stellar atmosphere, the velocity field is calculated consistently such that the quasi-hydrostatic density stratification is fulfilled (Sander et al. 2015). In the supersonic region, the wind velocity field v(r) is prescribed assuming a so-called β-law v(r) = v_∞(1−R_*/r)^β (Castor et al. 1975), where v_∞ is the terminal wind velocity.

In PoWR models, wind inhomogeneities are accounted for in the microclumping approach that assumes optically thin clumps (Hillier 1991; Hamann & Koesterke 1998) with a void interclump medium and described by a so-called density contrast D(r). The matter density in the clumps is enhanced by a factor D = 1/f_V, where f_V is the fraction of volume filled by clumps. In the current study, we account for a depth-dependent clumping factor.

The models account for complex atomic data of H, He, C, N, O, Mg, Si, P, S, and the iron-group elements. Iron-group elements with millions of lines are included in the PoWR code through the superlevel approach (Gräfener et al. 2002).

The PoWR code can account for ionization due to X-rays. The X-ray emission is modeled as described by Baum et al. (1992), assuming that the only contribution to the X-ray flux is coming from free–free transitions. Since the current generation of PoWR models is limited to spherical symmetry, the X-rays are assumed to arise from an optically-thin spherical shell around the star. The X-ray emission is specified by three free parameters, which are the fiducial temperature of the X-ray emitting plasma T_X, the onset radius of the X-ray emission R₀ (R₀ > R_*), and a filling factor X_fill, describing the ratio of shocked to nonshocked plasma. We set the onset radius to the orbital distance between the donor star and the BH. To obtain the fiducial temperature which describes the spectral shape in the X-ray regime, we fit the X-ray spectra with the brems model. The resulting temperatures are shown in Table 3. As the next step, the X-ray filling factor is adjusted such that the wavelength integrated X-ray flux from the observations are reproduced by the model. This is a similar procedure as adopted in Hainich et al. (2020).

4. Stellar parameters

We primarily used UV spectra taken at the time of eclipse to measure accurate stellar parameters of the donor star since in this phase we see the part of the star which is largely undisturbed by the BH. We also re-analyzed the optical spectra from Orosz et al. (2007, see Figs. A.1–A.3). We started the model calculations with the parameters provided by Orosz et al. (2007) and adjusted them to get the best fit. The best-fitting model is selected after several iterations of visual inspection and systematic variation of the parameters. We simultaneously reproduce the normalized line spectrum and the spectral energy distribution (SED) of M33 X-7 using our best-fit PoWR model, as shown in Figs. 5, 6, and A.1. Details on the wind parameters are described in Sect. 6.

Fig. 5. Spectral energy distribution of M33 X-7. The model SED (red solid line) is adjusted to fit the photometry in UV and optical bands (blue boxes).

Fig. 6. Spectral fit for M33 X-7 UV spectra at the X-ray eclipse. Normalized HST spectra (blue solid line) is overplotted with the PoWR model (red dashed line). The atmospheric parameters and abundances of this best-fit model are given in Tables 4.

The observed spectra are best reproduced by a model with T_* = 31 kK and log g_* = 3.4, which is in line with the spectral class of a late O giant. To constrain the temperature, we used photospheric diagnostic lines such as C III and Si III in the UV along with He I and He II lines in the optical spectra. Our estimated values of T_* and log g_* are lower by 4 kK and 0.35 dex, respectively, than previously suggested by Orosz et al. (2007). How the variations in T_* and log g_* affect different line profiles is depicted in Figs. A.1–A.3). Hotter models would result in stronger He II lines in the optical and weaker C IV and Si IV lines in the UV. We measured the surface gravity using the pressure-broadened wings of the Balmer lines. Since the ionization balance also reacts on log g_*, we readjusted the temperature accordingly. Overall, our best-fit model agrees well with the observed UV and optical lines. Subsequently, the emergent spectrum of the PoWR model was convolved with a rotational velocity of 250 km s⁻¹ and macro-turbulent velocity of 150 km s⁻¹ to reproduce the shape of the observed line profiles. Taking into account the updated inclination of the system, the star rotates with v_rot/v_critical ≈ 0.55.

The luminosity, color excess E_B − V, and the extinction-law parameter R_V were determined by fitting the model SED to the photometry and flux calibrated UV spectra (see Fig. 5). We have used UBV broad-band magnitudes from Pietsch et al. (2006) and R magnitude from Massey et al. (2006). The extinction A_V = R_V × E_B − V derived by the final SED fit is in agreement with that of Orosz et al. (2007). Subsequently, the model flux is scaled with the adopted distance modulus of 24.62 mag (see Orosz et al. 2007, for more details). The foreground Galactic extinction is set to 0.052 mag (according to the NED extinction calculator, Schlegel et al. 1998), and the intrinsic extinction in M33 is negligible according to the detailed dust maps by Hippelein et al. (2003). The derived luminosity was found to be lower than obtained by Orosz et al. (2007).

The Iron forest present in the FUV is most useful to constrain the metal content of the stars. Previous studies assumed a low metallicity of Z = 0.1 Z_⊙ for the system. However, we find that the iron forest in the UV range is best reproduced with models assuming an LMC-like metallicity of Z = 0.5 Z_⊙ (see Fig. 7). Considering the galactocentric distance of ≈2 kpc from the nucleus of M33, this is in agreement with ISM abundances reported in the M33 star-forming regions (Magrini et al. 2007; Alexeeva & Zhao 2022).

Fig. 7. Iron forest in the FUV range observed in the M33 X-7 spectra (black line). Over plotted in blue, red, and green are synthetic spectra obtained for models with the same parameters but with different iron abundances, ranging to 1, 1/2, and 1/5 Fe⊙ respectively.

We adjusted the abundances of individual elements to reproduce the observed strength of their respective absorption lines. The α-element abundances are also in agreement with LMC-like metallicity. We found evidence for nitrogen enrichment and carbon deficiency in the observed spectra. We needed to reduce the hydrogen mass fraction to 60% and enrich helium up to 40% to account for the strength of the He I lines in the optical spectra. The He enrichment could be a result of the fast rotation of the donor. Alternatively, the outer layers might have been stripped by interaction with the BH companion. The UV spectra also show indications for CNO processing in the star. In order to account for the line strengths, we had to reduce the carbon abundance and increase the nitrogen abundance by a factor of five from baseline values. Individual element abundances and parameters derived from the spectral fitting are tabulated in Table 4. For elements Mg, P, S, Al, and Ne, we adopted typical LMC abundance values derived from OB stars (Trundle et al. 2007).

The masses of the components in M33 X-7. The reduction in luminosity, temperature, and gravity resulted in a much lower spectroscopic mass of ≈38 M_⊙ for the O star than previously found in Orosz et al. (2007). The spectroscopic mass is calculated from log g_* and R_* ($g_{\ast }=G{M}_{\ast }{R}_{\ast }^{-2}$), and the uncertainties in these parameters propagates into the final uncertainty of the donor mass. The newly estimated mass falls in a similar range as that of O star donors in other BH HMXBs such as Cyg X-1, and LMC X-1 (Miller-Jones et al. 2021; Orosz et al. 2009).

Considering the orbital parameters derived by Orosz et al. (2007), we estimated the BH mass, the orbital radius and the semi-eclipse angle for different values of the inclination angle. We found that an inclination angle of 65.5° best matches the observed semi-eclipse angle of 26.5° (Pietsch et al. 2006). If we take into account the effects of an extended wind from the O-star, the donor radius slightly increases, therefore reduce the inclination angle to 65°. In comparison, the inclination angle derived in Orosz et al. (2007) was 74.6°. The newly derived BH mass is ≈11 M_⊙, which is much lower than the previous estimate of ≈15.65 M_⊙. The scaling between the estimated masses of the donor and the BH along with their respective uncertainties are illustrated in Fig. 8.

Fig. 8. Mass plane diagram for components in the M33 X-7. Masses derived in this work are marked in red whereas gray points represent previously determined values from Orosz et al. (2007). The uncertainty in masses is shown as error bars. The solid lines represent loci of constant orbital inclination angle of 65° derived in this work. The dashed line corresponds to the inclination angle of 75° derived by Orosz et al. (2007).

The mass ratio of the system is now ≈3.3. Based on the revised masses, we recalculated the orbital separation to 35.3 R_⊙, whereas the radius of the donor is ≈20.5 R_⊙. The distance to the BH is always the same as the orbital separation, since the system has a nearly circular orbit. According to our results the BH is located much closer to the O giant than suggested by Orosz et al. (2007). Abundances and stellar parameters of the O star indicate that it is close to the main-sequence turn-off, overfilling its Roche-lobe radius.

5. Orbital modulations in the UV wind lines

In HMXBs, the X-ray source produces an extended zone of high ionization that moves around the orbit. When this region is in front of the donor (ϕ = 0.5), the stellar wind gets over-ionized by X-ray irradiation. Therefore the intensity of the UV resonance lines is reduced due to the depletion of the absorbing ion species, while the opposite occurs at X-ray eclipse. These orbital variations were first suggested by Hatchett & McCray (1977) and they are expected to be observable in the ultraviolet resonance lines of ions such as C IV, Si IV, or N V.

In Fig. 9 we present the UV spectra of M33 X-7 taken at three different orbital phases. We identified prominent P Cygni profiles in the resonance lines of C III at 1176 Å, Si IVλλ1393-1403 and C IVλλ1548-1551 at X-ray eclipse. By comparing UV and optical spectral features to that of normal O stars, we re-classify the donor to an O9 II spectral type. We notice a significant change in these P Cygni profiles at phases when the compact object is in the line of sight or in egress compared to that during eclipse. Such phase-dependent profiles variation of UV resonance lines across different orbital phases can be attributed to the X-ray ionization of the wind.

Fig. 9. Comparison of HST/COS spectra of M33 X-7 taken during different orbital phases. Line colors indicate different orbital phases as given in the legend. Upper panel: the full spectra range and lower panels focus on wind line profiles. The observed line variations are due to H–M effect (see Sect. 5).

The Hatchett-McCray effect (HM effect) is not observed for all HMXB systems since the orbital configuration influences the appearance of the wind line profiles. We report a strong HM effect in M33 X-7 (compared to Galactic HMXBs), in both Si IV and C IV lines, indicating that the donor wind is significantly disturbed by the accreting BH. During the eclipse, the donor star wind profiles are much stronger and broader whereas the other two phases display a large reduction in the extent and depth of the blue-shifted absorption component. Although variations at Si IV and C IV are most prominent at high velocities (−500 to −1500 km s⁻¹), we found notable changes also at low velocities (−100 to −500 km s⁻¹). The variation between phase 0.33 and phase 0.53 is mostly negligible, except at lower velocities, suggesting that the wind ionization extends all the way toward the photosphere of the O-star when the BH is in the line of sight, leaving only a “shadow wind” (see e.g., Blondin 1994) unaffected.

The redward emission component of the P Cygni lines is significantly reduced or absent in the eclipse spectra compared to that found in typical O giants. We find evidence for an increase in the emission component during phase 0.33 and phase 0.53. The He II at 1640 Å emission is also found to be stronger when the BH is in the line of sight. The reduction in the strength of the red emission component in the UV resonance lines is possibly due to the loss of ions by strong X-ray photoionization, because the red emission component forms in the opposite part of the star facing the BH.

6. Modeling the donor wind

Winds of massive stars are characterized mainly by their mass-loss rate Ṁ, terminal wind velocity v_∞, velocity law exponent β, and clumping factor D. In this work, we used UV spectra taken at ϕ = 0.95 to yield parameters of the O star wind and spectra taken at ϕ = 0.53 and 0.33 to understand the wind quenching by X-rays.

The main diagnostic lines used to determine the terminal velocity are Si IVλλ1393–1403, and C IVλλ1548–1551. The maximum blue edge velocity in the unsaturated P Cygni lines is v_edge ≈ 1900 km s⁻¹. We find a narrow absorption feature (narrow absorption component, NAC) at a blueshifted velocity of ≈1500 km s⁻¹. Based on this, we adopted the terminal velocity to be 1500±100 km s⁻¹, which is also in agreement with v_∞ ≈ 0.8 v_edge (Howarth & Prinja 1989; Prinja et al. 1990). The ratio of terminal velocity to the escape velocity of the donor is found to be v_∞/v_esc ≈ 2, lower than for typical O stars (Lamers et al. 1995; Kudritzki & Puls 2000).

A microturbulent velocity of 10 km s⁻¹ in the photosphere is derived based on the optical lines, and is also in agreement with the UV spectra. The turbulent velocity grows proportional to the terminal velocity in the wind regime. In the outer part of the wind, it reaches typically up to a value of ξ(R_max) = 0.1v_∞. However, for M33 X-7 we measure a higher turbulence velocity, with a ratio ξ(R_max)/v_∞≈ 0.3. The wind turbulent velocities derived here are larger compared to the previous studies of OB stars in the LMC and the Galaxy (e.g., 0.1 v_∞ by Kudritzki & Puls 2000, 0.14 v_∞ by Herrero et al. 2001). A UV spectroscopy study of early B supergiants in M33 by Urbaneja et al. (2002) also reported high turbulent velocities in the range of 0.25–0.35v_∞.

The adopted microclumping is depth-dependent, assumming that clumping begins at the sonic point, and increases outward. The maximum value of the clumping factor is found to be D = 40, which is required to reproduce the UV lines, especially P Vλλ1118–1128 and N IV at 1718 Å. Our analysis suggested that the clumping factor reaches the maximum value only beyond 3R_*. Compared to this distance, the location of the BH is much closer (1.75R_*) to the photosphere of the donor star. Near the location of the BH, the wind clumping is lower (D = 10).

One noticeable feature in the UV spectra are the weak unsaturated P-Cygni lines with negligible or no emission component. Whereas, OB supergiants with similar spectral types typically show Si IV and C IV resonance lines with saturated blue absorption and strong emission in the red part of the line. We tried different approaches such as adjusting the mass-loss rate, the velocity law, the outer boundary of the wind, and including X-rays, to avoid the red emission part and to improve the fit quality to the observed UV lines. For that, we computed models with log Ṁ in the range of −6.0 to −9.0. To model the supersonic part of the wind, we tested β values in the range 0.6 to 3.0. A comparison of different possible solutions that can reproduce the observed spectra and our final best-fit parameters are discussed below.

As a first attempt, we tried to change the wind parameters in the model without incorporating X-rays. A model with a low mass-loss rate (log Ṁ ≈ −7.2) assuming a shallow wind velocity law (β = 2) is found to be in good agreement with the observations at phase 0.95 (orange line in Fig. 10). The observed profile can also be reproduced by models with a moderate mass-loss rate (≈ − ​6.6) and adopting a velocity law exponent of β = 0.8 if we assume the wind to be shut off by X-rays beyond 900 km s⁻¹. On the other hand, the theoretically predicted mass-loss rate from Vink et al. (2000, 2001) for the derived donor parameters is much higher, log Ṁ ≈ −6.2. Assuming such a higher Ṁ and β = 0.8 in the model would result in much stronger P Cygni profiles compared to the HST spectra taken during the X-ray eclipse (see the black line in Fig. 10). Models with a lower mass-loss rate would make the Si IV P Cygni lines much weaker, but not have much impact on the C IV profile.

Fig. 10. Comparison of observed Si IV (left panel) and C IV (right panel) line profiles (blue solid) at eclipse to the models calculated with different wind parameters and X-rays. Various colors represent changes in the adopted model parameters as given in the legend. Please see text for more details.

Compared to eclipse, the observed spectra taken at ϕ = 0.53 and 0.33 do not show any wind line features. Both Si IV and C IV appear to be mostly photospheric (see Fig. 11). To mimic the observed spectral features at this phase without accounting for the effect of X-ray photoionization, we need to adopt a very low mass-loss rate of log Ṁ ≲ −8.6.

Fig. 11. Same as Fig. 10, but at inferior conjunction

6.1. Impact of X-rays on the UV spectra

X-ray irradiation from the BH accretion disk makes a significant impact on the ionization stratification of the donor wind. While the photospheric absorption lines are unaffected, the appearance of spectral lines that originate in the wind, most notably in the UV range, changes significantly. Even adding a small amount of X-rays (L_x ≲ 10³² erg s⁻¹) in the model to account for wind intrinsic X-ray emission of the star, changes the ionization structure of the donor (e.g., Oskinova et al. 2011). However, the line profiles can get either stronger or weaker depending on the wind density near the BH and the X-ray intensity. To understand the impact of X-ray photoionization on the observed UV spectra, we have calculated models with different X-ray luminosities (10^{32 − 38} erg s⁻¹). In our spherically symmetric models, the X-rays are assumed to be generated from an optically-thin spherical shell around the star, whereas the geometry of the HMXB is more complex. In this system, the BH is behind the donor during the eclipse and in the foreground during other phases. On top of this, the radius of the donor star is large enough to fill its Roche-lobe (see also Sect. 7). Nevertheless, this is one of the first attempts to study the detailed wind driving of a massive star in a BH HMXB.

The effects of the X-ray field on the emergent spectra are illustrated in Fig. 12. The mass-loss rate is fixed to log Ṁ ≈ −6.2 for all cases. The synthetic spectra generated by models with no X-rays and moderate X-ray irradiation (L_x ≲ 10³⁴ erg s⁻¹) show strong saturated Si IV and C IV P Cygni profiles. Including X-rays makes a pronounced effect on wind lines of highly ionized species like O VI and N V. These lines are absent in models which do not account for X-rays, but the inclusion of X-rays causes sufficient photo- and Auger ionization to populate the O VI and N V states. The overall wind acceleration remains about the same for both cases (see black and cyan lines in Fig. 13). A further increase in the X-ray field to L_x ∼ 10³⁶ erg s⁻¹ affects the emerging UV spectrum drastically. The strength of the wind lines (such as Si IV, C IV, C III, and P V) reduces significantly whereas O VI becomes stronger. A model including a strong X-ray field on the order of L_x ∼ 10³⁸ erg s⁻¹ considerably reduces the strength of Si IV and C IV lines, making them appear to be photospheric lines.

Fig. 12. Comparison of model SEDs (left upper panel) and emerging UV spectra (lower panel) for different X-ray illuminations. The model SEDs are compared with XMM-Newton and HST spectra taken at phase ϕ = 0.53. A zoom on C IV and Si IV line profiles are displayed in the upper panels (middle and right).

Since the important wind-driving ions are depopulated, the line acceleration in the outer wind decreases correspondingly. As shown in Fig. 13, strong X-ray illumination causes a strong decline in the wind accelerations from the sonic radius onward and hence stalling of the wind. This is because photoionization destroys the ions that are capable of absorbing photons in the UV. The normalized wind acceleration drops to unity for models with L_x ∼ 10³⁶ erg s⁻¹. The wind driving drops further below unity for models irradiated by L_x ∼ 10³⁸ erg s⁻¹ (see red line in Fig. 13). The total wind acceleration is a factor of seven lower at the side of the donor facing BH and a factor of three lower in the shadow region compared to that of an unperturbed O-star. Our investigations thus suggest a breakdown of wind toward the BH in M33 X-7. For this model, the blue O VI and N V absorption components only reach up to 1000 km s⁻¹ instead of 1500 km s⁻¹ measured at eclipse. The decrease of the wind terminal velocity due to X-ray irradiation is supported by previous HMXB studies (see e.g., Watanabe et al. 2006; Sander et al. 2018; Krtička et al. 2018).

Fig. 13. Comparison of total wind acceleration normalized to the gravitational acceleration for models including different amounts of X-rays. The location of the sonic and the BH radius are marked with dashed and solid gray lines.

6.2. Wind driving

In order to have a detailed picture of wind driving in these situations, we need to look at the contributions of the various elements in different ionization stages. This is illustrated in Figs. 14 and 15 by showing the radiative acceleration contribution from different ions and how this is affected by X-rays. While these are not dynamically consistent models, the changes due to the X-ray illumination are representative nonetheless. For a model without including X-rays (Fig. 14 upper panel), the leading ions are Fe IV, N III, and C III. O III becomes dominant in the far outer part of the wind. Near the sonic point, iron group opacities are the leading contributors to the wind along with the electron scattering and gas pressure. Beyond R > 1.1R_*, CNO and other higher elements become important agents in wind driving. Higher ions such as Si IV, C IV, S IV and N IV also make a significant contribution to the radiative acceleration. The contribution by C, N, O, and iron lines outweigh the free electron (Thomson) scattering contribution in the outer wind.

Fig. 14. Contributions to the radiative acceleration from the different ions for models with no X-ray (top) and including Lx ∼ 1036 erg s−1 (bottom). Symbols and colors representing different elements and their ionization state is shown in the figure. The plot is limited to the ions that contribute at least ∼1% to the radiative acceleration. Line acceleration due to scattering by free electrons (magenta solid curve) and gas pressure (orange dashed curve) and total radiative acceleration (red dashed curve) are plotted for comparison. The radius at which the BH is located is marked by a gray solid line and the location of the sonic point marked by gray dotted line.

Fig. 15. Same as Fig. 14, but for model including Lx ∼ 1038 erg s−1.

In contrast, the model with L_x ∼ 10³⁶ erg s⁻¹ (Fig. 14 lower panel) depopulates all the lower ions, especially beyond the radius at which the BH is located. In the inner part, the picture remains the same as for the model with no X-rays. Beyond R > 1.1R_*, the acceleration by Fe IV reduces rapidly and becomes unimportant. In this case O VI and N V are the leading contributors in the outer wind as we see can from the emerging UV spectra in Fig. 12. The total radiative acceleration by these ions is comparable to that of free electrons. Si IV and C IV also make a noticeable contribution to the wind driving, but their influence is a factor of three lower compared to the model with no X-ray ionization. Increasing the X-ray field to L_x ∼ 10³⁸ erg s⁻¹ makes a remarkable further difference as shown in Fig. 15. Strong X-ray photoionization results in the super-ionization of lower species hence only a few highly ionized species are contributing to the wind driving. O VI is still the leading ion in this case, and other higher ions such as Fe VII, N IV and O VII are also important. However, the contributions from these ions are lower than that by electron scattering. In comparison to the other two scenarios, here we can see a completely different picture in the inner part of the wind, where a plethora of ions contributing to the wind driving including the iron group components, are missing. It should be noted that apart from CNO and Fe, elements like P, S, Si, Mg, and Ne also make important contributions in wind driving. Excluding higher elements and higher ionization stages would significantly affect the wind structure and the emerging spectra.

As the calculation of hydrodynamically consistent atmosphere models is beyond the scope of the present paper, one has to be careful with drawing deeper conclusions about the absolute radiative acceleration in the deeper wind layers. However, our studies give us very valuable insights regarding relative changes of the wind acceleration between the different situations. Without X-rays, the acceleration in the outer wind is approximately constant with 1 < Γ_rad ≈ 3, which is quite similar to the result for the much cooler Vela X-1 donor obtained with a dynamically consistent model (Sander et al. 2018). With the amount of X-rays measured during eclipse in M33 X-7, the situation reduces to Γ_rad ≈ 1. Taking this result literally would imply that there is hardly any wind acceleration, but as we do see wind signatures in the UV spectrum during eclipse we know that Γ_rad > 1 must be fulfilled at least in an inner region where the wind would be launched and which we do not model dynamically consistent in this work. The consistent model from Sander et al. (2018, their Fig. 6) reflects such a situation where even for an illumination of L_X ≈ 10³⁷ there is still a meaningful radiative acceleration close to the star before a complete breakdown occurs further out. Moreover, the inner acceleration in our M33 X-7 models is not affected when comparing the model without X-rays to those using the L_X from the eclipse. This changes drastically, when looking at the strong irradiation case where the full L_X ≈ 10³⁸ from the BH are irradiated. Now, even the acceleration in the inner region decreases by about 40%. Close to the BH, the acceleration further decreases by another 60% compared to the eclipse situation. With Γ_rad < 1 now, the acceleration of the wind breaks down. For an isolated star, this would lead to a significant deceleration of the wind, as also discussed by Krtička et al. (2018, 2022) for a wider range of system parameters. On top, the presence of the BH and its gravitational potential will need to be considered here, which we will discuss in Sect. 7.

6.3. X-ray luminous HMXB

The typical wind-intrinsic X-ray emission of an O-giant donor is L_x ∼ 10³² (Nebot Gómez-Morán & Oskinova 2018). However, this is much lower than the observed X-ray luminosity at eclipse estimated by Pietsch et al. (2004) using XMM-Newton observations. They derived L_x = 2 − 3 × 10³⁶ erg s⁻¹ at phase ϕ ≈ 0. Such a high X-ray luminosity could arise from the photoionized region around the accretion source if the inclination angle of the system is low or if the size of the donor is too small to fully cover the photoionization zone. However, M33 X-7 is an eclipsing system with a high inclination angle (≈75°) and the donor has a large radius over-filling its Roche-lobe. Typically in Roche-lobe overflow systems, due to the high X-ray luminosity, the stellar wind is interrupted by the photoionization. Thus, the wind can only develop in the X-ray shadow behind the OB supergiant, a so-called “shadow wind” region (e.g., Blondin 1994; Kaper & van der Meer 2005). In M33 X-7, the strong X-ray luminosity near phase ϕ ≈ 0 suggests the presence of a significant ionization zone in the line of sight of the donor with the shadow wind region being very limited. Such an extended Strömgren zone could be formed due to the scattering of X-rays by ions in the stellar wind. As a result, the absorption and emission components of UV wind lines are diminished. By incorporating such a high amount of X-rays we are able to reproduce the observed UV spectra at the eclipse as shown in Fig. 6.

As discussed in Sect. 2.2, we estimated the X-ray luminosity of the system at phase ϕ = 0.33 and 0.5 to be around L_x = 10³⁸ erg s⁻¹. Hence, by including X-rays with L_x ∼ 10³⁸ erg s⁻¹ and adopting the same mass-loss rate (log Ṁ ∼ −6.2), our model SED and spectra are in agreement with the XMM-Newton and HST observations at phase ϕ = 0.33 and 0.53 (see Figs. 12, 16, and 17). Since the X-ray and UV spectra at phase 0.33 and phase 0.53 are similar and the UV wind lines in these phases suffer large reductions in the absorption strength at all velocities, we could not find any evidence for an accretion wake or focused wind stream toward the BH in this system.

Fig. 16. Same as Fig. 6, but spectra taken at inferior conjunction (ϕ ≈ 0.53) where BH is in front of the donor star.

Fig. 17. Same as Fig. 6, but spectra taken at Quadrature (ϕ ≈ 0.33).

HMXBs with high X-ray luminosities such as SMC X-1 and LMC X-4 (OB supergiant + neutron star) show similar orbital-dependent variations in the UV spectra. They show a similarly high X-ray luminosity around deep eclipse (∼1% of peak). As a result, their corresponding UV spectra show strong P Cygni profiles of O VI and N V (Vrtilek et al. 2005; Sonneborn et al. 2007). These systems are also reported to have a reduction in the strength of O VI and N V due to two orders of magnitude higher X-ray luminosity in other phases. This is exactly what is predicted by the models for M33 X-7 in Fig. 12. Unfortunately, the UV spectra of M33 X-7 do not cover these O VI and N V lines. In conclusion, our models are capable of explaining the UV wind line variation of high luminosity HMXBs.

Moreover, we varied Ṁ in models including X-rays to check its impact on the emergent spectral profiles. However, these results either over or under-predict the strength of the UV wind lines compared to the observed spectra. By fixing mass-loss rate and wind velocity, we have also tested models with different β values. The standard β = 0.8 gives the best-fit model whereas higher β values could not reproduce the observed UV line profiles. A low β value implies a strong acceleration of the wind, hence the wind velocity near the BH (at 1.75R_*) could reach up to 760 km s⁻¹ (see Fig. 18). In comparison, the orbital velocity (2πa/P) at this radius is only 510 km s⁻¹. However, the actual wind velocity near the BH might very well be comparable to the orbital velocity or even lower as we assume a prespecified wind velocity field in this work. In this approach, the effect of the X-ray illumination from the BH on the spectrum formation is taken into account, but there is no immediate coupling of the velocity field to the accompanying lower radiative acceleration derived in Sect. 6.2.

Fig. 18. Velocity field for M33 X-7 donor adopting β = 0.8. The location of BH and the corresponding wind velocity and orbital velocity are marked.

In conclusion, we prefer our final best-fit solutions to M33 X-7 using X-ray irradiated models over other possible solutions (see Figs. 10 and 11) since they simultaneously reproduce the X-ray, UV, and optical observations. Moreover, we are able to account for phase-dependent variations in the spectra with the same wind parameters. The wind parameters of our best possible solution are listed in Table 4.

6.4. He II emission

One of the main sources of the He II emission seen in high redshift galaxies is thought to be HMXBs. Strong X-ray photoionization in HMXBs causes an excess flux beyond the He ionizing edge. This is demonstrated in our model SEDs in Fig. 19. The number of He II ionizing photon flux provided by the O star donor model alone (black lines in Fig. 19) is substantially lower than those irradiated with X-rays. The He II ionizing photon flux generated by M33 X-7 is found to be approximately seven orders of magnitude higher than that from a normal O star of the same spectral type.

Fig. 19. Spectral energy distribution for models with different X-ray luminosities. The corresponding number of He II ionizing photons are displayed in the legend.

Maggi et al. (2011) have reported nebular He II 4686 Å emission in the region. They inferred the luminosity of the line to be ∼2 × 10³⁵ erg s⁻¹. Using a basic photoionization equation from Osterbrock & Ferland (2006), we calculate the He II photon flux to be $Q_{\mathrm{He}\mathrm{I}\mathrm{I}}=1.02\times {10}^{48}\frac{L(\mathrm{He}\mathrm{I}\mathrm{I}4686)}{{10}^{36}\mathrm{erg}{\mathrm{s}}^{-1}}=2\times {10}^{47}\mathrm{ph}{\mathrm{s}}^{-1}$. This is in good agreement with our model predictions for M33 X-7 considering an X-ray irradiation of L_x = 10³⁸ erg s⁻¹ out of eclipse.

Interestingly, the UV spectra does not show an absorption line at He II 1640 Å despite this being predicted by the model (see Fig. 20). This absence of absorption suggests that there is possibly an emission component filling the line. This feature is also variable with orbital phases, with the maximum “emission” occuring when BH is in our line of sight (green and blue lines in Fig. 20). Such He II emission features are not expected for a typical late O giant star. Although we account for a maximum X-ray luminosity of log(L_X/erg s⁻¹)≈38 in the model which predicts He II ionizing photons as high as 10⁴⁷ s⁻¹, we were unable to reproduce these observed emission features. Broad He IIλ1640 emission features were detected in another HMXB HD 153919 (alias: 4U1700-37) with a neutron star companion (Kaper et al. 1990). These lines are expected to originate from EUV emission lines close to 256 Å, which is Raman scattered by He II ions in the stellar wind. Kaper et al. (1990) show that these emission lines originate close to the compact object, causing the orbital modulation of the Raman scattered lines. However, HD 153919 also shows a wind-instrinsic P Cygni profile in He IIλ1640 (Hainich et al. 2020), which never appear in our M33 X-7 UV spectra, even at X-ray eclipse.

Fig. 20. Comparison of UV spectra of M33-X7 around He II 1640 Å taken during different orbital phases. Line colors indicate different orbital phases as given in Fig. 9. The dashed black line represent our best fit model including Lx = 1038 erg s−1, corresponding to ϕ = 0.53.

7. Accretion rate

Accretion of matter by the compact object is the most efficient way of producing X-rays. In wind-fed supergiants, X-ray accretion is powered by the stellar wind of the donor. Accretion of matter onto a star moving through a medium was first described by Bondi & Hoyle (1944) and applied to HMXBs by Davidson & Ostriker (1973). The mass accretion rate via the Bondi-Hoyle formula can be written as follows:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{acc}}=\pi R_{\mathrm{acc}}^2\rho v_{\mathrm{rel}}\end{array}$$(1)

Here ρ is the wind density near the BH and R_acc is the accretion radius.

In an HMXB, the compact object has a velocity relative to the stellar wind flow:

$$\begin{array}{c}\hfill v_{\mathrm{rel}}^2=v_{\mathrm{wind}}^2+{[v_{\mathrm{orbit}}-(R_{\ast }/d_{\mathrm{BH}})v_{\mathrm{rot}}]}^2.\end{array}$$(2)

Here v_wind and v_orbit are the wind velocity and orbital velocity near the BH which are 760 km s⁻¹ and 513 km s⁻¹ respectively. The eccentricity of the orbit is ignored in this case because it is essentially zero for this system. Since the O star is rotating fast while the wind is leaving its surface, there is an additional tangential component to its velocity. Plugging in our derived values from the analysis, the estimated velocity is v_rel ≈ 845 km s⁻¹. The derived relative velocity is comparable to 0.56 v_∞. Stellar wind matter is accreted onto the compact object if it approaches within an accretion radius, which can be approximated to

$$\begin{array}{c}\hfill R_{\mathrm{acc}}={\displaystyle \frac{2G{M}_{\mathrm{BH}}}{v_{\mathrm{rel}}^2}}.\end{array}$$(3)

Using the velocity calculated above, the accretion radius in M33 X-7 is ∼4 × 10⁶ km = 5.7 R_⊙. This is roughly 0.5R_* from the photosphere of the donor star. The density of the wind at the distance of the BH can be approximated to

$$\begin{array}{c}\hfill \rho ={\displaystyle \frac{{\dot{M}}_{\mathrm{donor}}}{4\pi d_{\mathrm{BH}}^2{v}_{\mathrm{wind}}}}.\end{array}$$(4)

By substituting the derived wind parameters, we get a density of ρ ≈ 6.7 × 10⁻¹⁵g cm⁻³, and based on Eq. (1), we get Ṁ_acc = 2.8 × 10¹⁷ g s⁻¹. The expected X-ray luminosity can then be estimated using the relation L_X = ηL_acc with the accretion luminosity

$$\begin{array}{c}\hfill L_{\mathrm{acc}}={\displaystyle \frac{G{M}_{\mathrm{BH}}{\dot{M}}_{\mathrm{acc}}}{R_{\mathrm{BH}}}},\end{array}$$(5)

and η being an accretion efficiency parameter.

In BHs, the accretion disk luminosity will be created while the material falls in from the outermost radius to the innermost stable circular orbit. For a nonrotating BH, the inner edge of the accretion disk should be at three times the Schwarzschild radius $R_{\mathrm{sch}}={\displaystyle \frac{2G{M}_{\mathrm{BH}}}{c^2}}$. For a BH of 11 M_⊙, the accretion radius would be around 32 km. For a rapidly spinning BH, the inner edge of the accretion disk can be considerably closer to the event horizon than for a nonrotating BH. The accretion efficiency would be in the range of ∼6 − 40%. Therefore, the accretion luminosity can be in the range 1 − 8 × 10³⁷ erg s⁻¹. Previous studies have estimated the spin of M33 X-7 to be 0.84 (Liu et al. 2008, 2010). However, the considerable reduction in the mass of the BH can lower the spin parameter. Considering the effect of the BH mass on the spin parameter (see Fig. 3 in Liu et al. 2008), the updated spin of M33 X-7 is around a_* ≈ 0.6. Therefore, the corresponding accretion luminosity is only ∼4 × 10³⁷ erg s⁻¹.

The estimated accretion luminosity corresponds approximately to the observed X-ray luminosity for M33 X-7 at eclipse (Pietsch et al. 2006, 2004) but is a factor of four lower than our estimated values for the other two phases out of eclipse (Sect. 2.2). In comparison, the Eddington luminosity of the system is ∼1.4 × 10³⁹ erg s⁻¹. So the maximum observed luminosity is only ≈12% of the Eddington value, which excludes super-Eddington-accretion.

Possible explanations for the observed high X-ray luminosity would be excess mass transfer due to classical Roche-lobe overflow (RLOF) in the system or so-called wind-RLOF, where the radius at which the wind is accelerated is comparable to the Roche-lobe radius so that the wind material fills the donor’s Roche-lobe and is transferred to the compact companion through the inner Lagrangian point (Mohamed & Podsiadlowski 2007). The different modes of mass supply likely yield differences in the observed X-ray luminosity of X-ray binaries. Typically, HMXBs which accrete only via stellar winds show a constant X-ray radiation with luminosities in the range 10^{35 − 37} erg s⁻¹, while systems associated with RLOF have higher luminosities (L_x ≳ 10³⁸ erg s⁻¹).

In M33 X-7, the O supergiant companion is in a stage where it is already filling the traditionally defined Roche lobe. However, the presence of a radiative acceleration Γ_rad can significantly alter the Roche potential (e.g. Howarth 1997; Dermine et al. 2009). To get a better idea about this, we can use our insights gained from the radiative acceleration study in Sect. 6.2. As discussed in Dermine et al. (2009), there is no longer a Roche lobe for the donor star in regions with Γ_rad > 1. We know from UV spectroscopy that this must be the case for M33 X-7, at least in some region not facing the BH. Hence, some of the material will directly escape from the system and is not re-directed toward the BH. For Γ_rad > 1, there is no longer a clearly defined Roche lobe and thus the destinction between wind capturing and RLOF actually vanishes.

For M33 X-7, the situation is likely even more complex. Even the concept of wind-RLOF (Mohamed & Podsiadlowski 2007) does not really apply here as it would require Γ_rad ≲ 1 to define a Roche lobe that could be filled by a slow wind as we typically find them in late-type stars. Instead, as our acceleration study shows, our donor star likely has a region with Γ_rad ≪ 1 facing the BH. In this case, there is a definable modified Roche potential with generalized radius, which will be smaller than the one for the unperturbed Roche lobe (see Dermine et al. 2009). Since the derived donor star parameters for M33 X-7 are already sufficient to yield a star that would fill the standard definition of a Roche Lobe, the further reduction of this criterion for a region towards the BH paves the way for an additional “easy” overflow-like mass transfer, that is further boostered by the diminished radiative acceleration of the material due to the X-rays from the BH. In contrast to the donor, a meaningful Roche lobe can still be defined for the accretor, which would play a major role in capturing a considerable (> 10%) fraction of the wind as indicated in recent hydro-dynamical simulations by El Mellah et al. (2019).

Based on our atmosphere modeling, the maximum wind velocity outside eclipse should be v_∞ ≈ 1000 km s⁻¹, which is similar to other high-L_X HMXBs. Hence, we can assume that the wind speed at the distance of the BH is lower than in our eclipse model (cf. Fig. 18), i.e., close to the orbital speed. Subsequently, the accretion luminosity in the Bondi-Hoyle scenario already increases to ∼1.4 × 10³⁸ erg s⁻¹, in better agreement with the X-ray observations. With the additional mass overflow in the highly ionized region, we have a plausible “wind overflow” scenario for the observed L_X, even taking into account the efficiency factor between L_acc and L_X. While our scenario likely explains the higher accretion luminosity compared to the classical Bondi-Hoyle wind-fed scenario, the condition of Γ_rad > 1 in other directions will avoid a standard RLOF where all of the donor’s mass loss is redirected toward the BH. Consequently, our derived L_X is on the lower end compared to systems that are considered to be powered by standard RLOF.

Recent studies of other BH HMXB systems, namely Cyg X-1 (Orosz et al. 2011; Miller-Jones et al. 2021) and LMC X-1 (Orosz et al. 2009; Hyde et al. 2017), also suggest that the donor stars in these systems are close to filling their Roche lobe. Yet, a detailed consideration of the radiative acceleration has so far not been done. Depending on the impact of the X-ray illumination of the donor atmosphere, these systems will either have a similar “wind overflow” like M33 X-7 or a standard RLOF in case of strong wind inhibition.

8. Formation and evolution of the HMXB

M33 X-7 has been identified as an evolutionary challenge, considering the previously reported very high mass donor star along with a heavy BH in a tight orbit. The updated stellar and binary parameters allow us to gain new insights into the evolution of this system. In the first step, we explore the formation channel of the system starting from a binary with two stars on the zero-age main-sequence (ZAMS) to a BH – O star binary and then follow its path to later evolution.

As a first approach, we compare the updated parameters of the O star donor to that of single-star evolutionary tracks. For this purpose we used tracks provided by Brott et al. (2011) at LMC metallicity along with the Bayesian statistic tool “The BONN Stellar Astrophysics Interface” (BONNSAI³, Schneider et al. 2014). The results from the analysis are summarized in Table 5. Single-star evolutionary tracks suggest lower masses for the O star donor given its observed luminosity, i.e., the donor is underluminous compared to a normal single star of the same mass. The observed wind mass-loss rate of the donor is in good agreement with the prediction of Vink et al. (2000, 2001) which is used in the stellar evolutionary models. The single star models could account for the fast rotation of the star by assuming an initial velocity of 280 km s⁻¹. Abundances are also reproduced within the uncertainty limit except for helium.

With our knowledge of the currently observed properties of the M33 X-7 system, we tried to reconstruct its evolution from the moment the system was comprised of two ZAMS stars through the intermediate phases such as mass transfer from the primary to the secondary, and the formation of the BH, followed by the currently observed HMXB phase and the final evolution of the system.

8.1. Modeling the binary system with MESA

We performed detailed binary evolution calculations to explore the possible evolutionary history of M33 X-7. The tracks were computed using the stellar evolution code MESA (“Modules for Experiments in Stellar Astrophysics”, Version No. 10398) with a physics implementation as described in Paxton et al. (2011, 2013, 2015, 2018). We adopted LMC metallicity (Z = 0.006) for the calculation, which is in line with the derived metallicities based on the spectroscopic analysis. To make our models comparable to the Brott (2011) LMC models used for the single star evolution in BONNSAI, we modified MESA as described in Marchant et al. (2016).

Stellar winds were implemented as in Brott et al. (2011), with mass-loss rates for hot hydrogen-rich stars modeled as in Vink et al. (2001). For stars with surface hydrogen of X_H < 0.4, we use the mass-loss rate of Hamann et al. (1995) divided by a factor of 10 to account for clumping (Yoon et al. 2010). We interpolate these two mass-loss rates for surface hydrogen fractions between 0.7 and 0.4. For temperatures below the bi-stability jump, the mass-loss rate was taken as the maximum between the Vink mass-loss rate and that of Nieuwenhuijzen & de Jager (1990).

Convection was modeled using the mixing-length theory (Böhm-Vitense 1958) with a mixing-length of α_ML = 1.5. Semiconvective mixing is modeled with an efficiency of α_SC = 1.

Rotational mixing and angular momentum transport were modeled as diffusive processes following Heger & Langer (2000), including the effects of Eddington-Sweet circulations, the Goldreich-Schubert-Fricke instability, secular and dynamical shear, and sensitivity to composition gradients. We also take into account the transport of angular momentum due to the Spruit-Tayler dynamo (Spruit 2002).

Stars are assumed to be synchronized with the orbital period at the ZAMS. The evolution of orbital angular momentum considers the effects of mass loss, gravitational wave radiation, and spin-orbit coupling, as described in Paxton et al. (2015).

8.2. Evolutionary status of M33 X-7

To investigates possible scenarios of M33 X-7’s evolutionary history, we considered progenitor masses in the range of M₁ = 38 − 60 M_⊙ for the primary, M₂ = 18 − 40 M_⊙ for the secondary, and initial orbital periods from 1.8 to 5 days. Guided by the mass derived from the spectroscopic analysis of the donor, which is higher than what is expected from evolutionary models of single stars, we searched for binary models where the secondary accretes mass from the BH progenitor. Given the short orbital period, we considered only binary systems that start their life in a tight orbit, hence undergoing a phase of mass transfer during the main sequence phase of the primary. In our model calculation, we allow both stars to evolve up to core carbon depletion. After core carbon depletion, a star is assumed to directly collapse into a BH without losing mass or the orbit receiving a kick due to a possible supernovae explosion. In this approach, we attribute the final CO mass to the compact companion and treat it as a point mass. Subsequently, we model the evolution of the secondary together with the compact object in the binary system.

While exploring the parameter space, we discovered that (1) models with high initial masses and short orbital periods are predicted to undergo overflow through the second Lagrangian point L2 during core-H burning and are expected to merge; (2) models with lower mass ratios (q ≲ 0.5) and short periods (P < 5 days) had numerical issues during the fast Case A mass transfer phase; and (3) models with mass ratios going toward unity (q ≳ 0.85) have a reverse mass transfer from the secondary to the primary while both stars are still core-H burning, leading to an unstable situation. These models terminate before forming a BH + O star binary, and are therefore ruled out, limiting the possible parameter space. For some models that successfully reach the HMXB stage we found that they form a convective envelope during mass transfer on the compact object, leading to unstable situations followed most likely by a common envelope phase and a merger. We think that these models are unlikely to represent the true nature of this system which is believed to be currently in a stable mass transfer phase.

Based on our tests, the resulting models exhibit the following properties at the HMXB phase: M_OB ≈ 30 − 40 M_⊙, M_BH ≈ 8 − 13 M_⊙, and P ≈ 4 − 12 days. Our best possible evolutionary model which matches the current parameters of M33 X-7 is shown in Fig. 21 (Binary model 1). The parameters derived from this model are compared to the observations in Table 5. This model suggests the initial masses of the stars to be M₁ = 45 , M₂ = 32 M_⊙ with a period of 2.7 days. Our best model satisfies the current observational properties (T_*, log g, log L, R_*) of the donor star. Moreover, the evolutionary masses suggested by the model for the O-star donor and the BH are in good agreement with the observations.

Fig. 21. Best matching MESA evolutionary tracks reproducing the HRD position of O-star donor in M33 X-7. The current position and errors of the O star are taken from the spectroscopic analysis. Initial masses and periods are labeled on the plot.

The initial masses of the M33 X-7 progenitors result from our tests are comparable to typical main-sequence O-stars and are much lower than the predictions from previous studies (Abubekerov et al. 2009; Valsecchi et al. 2010; de Mink et al. 2010; Bogomazov 2014). According to the evolutionary model, the donor is overflowing its Roche lobe and transferring mass to the BH at the currently observed stage as a HMXB. As a result, the total mass-loss rate of the system is an order of magnitude higher in this model, despite the wind mass-loss rate being in agreement with the empirical estimation. This would qualitatively be in line with out findings in Sect. 7, where we conclude an enhanced mass transfer due to the strong X-ray illumination and the perturbance of the Roche potential due to the presence of a donor wind. Of course such a detailed consideration of the wind is beyond the MESA treatment where standard Roche-lobe overflow is assumed instead.

Our suggested binary model 1 does not reproduce the fast rotation of the secondary star (the current O-star donor). The current period of the system suggested by the model is slightly higher than the observed period. This is mainly because the initial orbit got widened as a result of the strong mass-loss rate of the primary during the assumed Wolf-Rayet phase. When reducing the initial period, we find no stable solution during the first mass-transfer phase. However, adopting a reduced mass-loss rate in the WR stage compared to the standard prescriptions as recently derived for lower metallicities by theoretical calculations (Sander & Vink 2020; Higgins et al. 2021) might be able to resolve this conflict.

Considering the predicted abundances, we find that He is by far not as enriched in our model as derived from our observations. There are major effects that can lead to altered abundances, mixing and losing mass, so that the inner enriched layers are exposed to the observer. As the fast rotation rate of the OB star is not reproduced either, these two findings are likely linked and the evolutionary abundances suffer from missing the additional rotational mixing.

Figure 22 illustrates possible evolutionary sequences of M33 X-7 based on the results from our best-match model. The progenitor comprises of a primary with 45 M_⊙ (BH progenitor) and a secondary with 32 M_⊙ (O-star donor progenitor). At this point, both stars are on the main sequence. During the first ≈3 Myr the system remains detached. After core-hydrogen exhaustion, the primary leaves the main sequence and starts to expand quickly. When its radius approaches the Roche lobe, a rapid phase of mass transfer initiates to the less massive secondary which still resides on the main sequence. The secondary accretes most of the transferred mass and acquires large angular momentum so that it can spin up to critical velocity. During this phase, the orbit shrinks to 2.5 days. Even though the secondary spun up during the mass transfer, the rotation rate drops quickly to the earlier values after this stage due to tidal locking. After the rapid interaction the system experiences episodes of slow mass transfer until ≈4 Myr. The mass transfer ends when the primary’s hydrogen envelope is partially stripped, with only ≈30 M_⊙ left. The primary will evolve to the core-He burning phase, transitioning to a Wolf-Rayet star with a strong stellar wind. At this point, the mass of the He core would be around 15 M_⊙. Due to its strong mass loss, the orbital period would increase to 5.5 days. The Wolf-Rayet phase will last for 0.6 Myr. The remaining He core of the primary eventually collapses into a BH at 5.6 Myr.

Fig. 22. Illustration of possible evolutionary scenario for M33 X-7 starting from the ZAMS massive binary to HMXB stage followed by common envelope merging.

At around 5.76 Myr, the scenario described above forms a HMXB, resembling the properties of M33 X-7. The total duration of this phase lasts for ≲0.2 Myr, implying that this stage has a reasonable probability to be observed. At this stage, similar to our observations, the O-star donor has already filled up its Roche lobe, and the BH is accreting this material (though see Sect. 7).

According to our model, the mass-transfer will become unstable in a few 10 000 yr and the system will merge. However, a more fine-tuned model (e.g., with different mixing) might lead to a stable solution. Nonetheless, since the mass ratio is large and the orbital period is short, it does not seem to be unreasonable that the system will eventually merge⁴. Under this assumption, our system will not be a progenitor of a double-BH binary. Instead, the BH and the He core will form a common envelope with an in-spiral of the BH eventually resulting in a merger, leading to the formation of a more massive BH. This could be accompanied by a γ-ray burst.

As a slightly different alternative we show a second binary model (Binary model 2) that also reproduces the observed parameters quite well while having a different fate. The results are shown in Table 5 and the corresponding evolutionary tracks are plotted in Fig. 23. In this scenario, we have discarded our assumption of tidal synchronization in the binary system. However, it is most likely that the short period binaries would evolve into tidal locking very soon, independent of the initial stellar spins. In binary model 2, we assume an initial rotation velocity of 300 km s⁻¹ for both stars. Due to the change in our assumptions on rotation, now the secondary is only accreting a small fraction of mass during the first mass transfer. Contrary to the earlier model, the rotation rate of the secondary reduces only gradually after the spin up. Hence, this model resembles the observed fast rotation of the O-star donor at the current HMXB phase. Subsequently, the primary quickly proceeds to the Wolf-Rayet phase at around 4.5 Myr with a He-core mass of ≈20 M_⊙. At 5 Myr the primary collapses to a 12 M_⊙ BH. At this stage, the orbit has already widened due to the strong wind mass loss. The HMXB phase in this model happens earlier than in the previous model. However, the period of the system is almost a factor two higher than that of M33 X-7. In addition, the mass of the donor is lower whereas the BH mass is slightly higher compared to our empirical estimations. This translates to a lower mass ratio for the system (q ≈ 2.8). Because of this slight decrease in the mass ratio and increase in the orbital period, the second model predicts a stable mass transfer phase after the donor fills its Roche lobe at 5.4 Myr, showing the complexity of stellar evolutionary modeling. At this point, the orbital period shrinks down to two days. Subsequently, the secondary starts to contract and evolves further as a Wolf-Rayet star. Throughout these final phases, the primary BH accretes and gains around 1 M_⊙ via stellar winds. At 6.3 Myr, the secondary will reach core carbon burning exhaustion. The final system consists of a BH-BH binary with 13 and 9 M_⊙ in a four-day orbit. Based on their large separation, these BHs will not merge in a Hubble time.

Fig. 23. Similar to Fig. 23, but an alternative evolutionary track assuming no tidal synchronization.

On the elusive detection of M33 X-7-like systems. Although the initial masses of the M33 X-7 progenitors are comparable to that of typical O-star binaries, such BH HMXB systems are rare. Using binary evolutionary model calculations, we investigated the formation of OB + BH binaries with a range of masses and periods. According to these models, the duration of the OB + BH system is only ≈2 − 5% of the total lifetime of the binary before both stars become compact objects. Moreover, such a system would be observed as X-ray binaries only when an accretion disk forms around the BH. The accretion disk formation criterion for wind-fed systems mainly depends on the wind and orbital velocity near the BH radius and on the mass ratio (Illarionov & Sunyaev 1975; Shapiro & Lightman 1976). Thus, BH HMXBs are mostly limited to short-period systems with donors evolved from the main sequence. With the same observed parameters as of M33 X-7, a system would not be detectable in X-rays if the orbital period is > 15 days. Besides, if the donor is a main-sequence star with a smaller radius and a higher terminal wind velocity than in M33 X-7, the accretion disk will only form when the orbital period is less than its current value.

9. Conclusions

We present a detailed spectroscopic analysis of the only known eclipsing BH HMXB system, M33 X-7. For this study we obtained the first UV spectra of the system accompanied by X-ray observations. The observations were carried out at three key orbital phases, tracing the BH eclipse, egress, and inferior conjunction using the HST. At the latter two phases, X-ray observations were obtained simultaneously using XMM-Newton. These phase-resolved observations shed light on the interaction of the stellar wind with the BH. The UV resonance lines show the Hatchett-McCray effect with a large reduction in absorption strength when the BH is in the foreground due to the strong X-ray ionization. We performed a detailed analysis of X-ray, UV, and archival optical spectra using stellar atmosphere models, arriving at the following conclusions:

• The donor is much cooler and dimmer than previously predicted. The mass-accreted O star shows signs of He enrichment and fast rotation.

• Our detailed analysis considering multi-wavelength data suggests a large reduction in the mass of the system compared to previous results. The mass of the donor star is reduced from 70 M_⊙ to 38 M_⊙. Correspondingly, the BH mass changes from 15.6 M_⊙ to 11.4 M_⊙ and spin of the BH decreases to a_* ≈ 0.6.

• Based on the newly derived parameters of the donor and mass of the BH, the orbital separation has come down to 35 R_⊙. The comparison of the radius of the star with a standard Roche-lobe approximation suggests an overfilling factor of f ≈ 1.2.

• The strong X-ray radiation of M33 X-7 ionizes much of the wind and leaves a large Strömgren zone. These extreme ionization conditions prevent the formation of a normal radiatively driven wind from significant parts of the donor. Our one-dimensional calculations confirm that the photoionization by the X-ray radiation can significantly change the ionization structure and diminish the wind accelerations. Our investigations on wind driving and the impact X-ray irradiation in M33 X-7 can also be applied to other high-luminosity HMXB systems in general.

• We estimated the wind parameters of the donor based on the UV spectra taken at X-ray eclipse. The derived mass-loss rate is in good agreement with the theoretical prediction assuming a depth-dependent microclumping. By incorporating X-rays corresponding to the observations at different orbital phases, we were able to reproduce the observations at three different phases including eclipse with the same mass-loss rate.

• The observed X-ray luminosity derived from our analysis is ∼2 × 10³⁸ erg s⁻¹, which is about 12% of the BH Eddington luminosity. The accretion luminosity derived using the Bondi-Hoyle calculation is found to be much lower than the observed X-ray luminosity. Fostered by the wind inhibition due to the X-ray photoionization, we instead have a “wind overflow” scenario where the BH not only captures a part of the wind, but it is further fed by an additional mass overflow along the inner Lagrangian point. In other directions, the presence of a significant stellar wind significantly perturbs the Roche potential, allowing material to escape from the system. As a consequence, the classical distinction between wind-fed and Roche-lobe overflow systems becomes meaningless for our system. The observed X-ray luminosity is exactly in line with such a “wind overflow” scenario.

• The X-rays contribute to the formation of the He II emission region around the system emitting ∼10⁴⁷ ph s⁻¹.

Taking the parameters of M33 X-7 from the observations into account, we attempted to trace its evolution using binary evolution models. We modeled the system from a ZAMS binary system until its current state. In addition, we proposed two possible scenarios for its future evolution. The system is transitioning toward an unstable mass transfer phase, resulting in a common envelope with a deep spiraling of the BH to the envelope of the massive donor. Since the mass ratio is q ≳ 3.3 and the period is short, the common envelope phase is likely to result in the merger of two companions resulting in a heavier BH.

Acknowledgments

The authors are grateful to the referee for the detailed and constrictive report which helped to improve the paper, and also for providing the optical spectroscopic data. We thank O. Kargaltsev for the useful discussions. We thank the Swift telescope team for approving our ToO request. This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester. VR acknowledges support by the Deutsches Zentrum fsür Luft- und Raumfahrt (DLR) under grant 50 OR 1912. VR and AACS acknowledge support by the Deutsche Forschungsgemeinschaft (DFG – German Research Foundation) in the form of an Emmy Noether Research Group (grant number SA4064/1-1, PI Sander). DP acknowledges financial support by the Deutsches Zentrum für Luft und Raumfahrt (DLR) grant FKZ 50 OR 2005. This publication has benefited from discussions in a team meeting sponsored by the International Space Science Institute at Bern, Switzerland. This research made use of the VizieR catalog access tool, CDS, Strasbourg, France. The original description of the VizieR service was published in A&AS 143, 23.

References

Appendix A: Additional figures

Fig. A.1. Comparison of synthetic spectra in the blue optical range obtained for models with different stellar temperatures and surface gravities. The mean optical spectrum of M33 X-7 (Orosz et al. 2007) is shown as a black solid line. The deep cores in some of the Balmer lines are likely due to an over-subtraction of nebular emission (see supplementary information in Orosz et al. 2007). In contrast to the other He II lines, He II 4686Å appears as a weak emission in the observation. This line might be filled up with nonstellar emission, e.g. from the accretion disk.

Fig. A.2. Same as figure A.1, but zooming on to Hδ comparing models with different surface gravities.

Fig. A.3. Same as figure A.1, but zooming on to multiple He I and He II lines.

All Tables

All Figures

	Fig. 1. Sketch of the binary orbit in M33 X-7. The observer is located along the horizontal axis on the right. The thick lines span the range of orbital phases covered by the HST observations. The orbital separation is taken as ≈1.75 stellar radius as obtained in this work (see Table 4).
In the text

	Fig. 2. HST/COS aperture around M33 X-7. The background image is an RGB composite of HST images in F336W, F475W, and F814W bands. A nearby cluster falls within the FOV of the COS spectrograph.
In the text

	Fig. 3. XMM-Newton EPIC PN, MOS1, and MOS2 spectra of M33 X-7 (black, red, and green curves, respectively) obtained at phase ϕ ≈ 0.47. The error bars correspond to 3σ; the best-fit combined thermal plasma (apec) and diskpbb model is shown by solid lines. The model parameters are shown in Table 3. Lower panel: residuals between the data and the best-fit model.
In the text

	Fig. 4. Same as in Fig. 3, but for phase ϕ ≈ 0.33.
In the text

	Fig. 5. Spectral energy distribution of M33 X-7. The model SED (red solid line) is adjusted to fit the photometry in UV and optical bands (blue boxes).
In the text

	Fig. 6. Spectral fit for M33 X-7 UV spectra at the X-ray eclipse. Normalized HST spectra (blue solid line) is overplotted with the PoWR model (red dashed line). The atmospheric parameters and abundances of this best-fit model are given in Tables 4.
In the text

	Fig. 7. Iron forest in the FUV range observed in the M33 X-7 spectra (black line). Over plotted in blue, red, and green are synthetic spectra obtained for models with the same parameters but with different iron abundances, ranging to 1, 1/2, and 1/5 Fe⊙ respectively.
In the text

Fig. 8. Mass plane diagram for components in the M33 X-7. Masses derived in this work are marked in red whereas gray points represent previously determined values from Orosz et al. (2007). The uncertainty in masses is shown as error bars. The solid lines represent loci of constant orbital inclination angle of 65° derived in this work. The dashed line corresponds to the inclination angle of 75° derived by Orosz et al. (2007).

In the text

	Fig. 9. Comparison of HST/COS spectra of M33 X-7 taken during different orbital phases. Line colors indicate different orbital phases as given in the legend. Upper panel: the full spectra range and lower panels focus on wind line profiles. The observed line variations are due to H–M effect (see Sect. 5).
In the text

	Fig. 10. Comparison of observed Si IV (left panel) and C IV (right panel) line profiles (blue solid) at eclipse to the models calculated with different wind parameters and X-rays. Various colors represent changes in the adopted model parameters as given in the legend. Please see text for more details.
In the text

	Fig. 11. Same as Fig. 10, but at inferior conjunction
In the text

	Fig. 12. Comparison of model SEDs (left upper panel) and emerging UV spectra (lower panel) for different X-ray illuminations. The model SEDs are compared with XMM-Newton and HST spectra taken at phase ϕ = 0.53. A zoom on C IV and Si IV line profiles are displayed in the upper panels (middle and right).
In the text

	Fig. 13. Comparison of total wind acceleration normalized to the gravitational acceleration for models including different amounts of X-rays. The location of the sonic and the BH radius are marked with dashed and solid gray lines.
In the text

Fig. 14. Contributions to the radiative acceleration from the different ions for models with no X-ray (top) and including Lx ∼ 1036 erg s−1 (bottom). Symbols and colors representing different elements and their ionization state is shown in the figure. The plot is limited to the ions that contribute at least ∼1% to the radiative acceleration. Line acceleration due to scattering by free electrons (magenta solid curve) and gas pressure (orange dashed curve) and total radiative acceleration (red dashed curve) are plotted for comparison. The radius at which the BH is located is marked by a gray solid line and the location of the sonic point marked by gray dotted line.

In the text

	Fig. 15. Same as Fig. 14, but for model including Lx ∼ 1038 erg s−1.
In the text

	Fig. 16. Same as Fig. 6, but spectra taken at inferior conjunction (ϕ ≈ 0.53) where BH is in front of the donor star.
In the text

	Fig. 17. Same as Fig. 6, but spectra taken at Quadrature (ϕ ≈ 0.33).
In the text

	Fig. 18. Velocity field for M33 X-7 donor adopting β = 0.8. The location of BH and the corresponding wind velocity and orbital velocity are marked.
In the text

	Fig. 19. Spectral energy distribution for models with different X-ray luminosities. The corresponding number of He II ionizing photons are displayed in the legend.
In the text

	Fig. 20. Comparison of UV spectra of M33-X7 around He II 1640 Å taken during different orbital phases. Line colors indicate different orbital phases as given in Fig. 9. The dashed black line represent our best fit model including Lx = 1038 erg s−1, corresponding to ϕ = 0.53.
In the text

	Fig. 21. Best matching MESA evolutionary tracks reproducing the HRD position of O-star donor in M33 X-7. The current position and errors of the O star are taken from the spectroscopic analysis. Initial masses and periods are labeled on the plot.
In the text

	Fig. 22. Illustration of possible evolutionary scenario for M33 X-7 starting from the ZAMS massive binary to HMXB stage followed by common envelope merging.
In the text

	Fig. 23. Similar to Fig. 23, but an alternative evolutionary track assuming no tidal synchronization.
In the text

Fig. A.1. Comparison of synthetic spectra in the blue optical range obtained for models with different stellar temperatures and surface gravities. The mean optical spectrum of M33 X-7 (Orosz et al. 2007) is shown as a black solid line. The deep cores in some of the Balmer lines are likely due to an over-subtraction of nebular emission (see supplementary information in Orosz et al. 2007). In contrast to the other He II lines, He II 4686Å appears as a weak emission in the observation. This line might be filled up with nonstellar emission, e.g. from the accretion disk.

In the text

	Fig. A.2. Same as figure A.1, but zooming on to Hδ comparing models with different surface gravities.
In the text

	Fig. A.3. Same as figure A.1, but zooming on to multiple He I and He II lines.
In the text