© The Authors 2025

1 Introduction

Cygnus X-1 is one of the most well-studied high-mass X-ray binaries (HMXBs) and the only dynamically confirmed, per-sistent black hole (BH) HMXB in the Milky Way. It consists of a BH in a 5.6-day orbit around an O9.7 Iab supergiant star, HD 226868 (Walborn 1973; Brocksopp et al. 1999). Recent studies have yielded more precise orbital and distance measurements of the system, suggesting a BH mass of M_• ≈21 M _⊙, thereby establishing Cyg X-1 as the most massive BH identified in an X-ray binary (e.g., Miller-Jones et al. 2021).

To derive the BH mass in an X-ray binary, it is crucial to accurately estimate the mass of the optical companion. This requires a detailed multi-wavelength spectroscopic analysis of the donor star using stellar atmosphere models. Although Cyg X-1 has been extensively studied in X-rays, there are relatively few studies dedicated to the optical and UV spectral analysis of the donor star to derive its fundamental parameters. Most of these studies have employed simplified approaches, such as models containing only hydrogen and helium (Herrero et al. 1995), plane-parallel models (Caballero-Nieves et al. 2009; Shimanskii et al. 2012), and local thermodynamic equilibrium (LTE) models (Orosz et al. 2009; Miller-Jones et al. 2021), without considering any effects of stellar wind. However, massive O stars possess strong stellar winds that significantly affect their spectral features. When simplified approaches are used in spectral analysis, they directly influence the derived stellar parameters, such as surface gravities, and hence directly affect the mass estimation. Accurate BH mass estimates in Cyg X-1 are essential for constraining massive binary evolution and understanding the formation of massive BHs at solar metallicity. In this work, we present a detailed spectral analysis of the donor star in Cyg X-1, considering stellar wind, to derive more accurate donor and BH mass estimates.

Cyg X-1 is a prolific source of X-rays, which are generated as the BH accretes mass from the stellar wind of its companion star. This accretion process results in complex interactions between the stellar wind and the BH’s gravitational field, leading to various observable phenomena across the electromagnetic spectrum. Our current understanding of the wind of massive stars in BH binaries is mostly based on the Cyg X-1 system. However, the UV spectra of this prototypical system has not been analyzed with the current generation of expanding stellar atmosphere codes, and previous studies are limited to the Sobolev with exact integration (SEI) method (Gies et al. 2008; Vrtilek et al. 2008). Wind mass loss is a key factor in stellar evolution, and hence, a precise measurement of the donor star’s mass-loss rate is essential for understanding the evolutionary pathways of these systems. The presence of wind clumping adds complexity to these measurements. Furthermore, in HMXBs, the mass-loss rate and wind structure of the donor star significantly influences the accretion flow onto the compact object. Cyg X-1 exhibits distinct X-ray states, transitioning between the “low-hard” and “high-soft” states (Liang & Nolan 1984), which are defined by their hard and soft X-ray spectra, respectively. While the trigger for these transitions is still debated, changes in the wind accretion rate are a potential contributing factor (e.g., Gies et al. 2003). Cyg X-1, therefore, provides an important testbed for investigating these processes. Recent work on the similar HMXB system M33 X-7, using phase-resolved combined UV and optical spectroscopy coupled with sophisticated atmosphere models, has demonstrated the potential for significant revisions in derived system parameters and wind properties (Ramachandran et al. 2022). This highlights the potential for new insights and advancements in our understanding of Cyg X-1 through a similarly detailed analysis.

Therefore, in this paper, we present a combined UV and optical spectroscopic analysis using detailed atmosphere models to address the current gap in our understanding of Cyg X-1. The distance measurement, orbital parameters, and X-ray luminosities of the system are adopted from the literature as listed in Table 1. The archival multiwavelength spectroscopic observations are described in Sect. 2. The details of the atmosphere modeling are given in Sect 3. The results of the spectral analysis, including fundamental properties, abundances, mass of the components, and stellar wind parameters, are presented in Sect. 4. Evolutionary implications of the system are discussed in Sect. 5, and a summary of this study is provided in Sect. 6.

2 Spectroscopy

Archival high-resolution UV and optical spectra, spanning various orbital phases and X-ray states, are used in this study. Specifically, observations near phases 0 and 0.5 were selected to sample distinct points in the binary orbit and provide a more detailed picture of the system. Two sets of optical spectra were used in this study. The low-hard state sepctra, was retrieved from the IACOB database¹ and obtained with the high-resolution FIbre-fed Echelle Spectrograph (FIES) on the Nordic Optical Telescope (NOT). We also utilized high-soft state optical spectra obtained with the High-Efficiency and high-Resolution Mercator Echelle Spectrograph (HERMES) on the 1.2 m Flemish Mercator Telescope (Raskin et al. 2011). The reduction of the HERMES spectra is described in Mahy et al. (2022).

For the UV, we primarily utilize spectra acquired with the Hubble Space Telescope (HST) using the Space Telescope Imaging Spectrograph (STIS). These HST UV observations were taken when Cyg X-1 was in the high-soft state (Gies et al. 2008). To gain a preliminary understanding of wind at the low-hard state (for which HST UV data do not exist), we used the low-resolution International Ultraviolet Explorer (IUE) far ultraviolet (FUV) spectra at phase 0. The details of the spectra used in this work are tabulated in Table 2. Additionally, we used archival FUV and near ultraviolet (NUV) IUE spectra to construct the spectral energy distribution (SED). In addition to the spectra, we used various photometric data (from UV to infrared) to build the SED, including U (Bregman et al. 1973); B, and V magnitudes (Zacharias et al. 2013); the Gaia G magnitude (Gaia Collaboration 2022), and JHK and WISE magnitudes (Cutri et al. 2003, 2012).

3 Stellar atmosphere modeling

Stellar atmosphere models have undergone significant advancements in recent decades. Despite being the subject of numerous studies, the current generation of stellar atmosphere models has never been used to conduct a detailed spectral analysis of Cyg X-1. The study conducted by Herrero et al. (1995) utilized atmospheric models consisting solely of hydrogen and helium, without taking into account the influence of metal line blanketing. Caballero-Nieves et al. (2009) employed plane-parallel models from the TLUSTY grid to study optical and UV spectra. Similarly, Shimanskii et al. (2012) also used hydrostatic, plane-parallel stellar atmosphere models to analyze the optical spectra. However, the plane-parallel approximation is suitable only for stars that do not have winds and where the spectrum is solely formed below the sonic point in a quasi-hydrostatic photosphere. Given that the donor in Cyg X-1 is an O supergiant with substantial stellar wind, a model atmosphere inherently accounting for the wind is highly recommended, as the deduced stellar parameters could otherwise be wrong due to wind emission filling up the absorption profiles. This can have an impact on the derived stellar parameters, such as log g. To get the parameters of the donor, Orosz et al. (2009) employed LTE models using the generalized model stellar atmosphere code PHOENIX without considering any effect of stellar winds. Miller-Jones et al. (2021) followed the same method used in Orosz et al. (2009), albeit with a revised distance and optimization of parameters. On the other hand, studies such as Gies et al. (2008) and Vrtilek et al. (2008) which focused on studying the wind of Cyg X-1, are limited to the SEI method to calculate UV wind line profiles.

Sometimes, HD 226868 is included in samples of stars analyzed by modern non-LTE models, which account for stellar winds. Mahy et al. (2022) carried out an optical spectroscopic analysis of a sample of Galactic SB1 systems, including Cyg X-1, using CMFGEN models. HD 226868 is also included in the sample of the automatized analysis by de Burgos et al. (2024) where optical spectra are compared to a large set of FASTWIND models with fixed abundances.

Therefore, in this work, we address these limitations by performing a combined UV and optical spectral analysis of Cyg X-1 at two key orbital phases (ϕ ≈ 0 and 0.5) using the state-of-the-art, non-LTE Potsdam Wolf–Rayet (PoWR) atmosphere model². 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 reliable 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).

A PoWR model is specified by the stellar temperature T_* ,the bolometric luminosity L, the surface gravity g_* , the mass-loss rate $\dot{M}$, the wind terminal velocity ν_∞, 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}_{*}^2{T}_{*}^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 = 100 R _* .

In the main iteration, thermal broadening and turbulence are approximately accounted for by assuming Gaussian line profiles with a Doppler width of 30 km s⁻¹. For the calculation of the emergent spectrum in the formal integral, the Doppler velocity is split into the depth-dependent thermal velocity and a microturbulence velocity ξ(r). We assume ξ(r) = max(ξ_min, 0.3 (r)) for the models, where the photospheric microturbulent velocity ξ_min is varied during the spectral analysis. Pressure broadening is also taken into account.

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 (r) is prescribed assuming a so-called β-law ν (r) = ν_∞ (1 ™ R_* /r)^β (Castor et al. 1975), where ν_∞ is the terminal wind velocity. In this work, we explore different values for β.

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 inter-clump 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 this study, we account for a depth-dependent clumping assuming that clumping begins at the sonic point, increases outward, and reaches a density contrast of D at a radius of R_D. We adjusted the values of D and R_D during the analysis.

The models account for complex atomic data of H, He, C, N, O, Mg, Si, P, S, and the iron-group elements. A list of the included ions as well as the number of accounted levels and line transitions is given in appendix Table A.1. Initially, chemical abundances are adopted from Asplund et al. (2009) and then adjusted to get a better match to the observed spectra. The iron group elements (Sc, Ti, V, Cr, Mn, Fe, Co, and Ni), with their numerous levels and line transitions, were treated using a super-level approach (Gräfener et al. 2002). This method combines levels and transitions into superlevels with pre-calculated transition cross-sections, assuming solar abundance ratios relative to iron (see appendix Table A.2).

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 non-shocked plasma. We set the onset radius to the orbital distance Article number, page 3 of 23 between the donor star and the BH. The X-ray temperature and filling factor are adjusted such that the wavelength-integrated X-ray flux from the observations is reproduced by the model. We follow a similar procedure as Hainich et al. (2020), adjusting input parameters to simulate different X-ray spectral states.

4 Spectral analysis and results

We performed the spectral analysis by iteratively fitting the observed spectra with PoWR models. To do that, we begin the analysis with Milky Way OB-star grid models calculated for solar metallicity with the stellar temperature T_* and the surface gravity log g_* as main parameters (see Hainich et al. 2019, for more details). Subsequently, we calculate tailored models with minor adjustments in the stellar parameters and surface abundances to match the observed spectra. This is followed by varying wind parameters in the models to reproduce the UV wind line profiles. More details about deriving the stellar and wind parameters as well as chemical abundances, are described in the following subsections. The final parameters derived in this work are given in Table 3.

4.1 Stellar parameters

The optical spectra obtained at various orbital phases and X-ray spectral states are found to be similar, except for changes in radial velocity and emission components in the Balmer lines, particularly in Hα. The Hα variability is depicted in Fig. 1, displaying substantial variations in both the absorption and emission components. Additional emission components from the focused stellar wind can be observed during low-hard spectral states (Yan et al. 2008). Therefore, we excluded Hα from the spectral analysis. Other Balmer lines also exhibit slight variations across different spectral states, becoming deeper in the high-soft state. On the other hand, all the helium and metal absorption lines and wings of the Balmer lines, are not variable (aside from radial velocity shifts). Hence we focused on these lines to derive stellar parameters and abundances of the donor.

The main diagnostic we use to constrain the temperature of the donor is the He and Si ionization balance based on He I/He II and Si III/Si IV line ratios. The following set of optical lines are used for this purpose: He I 3820, He I 3866, He I 4026, He I 4121, He I 4144, He II 4200, He I 4388, He I 4471, He II 4542, He I 4713, He I 4920, He I 5016, He II 5412, He I 5876, He II 6670, He I 7065, Si IV 4089, Si IV 4116, Si III 4553, Si III 4568, Si III 4575, Si III 5740. The pressure-broadened wings of the Balmer lines are the primary diagnostics for surface gravity. We considered Hϵ, Hδ, and Hγ since they are less impacted by wind emissions. Since the ionization balance also reacts to gravity, we simultaneously re-adjusted T_* and log g_* to achieve a good fit to the observed spectra. A comparison of how models with different temperatures and surface gravities impact the diagnostic line profiles is illustrated in Figs. 2 and 3. The optical spectra are best reproduced by a model with T_* = 29 kK and log g_* = 3.2, which aligns with the spectral classification as a late O-type supergiant. We noted that our final model better reproduces the depth of Balmer lines at high-soft state (Fig. B.2) compared to low-hard state (Fig. B.4).

Additionally, we varied the micro-turbulence velocity ξ in our models in the range of 10–30 km s⁻¹ to fit the overall line profiles of metal lines. The determination of ξ is crucial, so it was constrained along with T_eff and log g during our iterative procedure. An inappropriately chosen microturbulence can lead to substantial shifts in temperature and surface gravity (see, e.g., Nieva & Przybilla 2010). Consequently, several elements were analyzed to derive ξ and its impact on line profiles of metal lines are shown in appendix Fig. B.1. We found that the metal lines in the optical and UV spectra are reproduced well with ξ = 22 km s⁻¹. With the exception of the automatized solution by de Burgos et al. (2024), who list ξ = 30 km s⁻¹, this is much higher than in any other previous spectral analysis of Cyg X-1: TLUSTY grid models have a fixed ξ = 2 km s⁻¹ for B stars and ξ = 10 km s⁻¹ for O stars; Orosz et al. (2009) and Miller-Jones et al. (2021) adopted ξ = 2 km s⁻¹ in their grids; Shimanskii et al. (2012) has used a fixed value of ξ = 5 km s⁻¹; and Mahy et al. (2022) used a fixed value of ξ = 10 km s⁻¹ in their CMF-GEN grids. Fixing the micro-turbulence to a lower value will have an impact on the overall spectral fit and the derived stellar parameters and chemical abundances. The fine-tuning of microturbulence is performed along with stellar parameters, and chemical abundances in an iterative manner in this work. Our derived value is in line with the median value of 20 km s⁻¹ obtained by de Burgos et al. (2023) for O9 – B0.5 supergiants in the Galaxy.

For the spectral fit, we adopted the projected rotational velocity ( sin i) and macro-turbulent ( ν_mac) velocity from Simón-Díaz et al. (2017). Subsequently, these velocities, along with instrumental broadening, were used to convolve the model spectra to match the observations. Overall, our best-fit model agrees well with the observed UV and optical spectra.

The luminosity and color excess E_{B ™ V} were determined by fitting the model SED to the photometry and flux-calibrated UV spectra (see Fig. 4). We adopted a standard reddening law by Fitzpatrick (1999). Subsequently, the model flux is scaled with the adopted distance modulus of 11.73 mag (based on Miller-Jones et al. 2021).

Our final best-fit model successfully reproduces the observed spectra across both X-ray states and orbital phases. The spectral fits that correspond to high-soft state are in Figures B.4 (phase 0) and B.5 (phase 0.5), while the same for low-hard state are shown in Figures B.2 (phase 0) and B.3 (phase 0.5), respectively. Despite small variations in Balmer line depths, our model provides a good fit to all observations. Table 3 lists the best-fit stellar parameters and abundances derived in this work for the Cyg X-1 donor star. Our derived stellar parameters are generally consistent with previous literature estimates (Table 4), though some variations exist.

Fig. 1 Hα line variability observed at different orbital phases and X-ray spectral state of Cyg X-1 as described in the legend.

4.2 Chemical abundances

We adjusted the abundances of individual elements to reproduce the observed strength of their respective optical absorption lines. The full list of lines that were considered for the abundance determination are

• Helium: He I 3820, He I 3866, He I 4026, He I 4121, He I 4144, He II 4200, He I 4388, He I 4471, He II 4542, He I 4713, He I 4920, He I 5016, He II 5412, He I 5876, He II 6670, He I 7065.

• Carbon: C III 1247, C II 4267, C III 4647–4650, C III 5696, C IV 5801, C IV 5812.

• Nitrogen: N IV 1483–1487, N II 3995, N III 4097, N II 4237– 4242, N III 4512–4518, N III 4634–4641, N II 5667–5676.

• Oxygen: O IV 1339–1343, O III 1409–1412, O III 3963, O II 4070–4076, O II 4185–4190, O II 4417, O II 4649, O II 4676, O III 5592.

• Silicon: Si III 1294–1303, Si IV 4089, Si IV 4116, Si III 4553, Si III 4568, Si III 4575, Si III 5740

• Magnesium: Mg II 4481.

Individual elemental abundances by mass and by number are given in Tables 3 and 5 respectively. Compared to solar and average abundances of nearby B stars, we found evidence for nitrogen enrichment (by a factor of ∼ 5) and oxygen deficiency in the Cyg X-1 donor, while carbon remains the same (see relative abundances in Table 6). This nitrogen overabundance is consistent with that reported by Caballero-Nieves et al. (2009). We tested a range of hydrogen mass fractions from 0.4 to 0.73 in the study and found that models with X_H = 0.55 reproduce the intensity of hydrogen and helium lines in the observed spectra most accurately. This provides evidence for He enrichment at the surface of the Cyg X-1 donor. Mahy et al. (2022) and Shimanskii et al. (2012) also reported enrichments in He and N in their analysis.

The surface abundances of massive stars evolve due to mass loss, which exposes deeper layers, and mixing processes, which alter internal compositions. In the case of the Cyg X-1 donor, both effects likely contributed to the observed He and N enrichment. Adopting the inclination of the orbit for the donor’s rotation axis, νsin i implies a rotational velocity at the equator ≈ 240 km s⁻¹, suggesting the donor star is rotating at ≈ 0.6 ν_crit. Assuming tidal synchronization would suggest a higher inclination of i ≈ 33.7 (see Sect. 4.3) and hence ν_rot/ ν_crit ≈ 0.5. Such rapid rotation can lead to enhanced mixing. Additionally, interaction with the BH companion may have partially stripped the outer layers of the donor, further contributing to the enrichment.

Although we found evidence for CNO processed material on the surface of the donor, the overall sum of CNO abundance is in agreement with solar value. The chemical abundances in Cyg X-1 donor relative to solar are illustrated in Fig. 6. Notably, α-elements like Si and Mg are overabundant, by factors of 1.8 and 1.37 respectively (relative to solar abundances from Asplund et al. 2009). Using more recent solar abundance estimates (Magg et al. 2022) these factors slightly change to 1.5 for both Si and Mg. These α-element abundances are also significantly higher than those observed in B stars in the solar neighborhood.

Fig. 2 Comparison of observed optical spectra (blue) to synthetic spectra obtained for models with different stellar temperatures and surface gravities (see legend). The plot illustrates the effect of temperature changes on Si III to Si IV and He I to He II line ratios.

4.2.1 Fe abundance determination from FUV iron forest

The forest of photospheric Fe absorption lines in the FUV spectra of hot massive stars provides a direct method to constrain their metal content. Previous studies have assumed a solar-like iron abundance for the system (X_Fe = 1.2 × 10⁻³). In our models, we varied the Fe mass fractions within the range of 1.1– 2.4 × 10⁻³ (iron group mass fraction X_G = 1.2–2.6 × 10⁻³) and employed Bayesian analysis to determine the best-fit values (see Schösser et al. 2025, for details on the method). We focused on two regions of the FUV spectra, 1420–1520 Å and 1560–1700 Å which are dominated by Fe V and Fe IV lines, respectively.

Our results indicate that the iron forest in Cyg X-1 UV spectra is best reproduced with models assuming a slightly super-solar Fe abundance of 1.33 Z_⊙ (X_Fe = 1.6 × 10⁻³) as illustrated in Fig. 5. The derived Fe abundance is higher than the average Fe abundance of B stars in the solar neighborhood (Nieva & Przybilla 2012). In these models, we used a microturbulence value (ξ) of 22 km s⁻¹, derived from the optical metal lines. However, these iron line features are sensitive to the adopted ξ in the model (see, e.g., Hillier et al. 2003). To determine the possible uncertainty range in the Fe abundance, we varied ξ between 17 and 27 km s⁻¹ along with various Fe mass fractions. The final derived Fe abundance is ${1.33}_{-0.15}^{+0.19}{Z}_{\odot }$.

Several studies have reported exceptionally high Fe abundances in Cyg X-1, ∼ 2.0–5 times solar using reflection fits to X-ray spectra (Tomsick et al. 2014; Walton et al. 2016; Parker et al. 2015; Basak et al. 2017). Similarly, analysis of the Fe Kα emission line indicates Fe enrichment, with estimates of ∼ 1.6–6 times solar (Duro et al. 2011, 2016). These results suggest a substantial enrichment of Fe in the accretion disk or the surrounding material. While these value are significantly higher that our estimated Fe abundance, they are consistent with the general trend of Fe enrichment. The discrepancy may be attributed to the sensitivity to different ionization states of Fe in X-ray and UV, or the underlying physical models used in each analysis.

Fig. 3 Same as Fig. 2, but illustrating the impact of log g on wings of Balmer lines.

Fig. 4 Spectral energy distribution for Cyg X-1. The model SED (red solid line) is adjusted to fit the flux-calibrated HST spectra in the FUV (blue line), IUE spectra in the NUV (black line), and photometry from UV infrared bands (blue boxes).

Fig. 5 Iron forest in the FUV region. The blue line represents the observed HST spectra, while the best-fit model, obtained using Bayesian analysis, is shown in red for comparison.

Fig. 6 Elemental abundances in Cyg X-1 donor O star relative to hydrogen (upper panel) and iron (lower panel). The abundance ratios are relative to solar and in logarithmic scale. For comparison, the average chemical abundance of B stars in the Cyg OB3 (Daflon et al. 2001) and B stars in the solar neighborhood (Nieva & Przybilla 2012) are marked.

4.2.2 Initial supersolar composition or recent enrichment by supernova

In Cyg X-1, we found that the Fe, Si, and Mg abundances are approximately 1.3–1.8 times higher (by mass) than the corresponding solar values (0.2–0.3 dex higher by number, see Fig. 6). This exceeds the typical metal content of H II regions located approximately 7.7 kiloparsecs from the center of the Galaxy. The Cyg X-1 donor may either have been polluted by material from a supernova explosion, particularly from its companion star, or it may have formed from interstellar medium (ISM) material which was already enriched.

Cyg X-1 is located at the boundary of the Cyg OB3 associa-tion. A recent study by Rao et al. (2020) using Gaia astrometric data and radial velocities supports Cyg OB3 as the parent association of Cyg X-1, as its kinematic characteristics align with most stars in Cyg OB3. This is also consistent with previous works by Mirabel & Rodrigues (2003) and Mirabel (2017b,a), indicating that Cyg X-1 likely formed in situ within the association.

Daflon et al. (2001) measured chemical abundances of slowly rotating B stars in OB associations, including Cyg OB3. The average abundance of B stars in Cyg OB3 is listed in Table 5 for comparison. The Fe and Si abundances they reported are rather lower than solar values, while Mg abundance is comparable. However, the Fe and Mg abundances provided by Daflon et al. (2001) are based on analyses using LTE models and hence have higher uncertainties. A similar under-abundance of elements relative to solar is also reported for stars in the nearby association Cyg OB7 by the same study. Berlanas et al. (2018) conducted a study in Cyg OB2 and did not observe any significant rise in [O/H] and [Si/H] levels, leading to the conclusion that there is no evidence of chemical enrichment. While non-LTE Fe determinations of hot stars in the region and in particular of Cyg OB3 are lacking, taking the existing literature at face value makes it puzzling why only the Cyg X-1 system shows super-solar Fe, Si, and Mg abundances compared to the neighborhood.

Certain studies have identified significant variability in the overall metallicity of stars within the proximate Galactic interstellar medium. For example, De Cia et al. (2021) identified substantial disparities in gas-phase metallicity, with variations reaching up to tenfold, in the neighborhood of the Sun. This was established by examining neutral clouds in the line of sight of 25 O- or B-type stars located in or near the Galactic plane and within 3 kpc. In contrast, research by Arellano-Córdova et al. (2020, 2021) indicated that the gas in H II locations is uniformly distributed at a specific Galactocentric distance, based on chemical abundance studies of B-type stars, classical Cepheids, and young clusters. Esteban et al. (2022) also suggest that the chemical composition of the interstellar medium in the solar vicinity is very homogeneous.

The hypothesis that the BH progenitor polluted the donor star’s surface via a supernova explosion has been previously considered. Shimanskii et al. (2012) reported higher Ne and Si abundances relative to solar in the Cyg X-1 donor, which they suggested as evidence for a supernova event. Our analysis of the surface abundances in Cyg X-1 reveals enrichment in Fe, Si, and Mg, but not in CNO elements. However, this observed abundance pattern presents a challenge to a pure supernova origin. While Rodríguez et al. (2023) estimated a mean Fe yield of core-collapse supernovae (CCSNe) to be ≈ 0.058 M_⊙ , with Type Ic supernovae producing more than twice this amount, this is still significantly lower than the Fe yield of Type Ia supernovae ( ∼ 0.6 ™ 0.7 M_⊙ ; e.g., Howell et al. 2009). Based on the Fe yields of CCSNe, Weinberg et al. (2024) inferred corresponding yields of α-elements such as Mg, Si, and O. Considering their estimates, a 33% increase in Fe abundance would be accompanied by ay a ∼ 30–35% increase in Si and Mg, along with an approximately twofold increase in O abundance (see Table 1 and Figure 3 in Weinberg et al. 2024). Even in alternative scenarios like a failed supernova or a direct collapse to a BH, some O enrichment relative to Fe would still be expected. Therefore, the observed abundance pattern in Cyg X-1 is challenging to explain solely by supernova yields or stellar winds from the BH progenitor.

4.3 Mass of the components

The spectroscopic mass of the donor is calculated from log g_* and $R_{*}(g_{*}=G{M}_{*}{R}_{*}^{-2})$ and the uncertainties in these parameters propagate into the final uncertainty of the donor mass. The reduction in luminosity, temperature, and gravity results in a much lower spectroscopic mass of ≈ 29 M_⊙ for the O star than previously reported in Miller-Jones et al. (2021). Considering the orbital parameters in Table 1, we estimate the BH mass, the orbital radius, and the Roche lobe radius of the donor. The newly derived BH mass is ≈ 17.5 M_⊙ , which is lower than the previous estimate of ≈ 21 M_⊙ by Miller-Jones et al. (2021).

The inclination angle of the binary system Cyg X-1 is a crucial parameter affecting the BH mass estimate as demonstrated in the Fig. 7. While previous optical measurements indicate an orbital inclination of i = 27.5^° (Orosz et al. 2011; Miller-Jones et al. 2021), X-ray reflection and polarization studies suggest a significantly higher disk inclination, ranging from 35^° to 45^° (e.g., Krawczynski & Beheshtipour 2022; Zdziarski et al. 2024; Walton et al. 2016; Tomsick et al. 2014). This discrepancy is attributed to a possible disk warp or misalignment between the binary and inner disk axes. If the orbital inclination is higher, then it would lead to significantly reduced BH masses (for instance M_BH ≈ 11.5 M_⊙ at i ∼ 40^°).

While detailed light curve modeling, which could provide further constraints on the system’s geometry and BH mass, is beyond the scope of this paper, we estimate the inclination angle i using the assumption of tidal synchronization. We assume that the primary star’s rotational period is synchronized with the binary orbital period, a plausible scenario given the system’s short orbital period (see, e.g., Stoyanov & Zamanov 2009). The relationship between i, the observed ν sin i, the primary star radius R₁, and the orbital period P is given by

$i=\arcsin (\frac{v\sin i\cdot P}{2\pi R_1}).$(1)

Using ν sin i of 112 km s⁻¹ (see Table 3) in Eq. (1) yields an inclination of i ∼ 33.7^° and BH mass estimate of M_BH ≈ 13.8 M _⊙ . The correlation between the estimated masses of the donor star and the BH, along with their respective uncertainties, is illustrated in Fig. 7.

To further investigate the uncertainty in the BH’s mass, we conducted a Monte Carlo simulation by sampling the O-star mass within the 2σ range. This wider range of donor masses allowed us to investigate the impact of more extreme M_donor values and to ensure the robustness of our results against potential outliers or underestimated systematic errors. Given that the uncertainties in donor mass are asymmetric, we utilized a truncated normal distribution to accurately represent its probability distribution. The simulation adopts a uniform probability distribution of inclination in the range i ∼ 25^°–39^°. This range was estimated based on the assumption of synchronous rotation and reflects the spread of νsin i values in the literature ( νsin i ∼ 90–120 km/s). The resulting M_BH distribution (adopting a mass function as in Table 1), is shown in Figure 8. Applying kernel density estimation (KDE) to characterize this distribution reveals a peak at ≈14.2 M_⊙ with a 1σ range of masses 12.7 to 17.8 M_⊙. The inclination-driven uncertainty, as shown in the M_BH distribution, suggests that the 17.5 M_⊙ estimate (corresponds to i = 27.5^°), likely represents an upper limit. Notably, future astrometry data would provide independent and robust constraints on inclination and mass ratio, significantly reducing uncertainties in BH mass determination.

For M_BH ≈ 17.5 M_⊙ , the orbital separation is 47.8 R_⊙ while the radius of the donor star is R_* ≈22.3 R_⊙ . Eggleton (1983) defines R_eq as an equivalent radius that corresponds to the volume of the Roche lobe when approximated as a sphere. In this approximation, R_eq in Cyg X-1 is roughly similar to the donor radius with a Roche lobe filling factor of $f_{\text{Req}}=R_{*}/R_{\text{eq}}\approx {1.1}_{-0.2}^{+0.15}$. However, these radii are considerably smaller than the distance to the inner Lagrangian point (R_L1).

As further discussed in Sect. 4.5, the classical Roche lobe geometry of the donor star is modified in the presence of a radiation-driven stellar wind (see also Dermine et al. 2009). To visualize this scenario, we employed the formalism from Leahy & Leahy (2015) and constructed a 2D representation of the Roche lobe geometry in Cyg X-1, as shown in Fig. 9. We calculated various radii defining the Roche lobe’s shape: the distance to the inner Lagrangian point R_L1, the radii along the line connecting the stars (R_fr, R_bk), and the radii perpendicular to this line and the orbital plane (R_y, R_z). These radii represent distances to the Roche equipotential surface at specific points, allowing us to map the shape of the Roche lobe. This plot, generated using a mass ratio of q = M_BH/M_donor = 0.603 17.5/29, illust the relative sizes and positions of the donor star, the Roche lobe, and the position of the BH. This analysis yields $f_{\text{R}1}=R_{*}/R_{\text{L}1}\approx 0.85$. Even with a BH mass of ≈ 13.8 M _⊙ , this value only changes marginally to ≈ 0.835. The presumably still rather spherical donor star ( ν_rot/ ν_crit ≤ 0.6) extends beyond the R_eq and the radii in the y, z, and back directions (R_y, R_z, R_bk), yet remains smaller than the distance to the inner Lagrangian point R_L1 ≡ R_fr. This holds true across the entire BH mass probability distribution shown in Fig. 8 (see the Roche lobe filling factor distribution in various directions in appendix Fig. B.7). Nonetheless, the donar star is still close to a situation where overflowing mass transfer could happen. This result is consistent with prior studies (Gies & Bolton 1986b; Gies et al. 2003; Ziółkowski 2005). Similar behaviorior has been observed in other BH-HMXBs, such as Roche lobe overfilling in M33 X-7 (Ramachandran et al. 2022) and nearly Roche lobe filling in LMC X-1 (Orosz et al. 2009).

Fig. 7 Mass plane diagram for Cyg X-1, showing component masses and their uncertainties. BH masses are estimated for i = 27.5°(blue) and i = 33.7° (orange).The solid black line represents the mass relation for the adopted inclination of i = 27.5° derived from optical observations (Orosz et al. 2011; Miller-Jones et al. 2021). The dotted line corresponds to i = 33.7° derived assuming synchronous rotation. The dashed lines illustrate the variation in the mass relation for alternative inclination angles.

Fig. 8 Probability distribution of BH mass in Cyg X-1. The gray histogram and blue line depict the MBH distribution and its KDE, respectively. The shaded region shows the 1-σ range (16th and 84th percentiles). The dotted green line marks the KDE peak, and the dashed orange and blue lines represent the MBH values at inclinations of 27.5° and 33.7° (see legend).

Fig. 9 Illustration of the Roche lobe geometry for Cyg X-1, ing the O-star donor (light blue circle) and the BH (black do MBH = 17.5 M⊙). The dark blue lines and the mesh depict the ap mate Roche equipotential surface, with arrows indicating the dist RL1, Rbk, and Ry. The axes are scaled in R⊙.

4.4 Wind parameters

Winds of massive stars are parameterized mainly by their mass-loss rate $\dot{M}$ terminal wind velocity ν _∞, velocity law exponent β, and clumping factor D. To determine these parameters using spectral analysis, in this work, we used UV spectra taken at ϕ∼ 0 to yield parameters of the O star wind, while the spectra taken at ϕ ∼ 0.5 are employed to study the wind quenching by X-rays. The derived wind parameters at the high-soft and low-hard states are listed in Table 7.

The main diagnostics used to determine the terminal velocity is the blue edge of prominent P Cygni lines such as N V λλ1238– 1242, Si IV λλ1393–1403 and C IV λλ1548–1551. The maximum blue edge velocity from the unsaturated P-Cygni lines is ν _edge ≈1300 km s⁻¹. A terminal velocity of ν_∞≈ 1200 km s⁻¹ shows better agreement with the observations in the high-soft state as shown in Fig. 10. This is in agreement with Gies et al. (2008), where they used the SEI method to derive the wind velocity. On the other hand, the IUE spectra indicate higher wind velocities at low-hard state (Fig. 11). While the low resolution of the IUE data precludes a precise measurement, the modeling of the C IV profile suggests ν_∞ ≈ 1600 km s⁻¹, and the N V profile indicates an even higher value of ν_∞ ≈ 1800 km s⁻¹. Even the hard-state values differ significantly from the commonly used value of ≈ 2100 km s⁻¹, which was originally adopted in Herrero et al. (1992) based solely on optical spectra. Davis & Hartmann (1983) also reported a value of 2300 km s⁻¹, but this was based on examining IUE line widths alone, without any modeling. With our obtained terminal velocity, the ratio of ν_∞ to the escape velocity of the donor at high-soft state is found to be ν_∞ / _esc ≈ 1.7, lower than the average value for typical O stars (Lamers et al. 1995; Kudritzki & Puls 2000). In contrast, the low-hard state exhibits a ratio of ≲2.6, consistent with typical O star values.

For our spectral analysis, we computed models with log $\dot{M}$ in the range of –6.0 to –7.0. To model the supersonic part of the wind, we tested β values in the range 0.8 to 3. Since there are no simultaneous UV and X-ray observations of the system, we varied L_X ∼ 4–9 × 10³⁶ erg s⁻¹ in the model to match UV spectra taken at ϕ∼ 0 in the high-soft state. The P Cygni line profile shapes are best reproduced with β = 1.7. Our best-fit model for the high-soft state incorporates X-rays with L_X ∼ 7 × 10³⁶ erg s⁻¹ (see top panel in Fig. 10) and yields log $\log \dot{M}=-6.5$ or $\dot{M}\approx 3\times {10}^{-7}{M}_{\odot }{\text{yr}}^{-1}.\log \dot{M}$ between −6.6 and −6.4 reasonably reproduce the observed UV P Cygni profiles when varying L_X ∼ 4 − 9 × 10³⁶ erg s⁻¹. In contrast, to reproduce the IUE spectra at ϕ ∼ 0 and in the low-hard state, models with slightly higher mass-loss rates (log $\dot{M}$ between −6.5 and −6.3) and lower beta (β < 1.5) are required (see Fig. 11). To model the low-hard state, we use lower L_X values in the range of ∼5–9 ×10³⁶ erg s⁻¹.

One noticeable feature in the UV spectra are the weak emission components in Si IV λλ1393–1403 and C IV λλ1548–1551, whereas OB supergiants with similar spectral types typically show Si IV and C IV resonance lines with strong emission in the red part. This is likely due to the asymmetries in the spatial distribution of corresponding ions destroyed by photoionization in the hemisphere towards the X-ray source. At phase 0 (BH behind the donor), the less-disturbed wind along our line of sight contributes to strong blue-shifted absorption in the P Cygni profiles, while X-ray photoionization on the opposite side weakens the emission component (as seen in top panel in Fig. 10). If a focused wind component were present toward the BH, we would expect enhanced emission, particularly at phase 0. This enhancement is not observed; the emission is notably weak. A similar absence of a stronger emission component is also observed in other BH-HMXBs such as M33 X-7 (Ramachandran et al. 2022), while systems with neutron stars such as Vela X-1 or 4U 1700-37 show this emission component (e.g., Sander et al. 2018; Hainich et al. 2020) possible because of the lower X-ray luminosities of these systems.

The Hatchett–McCray effect is clearly observed in Cyg X-1 donor’s UV wind lines. X-ray photoionization strongly affects the wind on the X-ray facing hemisphere (phase 0.5), leading to higher ionization stages and decreased UV line intensity (Figure 10, bottom panel), compared to the less ionized wind observed at phase 0. The dramatic decrease in line strength at phase 0.5, especially at high velocities (Figure 10, upper and lower panels), indicates significant disruption of the donor wind by the accreting black hole and suggests that X-ray ionization extends to the photosphere. This decrease contradicts the expectation of increased absorption from a focused wind when the BH is in the line of sight (e.g., Gies & Bolton 1986a; Sowers et al. 1998; Tarasov et al. 2003; Miller et al. 2005), as discussed by Vrtilek et al. (2008). Furthermore, our final model, using the same wind and stellar parameters (except for an increased L_X ∼2 × 10³⁷ erg s⁻¹ consistent with observations, e.g., Wilms et al. 2006; Sugimoto et al. 2017) successfully reproduces the UV spectra at phase ϕ∼ 0.5 (see bottom panel of Fig. 10). Therefore, UV observations show no evidence of increased mass flow focused towards the BH, further strengthening the argument against a focused wind.

A microturbulent velocity of 22 km s⁻¹ in the photosphere is derived based on the optical lines, and is also in agreement with the UV spectra. Our adopted micro-turbulent velocity ξ in the formal integral grows proportional to the wind velocity. For typical OB star analyses, it is sufficient to put the scale factor such that the maximum is ξ(R_max) = 0.1ν_∞ . However, for Cyg X-1 we measure a higher turbulence velocity, with a ratio ξ(R_max)/ ν_∞ ≈ 0.3, in agreement with the reported values for M33 X-7 (Ramachandran et al. 2022). The wind micro-turbulent velocities derived for these HMXBs are systematically larger compared those found in previous studies of OB stars in the Galaxy (e.g., 0.1ν_∞ by Kudritzki & Puls 2000, 0.14 ν_∞ by Herrero et al. 2001).

The adopted microclumping is depth-dependent. Assuming a (maximum) clumping factor of D = 30 and R_D = 2R_* best reproduce the emission in He II at 4686 Å. Our analysis suggested that the clumping factor reaches this maximum value close to the location of the BH (2.1 R_*).

In this work, we derive donor mass-loss rates of $\dot{M}\approx 3\times {10}^{-7}{M}_{\odot }{\text{yr}}^{-1}$ in the high-soft state while in the low-hard state we measured $\dot{M}\lesssim 5\times {10}^{-7}{M}_{\odot }{\text{yr}}^{-1}$. These significantly lower estimates represent a more detailed and accurate determination of the mass-loss rate compared to previous works. For instance, Gies et al. (2003) derived $\dot{M}\approx 2-2.57\times {10}^{-6}{M}_{\odot }{\text{yr}}^{-1}$ across different X-ray states. However, their analysis was solely based on the Hα equivalent width, neglected the effects of clumping and X-ray ionization, and relied on stellar parameters from Herrero et al. (1995). Using the UV spectra taken during the high-soft state, Vrtilek et al. (2008) found a much higher mass-loss rate, $\dot{M}\approx 4.8\times {10}^{-6}{M}_{\odot }{\text{yr}}^{-1}$ using the SEI method.

Considering theoretical predictions, our derived wind mass-loss rate is a factor of 1.7–2.7 smaller than predicted by the Vink et al. (2000) recipe for the same mass, luminosity, and temperature of the star at solar metallicity. When take into account the derived high iron abundances, the emprically derived $\dot{M}$ is a factor of 2.2–3.5 lower than the “Vink mass-loss rate” at Z = 1.33 Z _⊙ . This discrepancy can be attributed to the effects of wind clumping. Inadequate consideration of wind clumping can result in substantial overestimations of mass-loss rates, with discrepancies ranging from factors of 2–10 (see, e.g., Puls et al. 2008; Fullerton et al. 2006). This aligns with a recent analysis of early-type O supergiants in the Galaxy by Hawcroft et al. (2021), who reported a mass-loss rate reduction by a factor of 3.6 compared to Vink et al. (2000) predictions, along with high clumping factors.

Fig. 10 Comparison of UV P-Cygni profiles observed with HST (blue) at orbital phase ϕ∼ 0 (top panel) and ϕ ∼0.5 (bottom panel) during high-soft state to the best-fit models (red). The models assume the same input parameters except for the incorporated amount of X-rays. The wind parameters used in the best-fit models are log $\log \dot{M}(M_{\odot }{\text{yr}}^{-1})$ and ν∞ =1200 km s−1. For superior conjunction ϕ ∼ 0 we adopted LX ∼ 7 × 1036 erg s−1 in the model, and at inferior conjunction ϕ∼ 0.5 we used higher amount of X-rays LX∼ 2 ×1037 erg s−1. The variability in the wind line profiles shows the impact of X-ray ionization.

Fig. 11 Same as Fig. 10 but IUE observations at orbital phase ϕ ∼ 0 during low-hard state are shown. We adopted a lower X-ray luminosity LX ∼ 3--5 × 1036 erg s−1 in the model. Adopted wind parameters in models are shown in legend.

4.5 Wind accretion versus Roche-lobe overflow

In wind-fed supergiant-HMXBs, X-rays are created by the accretion of wind material from the donor star onto the compact companion. The mass accretion rate via the Bondi-Hoyle formula is

${\dot{M}}_{\text{acc}}=\pi R_{\text{acc}}^2\rho v_{\text{rel}},$(2)

where ρ is the wind density near the BH and R_acc is the accretion radius. In Cyg X-1, the wind velocity near the BH (at 2.1 R_* ) reaches around ≈ 340 km s⁻¹ at high-soft state and up to ≲860 km s⁻¹ at low-hard state. The value of the velocity exponent β derived in this work ranges from 1 to 1.7, defining the wind structure and its influence on the wind velocity determined for the position of the BH. This higher β value indicates a more gradual acceleration of the wind, resulting in lower wind velocities near the BH orbit compared to the typically assumed value of β = 0.8, which would imply a steeper acceleration and higher wind velocities. The mass accretion rate onto the compact object is a sensitive function of the wind velocity law, following the relationship ${\dot{M}}_{\text{acc}}\propto v^{-4}$. This means that even small changes in the wind velocity can lead to significant variations in the mass accretion rate (e.g., Oskinova et al. 2012).

The orbital velocity (2πa/P) near the BH is 270 km s⁻¹. The relative velocity between the accreting BH and the stellar rel wind can be estimated using $v_{\text{rel}}^2=v_{\text{wind}}^2+v_{\text{orbit}}^2\approx 435\hspace{0.17em}\text{km}\hspace{0.17em}{\text{s}}^{-1}$ at high-soft state. This derived ν_rel is only 0.36 ν_∞. In the low-hard state, _rel ≲ 901 km/s or 0.5ν_∞. The density of the wind at the distance of the BH can be approximated to

$\rho =\frac{{\dot{M}}_{\text{donor}}}{4\pi d_{\text{BH}}^2{v}_{\text{wind}}}.$(3)

By inserting the derived wind parameters, we get a density of ρ≈ 4×10⁻¹⁵g cm⁻³ in the high-soft state and ≳1.5×10⁻¹⁵ g/cm³ in the low-hard state. The accretion radius, $R_{\text{acc}}=2G{M}_{\text{BH}}/v_{\text{rel}}^2$, is 2.4×10⁷ km in the high-soft state, yielding a mass accretion rate of ${\dot{M}}_{\text{acc}}=3.3\times {10}^{18}\hspace{0.17em}\text{g}\hspace{0.17em}{\text{s}}^{-1}$. In the low-hard state, the smaller accretion radius, R_acc ≈ 6×10⁶ km, results in a lower mass accretion rate, ${\dot{M}}_{\text{acc}}\approx 4.4\times {10}^{17}\hspace{0.17em}\text{g}\hspace{0.17em}{\text{s}}^{-1}$.

The accretion disk luminosity in BHs arises from material spiraling inward, releasing gravitational energy as it moves from the outermost radius to the innermost stable circular orbit R_BH. The accretion luminosity can be written as

$L_{\text{acc}}=\frac{G{M}_{\text{BH}}{\dot{M}}_{\text{acc}}}{R_{\text{BH}}}=\eta {\dot{M}}_{\text{acc}}.$(4)

Rapidly spinning black holes, like the one in Cyg X-1, have accretion disks with inner edges much closer to the event horizon. This significantly increases the accretion efficiency η and hence the L_acc. For Cyg X-1, estimates of the spin parameter range from a_* 0.9 (Zdziarski et al. 2024) to nearly 1 (Zhao et al. 2021). Our lower BH mass estimate suggests a potentially slightly lower spin (a_* ∼ 0.7–0.8), but high efficiency remains likely. For instance, with η = 0.2 the resulting accretion luminosity would be L_acc ∼ 6×10³⁸ erg s⁻¹ in the high-soft state and L_acc ∼ 4.5 × v 10³⁷ erg s⁻¹ in the low-hard state.

Although the calculated accretion luminosity is below the Eddington limit L_Edd ∼2 × 10³⁹ erg s⁻¹, it significantly exceeds the observed X-ray luminosity – by an order of magnitude in the high-soft state and a factor of 2 in the low-hard state. Such large differences between Bondi–Hoyle accretion luminosity and observed X-ray luminosity have been reported for other BH HMXB systems such as IC 10 X-1 and NGC 300 X-1 (Tutukov & Fedorova 2016). This discrepancy highlights the limitations of the Bondi-Hoyle accretion model, which simplifies the complex physics of accretion. Our analysis also uses a typical β-law assumption for the wind velocity field, potentially underestimating the wind velocity near the BH. Moreover, in the above calculations we used the BH mass close to its maximum value of 17.5 M_⊙ . Variations in accretion efficiency, radiative processes, absorption, scattering, and geometric effects can further contribute to discrepancies between theoretical and observed luminosities.

Our calculations indicate that the observed X-ray luminosity in Cyg X-1 could be fully explained by the accretion of the donor star’s stellar wind. If there would be additional mass flow in the direction of the BH due to the Roche lobe overflow of the donor, this should cause a much higher X-ray luminosity than currently observed. The latter phenomenon is observed in systems like M33 X-7, where the stellar wind of the donor is insufficient to account for the high X-ray luminosity, indicating the need for invoking additional mass flow towards the BH by Roche-lobe overflow or wind-RLOF (Ramachandran et al. 2022).

As we discuss and illustrate in Sect. 4.3, in Cyg X-1 the donor radius matches or even exceeds the equivalent radius of the Roche lobe. However, the conventional Roche lobe definition might not be applicable in a straightforward way to an O-type donor star with radiatively driven stellar wind. While we do not solve the wind hydrodynamics in this work and thus cannot examine the acceleration in the launching regime in more detail, the presence of a radiation-driven wind implies Γ_rad > 1, meaning that a conventional Roche lobe actually does not exist for the donor (see, e.g., Dermine et al. 2009). Even in the case of strong X-ray illumination, our models show Γ_rad ≈1 at the radius corresponding to the orbit of the BH. The acceleration in the outer wind is much higher than this. In addition, we see evidence of wind signatures in the UV spectrum during ϕ ∼0.5 such as C IV and N V resonance lines (see Fig. 10). We thus conclude from our modeling that there is no conventional Roche lobe overflow present in Cyg X-1 which is in line with our finding that the stellar radius does not yet reach the radius of the inner Lagrangian point (cf. Fig. 9).

Although the precise drivers of Cyg X-1’s X-ray state transitions remain unclear, our wind analysis points to a possible connection with wind-accretion feedback. In the high-soft state, the lower wind velocity and mass-loss rate result in a higher mass accretion rate and, consequently, a higher X-ray luminosity. Conversely, in the low-hard state, the significantly higher wind velocity and moderately higher mass-loss rate result in a lower mass accretion rate and, as a result, lower X-ray luminosity. This suggests a complex interplay where changes in the wind affect the accretion flow, which in turn influences the observed state transitions. Further investigation is needed to disentangle this complex relationship.

5 Evolution of the system

We explore the system’s future evolution channel towards binary BHs and investigate its possibility of becoming a gravitational wave source. We use 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) to compute the tracks. Given the substantial uncertainties in the system’s prior evolution (e.g., stable mass transfer, common envelope evolution, supernova kick uncertainties, potential dynamical capture), modeling it is beyond the scope of this work. Our focus is to assess the system’s current evolutionary state – specifically, whether it is wind-fed or currently undergoing Roche lobe overflow. We model the binary evolution, starting from a configuration with a ZAMS O-star donor (initial mass in the range of 29–35 M_⊙ ) and a BH companion (17.4 M_⊙ ) and a 5.5 day initial period. These initial parameters are chosen so that they are close to the estimated masses and period of the system from observations, and the main objective is to match the current state of the system. Evolution is then computed until either core carbon is depleted or until a merger between the donor star and the BH occurred. After core carbon depletion, a star is assumed to directly collapse into a BH without losing mass. Subsequently, we assign the final core mass to that of the compact companion.

Our models are computed with a solar metallicity of Z = 0.0142 (Asplund et al. 2009), where we take the relative metal mass fractions from Grevesse & Sauval (1998). In our calculation, the binary system is assumed to be synchronized with the orbital period at the ZAMS. The evolution of orbital angular momentum takes into account spin-orbit coupling, gravitational wave radiation, and mass loss, as explained by Paxton et al. (2015).

Convection was modeled using the standard mixing-length theory Böhm-Vitense (1958) with a mixing-length parameter α = 1.5, adopting the Ledoux criterion. Convective overshooting is implemented with a step-overshooting scheme, extending the convective core by α_ov = 0.345 pressure scale heights, thereby following Brott et al. (2011). The efficiency parameter α_SC = 1.0 is used to represent semiconvection, following the work of Langer et al. (1983). The effect of the centrifugal force was implemented as in Heger & Langer (2000). Rotation-induced composition-and angular-momentum transfer encompasses the GSF instability, secular and dynamical shear instabilities, and the impacts of Eddington-Sweet circulation.

We use MESA’s “Dutch” scheme for wind predictions, with cool star winds following Nieuwenhuijzen & de Jager (1990), hot-star winds applying Vink et al. (2000, 2001), and Wolf– Rayet stars following Nugis & Lamers (2000). The mass-loss rates in these prescriptions are most likely an overestimation, for example due to their assumptions of homogeneous winds. Our analysis of the donor wind suggests a high clumping factor and lower mass-loss rate than predicted by Vink et al. (2001) for the same stellar parameters. We thus computed models with scaling factors of 1, 0.5, and 0.2 for the wind mass loss in order to approximately counteract a possible overestimation of mass-loss rates.

To find the best matching models for the observed system parameters, we compare the stellar parameters of the O star donor derived from spectral analysis – such as effective temperature, luminosity, and surface gravity – along with the current orbital period of the system, to the evolutionary model tracks. Models starting with an O star donor of 34 M_⊙ and a BH companion of 17.4 M_⊙ with a 5.5-day period successfully reproduce the observed properties of the Cyg X-1 system within uncertainties for all three wind mass-loss scaling factors. The corresponding tracks are shown in the HRD in Fig. 12, and the best-matching parameters from these models are compared in Table 8. The current wind mass-loss rate of the O star donor is best reproduced by the evolutionary model that assumes a mass loss scaling factor of 0.5 (model “M2” in Table 8).

The current properties of the system match with a phase just before the first mass transfer from the donor to the BH via RLOF. At this stage, the donor star is still burning hydrogen in its core, indicating it is on the main sequence or slightly evolved from it. After the mass transfer, the donor’s mass is reduced to approximately 21 M_⊙ and continues hydrogen burning for another ∼1 million year. Despite this reduction in mass and luminosity, the donor retains the appearance of an O supergiant, along with surface He and N enrichment. When hydrogen is depleted in the core and helium burning begins, the next mass transfer phase is initiated. During this phase, the donor’s mass decreases further, eventually falling below 14 M _⊙ . The donor star rapidly contracts after that, becoming a Wolf-Rayet (WR) star.

The stellar wind during the WR phase is significantly influenced by the different mass-loss scaling factors adopted in our evolutionary models. For instance, a mass-loss scaling factor ( f_wind) of 1 and 0.5 leads to more extensive stripping of the stellar layers, driving the star to a more evolved WC phase. In contrast, a model with a lower mass-loss rate f_wind = 0.2 remains in the WN phase. Additionally, the orbit of the binary system widens due to the strong wind mass loss during the WR phase. The final core masses of the donor range from 9.5 to 13.3 M_⊙, and the final orbital period varies from 10.9 to 7.5 days for f_wind = 1 to f_wind = 0.2.

We calculated the accretion luminosity throughout the system’s evolution using Eq. (4), accounting for both Roche lobe overflow and wind mass transfer (Fig. 13). For _wind, we assumed a β law with a standard β-value of 1. The terminal velocity ν _∞ can be estimated from the effective escape velocity

$v_{\text{esc},\Gamma }=\sqrt{\frac{2GM}{R}(1-{\Gamma }_{\text{e}})}.$(5)

Empirical studies of OB stars in the Milky Way suggest ν_∞ / ν_esc,Γ ≃ 2.6 for stars with T_* ≥ 21 kK (Lamers et al. 1995; Kudritzki & Puls 2000). A similar ratio is also found for Galactic WN stars by Niedzielski & Skorzynski (2002). The accretion luminosity is calculated only when an accretion disk formed around the BH, requiring the accreted material to have sufficient angular momentum to exceed that of the innermost stable circular orbit (see equation 10 in Sen et al. 2021) Mass transfer is treated using the ‘roche lobe’ scheme. In this scheme, the mass transfer rate is dynamically adjusted to ensure the donor star remains within its Roche lobe, with the Roche-lobe radius computed as in Eggleton (1983). The observed X-ray luminosity matches the accretion luminosity of the evolutionary model just before the mass transfer phase (yellow region in Fig. 13). Soon the system will undergo RLOF and the luminosity will be increased ≳10³⁹ erg s⁻¹. The donor star continues to fill its Roche lobe until it finishes the second mass transfer phase and contracts towards the WR stage. During this time, the system will appear as an ultra-luminous X-ray source (ULX).

Initially, we assumed a highly non-conservative mass accretion scenario where the BH accretes at the Eddington limit, with the remaining material ejected from the system carrying the specific angular momentum of the accretor. This resulted in limited BH growth and consequently, wide post-mass transfer orbits. As a result, the predicted merger timescale is significantly longer, estimated to be approximately 51 Gyr. Three BH accretion scenarios – Eddington-limited, moderate super-Eddington, and totally conservative accretion – have been compared in a recent work by Xing et al. (2025). Considering their findings, the only scenario that could result in a potential BH merger for Cyg X-1 would be a highly conservative accretion scenario.

Motivated by this, we re-examined our best-fit model (M2, with f_wind = 0.5) using a fully conservative accretion scenario. While this change does not significantly affect the currently observed parameters of Cyg X-1 (which match a pre-Rochelobe overflow state, see Table 8), it drastically alters the future evolutionary outcome. As shown in Fig. 14, the final BH mass increases substantially (from 18.7 to ∼35.8 M_⊙ ) due to the enhanced accretion. Enhanced accretion would also result in much higher X-ray luminosities on the order of ≳10⁴⁰⁻⁴³ erg s⁻¹ during or after the mass transfer phase. The period of the system only slightly increases, and the mass of the donor remains the same. Moreover, Xing et al. (2025) suggest that only a highly conservative accretion onto BHs lead to extreme BH spins a_* ∼ 1. Our fully conservative model suggests that the system would become binary BHs with final masses of 35.8 M_⊙and ∼11.5 M_⊙ and a final period of ∼9.13 days, making it a strong gravitational wave progenitor that will merge in approximately 5 Gyr. It should be noted that this merger timescale will increase if we consider lower values for the current BH mass in the model. Nonetheless, our findings indicate that mass accretion efficiency plays a crucial role in the evolution of binary systems toward becoming merging black holes and potential gravitational wave sources, with Cyg X-1 as a prime example.

Fig. 12 Potential future evolution of the Cyg X-1 donor star. Here, the evolution starts when the primary is already a BH (point mass) of 17.4 M⊙, O star donor of initial mass 34 M⊙, and an initial period of 5.5 days. Three models are displayed here with mass-loss scaling factors of 1, 0.5, and 0.2 to demonstrate the effect of wind mass loss on evolution.

Fig. 13 Left panel: same as Fig. 12 but here the tracks correspond to a wind mass-loss scaling factor of 0.5, and the color scaling shows the accretion luminosity. Right panel: accretion luminosity (top) and Roche-lobe filling factor (bottom) as a function of normalized donor lifetime τ. The gray bar shows the range of observed X-ray luminosity of the system.

Fig. 14 Impact of mass accretion by the BH on the evolution of Cyg X-1.The panels show the evolution of accretion luminosity, orbital period, and component masses as a function of normalized time τ for MESA models with Eddington-limited (blue) and fully conservative (red) mass accretion.

6 Summary

We have presented a comprehensive analysis of Cyg X-1, utilizing archival high-resolution UV and optical spectra in conjunction with sophisticated atmospheric models. Notably, this is the first investigation to simultaneously analyze UV and optical spectra of Cyg X-1 while also incorporating X-rays to constrain the stellar and wind parameters. Our key findings are as follows:

• Our analysis yields notably lower masses for both the donor ( ≈29 M_⊙ ) and the BH. Assuming the inclination from the literature, we obtain a BH mass of ≈17.5 M_⊙ , which is most likely an upper limit. Estimating the inclination from tidal synchronization yields a BH mass of ≈13.8 M_⊙ , which is also close to the peak mass ≈14.2 M_⊙ predicted by our Monte Carlo simulation probing a range of inclinations; We found evidence for N and He enrichment at the surface of the donor. This is likely a result of the combined effects of material transfer from the BH progenitor and subsequent efficient mixing driven by rapid rotation;

• The analysis reveals that Cyg X-1 exhibits super-solar abundances of Fe, Si, and Mg, approximately 1.3–1.8 times higher than solar values. Notably, while there is evidence of CNO-processed material on the surface (nitrogen enrichment, oxygen depletion), the total CNO abundance remains consistent with solar values. This distinct abundance pattern sets Cyg X-1 apart from stars in its surrounding Cyg OB3 association;

• While a supernova from the BH progenitor could potentially enrich the donor star, the observed abundance pattern (specifically the simultaneous enrichment of Fe, Si, and Mg without a corresponding increase in total CNO) is difficult to reconcile with current models of supernova yields and pre-supernova wind yields. This implies that the system had a high initial metallicity and likely formed from an enriched ISM;

• Our analysis yields a terminal wind velocity of ≈1200 km s⁻¹ and a mass-loss rate of $\dot{M}\approx 3\times {10}^{-7}{M}_{\odot }{\text{yr}}^{-1}$ in the high-soft state. We found evidence for faster and stronger wind in low-hard state (ν_∞ ≲ 1800 km s⁻¹ and $\dot{M}\approx 5\times {10}^{-7}{M}_{\odot }{\text{yr}}^{-1}$). These estimates are significantly lower compared to the commonly adopted values in the literature for this system. This reduction is primarily attributed to several key improvements in our modeling approach: (1) the inclusion of X-ray photoionization effects, (2) the use of advanced, non-LTE atmospheric models with stellar winds (3) the simultaneous analysis of UV and optical wind features, (4) consideration of wind clumping, and (5) a detailed exploration of the wind velocity exponent. This derived mass-loss rate is lower than Vink et al. (2000) predictions by a factor of ≈ 2–3, further emphasizing the impact of clumping;

• The X-ray luminosity of Cyg X-1 is consistent with being driven primarily by wind accretion. Our Bondi-Hoyle calculations predict a higher accretion luminosity than observed, suggesting that not all accreted material contributes to the X-ray emission, and the persistent presence of wind features in the UV spectra, even under strong X-ray ionization, confirms the dominance of wind-driven mass transfer over RLOF. Furthermore, the UV observations are inconsistent with the presence of a focused wind toward the BH;

• We find that the radius of the donor star is larger than the equivalent spherical radius of its Roche lobe, adopting a classical description. However, the strong radiation force of the O-star may significantly change its effective potential and thereby modify the typical Roche lobe geometry. The donor star radius is smaller than the radius of the inner Lagrangian point for all plausible values of inclination and the BH mass. We thus conclude that there is no evidence for a significant focused direct overflow towards the black hole in Cyg X-1;

• The observed variations in wind velocity and mass-loss rate, coupled with corresponding changes in X-ray luminosity between high-soft and low-hard states, point to a complex interplay between the stellar wind, accretion, and X-ray states in Cyg X-1;

• MESA evolutionary models constrained by newly derived stellar parameters show that Cyg X-1 is currently in a pre-RLOF phase, consistent with the observed X-ray luminosity. However, these models predict that the system will undergo RLOF in the near future, leading to a significant increase in X-ray luminosity, potentially becoming a ULX. The subsequent evolution depends strongly on the wind mass-loss rate of the donor and mass accretion efficiency onto the BH;

• We consider Cyg X-1 as a potential gravitational wave source progenitor. While Eddington-limited mass accretion scenarios lead to limited BH growth and wider orbits (merging time > 51 Gyr), the models assuming fully conservative mass accretion predict a dramatic increase in the BH mass, a modest increase in the orbital period, and a final binary BH system that will merge within approximately ≳5 Gyr, making Cyg X-1 a potential gravitational wave source. This highlights the critical role of mass accretion efficiency in determining the system’s potential as a gravitational wave source.

Taken together, these results paint a complex and interconnected picture of the Cyg X-1 system, highlighting the crucial role of stellar winds, accretion processes, and initial conditions in shaping its evolution and ultimate fate.

Acknowledgements

VR, AACS, and MBP are supported by the Deutsche Forschungsgemeinschaft (DFG – German Research Foundation) in the form of an Emmy Noether Research Group – Project-ID 445674056 (SA4064/1-1, PI Sander). ECS acknowledges financial support by the Federal Ministry for Economic Affairs and Climate Action (BMWK) via the German Aerospace Center (Deutsches Zentrum für Luftund Raumfahrt, DLR) grant 50 OR 2306 (PI: Ramachandran/Sander). DP acknowledges financial support from the Research Foundation Flanders (FWO) under grant agreement No. 1256225N. B.K. gratefully acknowledges support from the Grant Agency of the Czech Republic (GAČ R 25-15910S). The Astronomical Institute of the Czech Academy of Sciences in Ondřejov is supported by the project RVO:67985815.

Appendix A Additional tables

Appendix B Additional plots

Fig. B.1 Impact of ξ on metal line profiles.

Fig. B.2 Optical spectra (blue) taken at orbital phase ϕ ∼ 0 and high-soft state compared to best-fit model spectra (red).

Fig. B.3 Optical spectra (blue) taken at orbital phase ϕ ∼ 0.5 and high-soft state compared to best-fit model spectra (red).

Fig. B.4 Optical spectra (blue) taken at orbital phase ϕ ∼ 0 and low-hard state compared to best-fit model spectra (red).

Fig. B.5 Optical spectra (blue) taken at orbital phase ϕ ∼ 0.5 and low-hard state compared to best-fit model spectra (red).

Fig. B.6 HST UV spectrum taken at ϕ ≈ 0 (upper) and ϕ ≈ 0.5 (lower) compared to best-fit model spectra.

Fig. B.7 Probability distributions of Roche-lobe filling factors ( fRL1, fRbk, fRy, fRz). Roche-lobe radii in various directions were calculated using the Leahy & Leahy (2015) formalism, with mass ratio (q) distributions derived from Monte Carlo simulations incorporating the BH mass probability distribution (Fig. 8). Shaded regions represent 1-sigma confidence intervals, and dashed green lines indicate distribution peaks.

References

All Tables

All Figures

	Fig. 1 Hα line variability observed at different orbital phases and X-ray spectral state of Cyg X-1 as described in the legend.
In the text

	Fig. 2 Comparison of observed optical spectra (blue) to synthetic spectra obtained for models with different stellar temperatures and surface gravities (see legend). The plot illustrates the effect of temperature changes on Si III to Si IV and He I to He II line ratios.
In the text

	Fig. 3 Same as Fig. 2, but illustrating the impact of log g on wings of Balmer lines.
In the text

	Fig. 4 Spectral energy distribution for Cyg X-1. The model SED (red solid line) is adjusted to fit the flux-calibrated HST spectra in the FUV (blue line), IUE spectra in the NUV (black line), and photometry from UV infrared bands (blue boxes).
In the text

	Fig. 5 Iron forest in the FUV region. The blue line represents the observed HST spectra, while the best-fit model, obtained using Bayesian analysis, is shown in red for comparison.
In the text

	Fig. 6 Elemental abundances in Cyg X-1 donor O star relative to hydrogen (upper panel) and iron (lower panel). The abundance ratios are relative to solar and in logarithmic scale. For comparison, the average chemical abundance of B stars in the Cyg OB3 (Daflon et al. 2001) and B stars in the solar neighborhood (Nieva & Przybilla 2012) are marked.
In the text

Fig. 7 Mass plane diagram for Cyg X-1, showing component masses and their uncertainties. BH masses are estimated for i = 27.5°(blue) and i = 33.7° (orange).The solid black line represents the mass relation for the adopted inclination of i = 27.5° derived from optical observations (Orosz et al. 2011; Miller-Jones et al. 2021). The dotted line corresponds to i = 33.7° derived assuming synchronous rotation. The dashed lines illustrate the variation in the mass relation for alternative inclination angles.

In the text

	Fig. 8 Probability distribution of BH mass in Cyg X-1. The gray histogram and blue line depict the MBH distribution and its KDE, respectively. The shaded region shows the 1-σ range (16th and 84th percentiles). The dotted green line marks the KDE peak, and the dashed orange and blue lines represent the MBH values at inclinations of 27.5° and 33.7° (see legend).
In the text

	Fig. 9 Illustration of the Roche lobe geometry for Cyg X-1, ing the O-star donor (light blue circle) and the BH (black do MBH = 17.5 M⊙). The dark blue lines and the mesh depict the ap mate Roche equipotential surface, with arrows indicating the dist RL1, Rbk, and Ry. The axes are scaled in R⊙.
In the text

Fig. 10 Comparison of UV P-Cygni profiles observed with HST (blue) at orbital phase ϕ∼ 0 (top panel) and ϕ ∼0.5 (bottom panel) during high-soft state to the best-fit models (red). The models assume the same input parameters except for the incorporated amount of X-rays. The wind parameters used in the best-fit models are log $\log \dot{M}(M_{\odot }{\text{yr}}^{-1})$ and ν∞ =1200 km s−1. For superior conjunction ϕ ∼ 0 we adopted LX ∼ 7 × 1036 erg s−1 in the model, and at inferior conjunction ϕ∼ 0.5 we used higher amount of X-rays LX∼ 2 ×1037 erg s−1. The variability in the wind line profiles shows the impact of X-ray ionization.

In the text

	Fig. 11 Same as Fig. 10 but IUE observations at orbital phase ϕ ∼ 0 during low-hard state are shown. We adopted a lower X-ray luminosity LX ∼ 3--5 × 1036 erg s−1 in the model. Adopted wind parameters in models are shown in legend.
In the text

	Fig. 12 Potential future evolution of the Cyg X-1 donor star. Here, the evolution starts when the primary is already a BH (point mass) of 17.4 M⊙, O star donor of initial mass 34 M⊙, and an initial period of 5.5 days. Three models are displayed here with mass-loss scaling factors of 1, 0.5, and 0.2 to demonstrate the effect of wind mass loss on evolution.
In the text

	Fig. 13 Left panel: same as Fig. 12 but here the tracks correspond to a wind mass-loss scaling factor of 0.5, and the color scaling shows the accretion luminosity. Right panel: accretion luminosity (top) and Roche-lobe filling factor (bottom) as a function of normalized donor lifetime τ. The gray bar shows the range of observed X-ray luminosity of the system.
In the text

	Fig. 14 Impact of mass accretion by the BH on the evolution of Cyg X-1.The panels show the evolution of accretion luminosity, orbital period, and component masses as a function of normalized time τ for MESA models with Eddington-limited (blue) and fully conservative (red) mass accretion.
In the text

	Fig. B.1 Impact of ξ on metal line profiles.
In the text

	Fig. B.2 Optical spectra (blue) taken at orbital phase ϕ ∼ 0 and high-soft state compared to best-fit model spectra (red).
In the text

	Fig. B.3 Optical spectra (blue) taken at orbital phase ϕ ∼ 0.5 and high-soft state compared to best-fit model spectra (red).
In the text

	Fig. B.4 Optical spectra (blue) taken at orbital phase ϕ ∼ 0 and low-hard state compared to best-fit model spectra (red).
In the text

	Fig. B.5 Optical spectra (blue) taken at orbital phase ϕ ∼ 0.5 and low-hard state compared to best-fit model spectra (red).
In the text

	Fig. B.6 HST UV spectrum taken at ϕ ≈ 0 (upper) and ϕ ≈ 0.5 (lower) compared to best-fit model spectra.
In the text

Fig. B.7 Probability distributions of Roche-lobe filling factors ( fRL1, fRbk, fRy, fRz). Roche-lobe radii in various directions were calculated using the Leahy & Leahy (2015) formalism, with mass ratio (q) distributions derived from Monte Carlo simulations incorporating the BH mass probability distribution (Fig. 8). Shaded regions represent 1-sigma confidence intervals, and dashed green lines indicate distribution peaks.

In the text