© The Authors 2024

1. Introduction

High-mass X-ray binaries (HMXBs) are a class of binary systems in which a neutron star or a black hole (BH) is fed by the stellar wind from a massive (≳10 M_⊙) OB companion or by the decretion disc surrounding a quickly spinning Be donor star. The stellar wind is line-driven (Castor et al. 1975) and can be accelerated up to terminal velocities of v_∞ ∼ 2500 km s⁻¹ (e.g. Martínez-Núñez et al. 2017). Internal shocks lead to the formation of overdense regions, or clumps, making the wind inhomogeneous and highly structured (e.g. Owocki et al. 1988; Feldmeier et al. 1997; Oskinova et al. 2012; Sundqvist & Owocki 2013). While several stellar wind models have been developed (e.g. Oskinova et al. 2012; Sundqvist & Owocki 2013; Sundqvist et al. 2018; El Mellah et al. 2018, 2020), the physical properties of the clumps, such as their shape, density, and ionisation structure, are yet to be precisely constrained. Investigating the structure of the stellar winds produced by massive stars is crucial for various purposes, including constraining models of stellar and galactic evolution, and understanding the accretion mechanism in wind-fed HMXBs.

As a matter of fact, stellar winds represent an important mechanism of stellar mass loss, with mass-loss rates up to 10⁻⁵ M_⊙ yr⁻¹ (Puls et al. 2008; Martínez-Núñez et al. 2017). This has an influence on the properties of the donor star (such as mass, luminosity, spin, and chemical composition) as well as its evolutionary timescales. Moreover, in binary systems interacting via mass transfer mediated by a stellar wind like most HMXBs, only a fraction of the expelled gas is accreted, while the rest escapes the system, causing a removal of angular momentum and a change in the orbital parameters. Whether this process leads to a shrinking or a widening of the orbit (e.g. Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; Paczynski 1976; El Mellah et al. 2020) depends on key parameters, such as the wind terminal velocity, the initial orbital separation, and the accretion efficiency (e.g. Saladino et al. 2018). The evolutionary pathways of such interactions can have significant impact on the rate of merging events and guide us in the quest for merger progenitors (e.g. Bulik et al. 2011; Belczynski et al. 2011; Neijssel et al. 2021).

The well-known HMXB system Cygnus X-1 (Cyg X-1 hereafter) is a good target to constrain the stellar wind structure and infer information on the physical properties of wind clumps. The system hosts a BH with M_BH = 21.2 ± 2.2 M_⊙ (Miller-Jones et al. 2021) in a ∼5.6 d quasi-circular orbit (Gies et al. 2003). The O9.7 Iab supergiant companion star HDE 226868, with $M_{*}=40.6_{-7.1}^{+7.7}{\mathrm{M}}_{\odot }$ and R_* = 22.3 ± 1.8 R_⊙ (Miller-Jones et al. 2021), launches fast stellar winds (v_∞ = 2100 km s⁻¹, Herrero et al. 1995). The orbital inclination of the system (i ∼ 27°, Miller-Jones et al. 2021) causes our line of sight (LOS) to intercept the clump-forming region, likely located close to the companion star’s surface (e.g. El Mellah et al. 2020).

In this system, the presence of the wind significantly modifies the spectral (e.g. Nowak et al. 2011) and timing (Lai et al. 2022) properties of the X-ray source. Indeed, X-ray photons undergo different levels of absorption depending on the amount and the ionisation state of the wind material intercepting the LOS. Wind absorption mostly affects the softest energy bands (E ≲ 1 keV), although its influence can be observed up to higher energies of E ∼ 8 − 10 keV (e.g. Hirsch et al. 2019; Grinberg et al. 2015). The strength of the absorption events strongly varies as a function of the orbital phase (e.g. Bałucińska-Church et al. 2000; Poutanen et al. 2008; Hirsch et al. 2019; Grinberg et al. 2020; Lai et al. 2022). In particular, the X-ray light curves show more intense and frequent absorption dips when the X-ray source is at superior conjunction (ϕ_orb = 0), that is when the BH is the farthest and the companion star the closest from the observer along the orbit (see Li & Clark 1974 and references therein). In Lai et al. (2022), we showed that the absorption dips significantly contribute to the X-ray variability of Cyg X-1 by increasing its observed fractional variability on timescales longer than ∼1 s, and suggested that the motion of clumps of the size of ∼10⁻⁴ R_* may be responsible for the enhancement of X-ray variability on these timescales.

Radiative hydrodynamical simulations show that the size of the clumps should increase as they move away from the star (Sundqvist et al. 2018), producing dips of increasingly long duration. Therefore, to fully investigate the properties of the clumps, we need to access a broad range of timescales, down to the shortest timescales (∼1 min or shorter) characterising the smallest clumps formed close to the base of the wind. However, the limitations of current instruments preclude wind-induced spectral changes occurring on very short timescales from being investigated via canonical spectral analysis.

An alternative, powerful approach is to use colour-colour diagrams (e.g. Nowak et al. 2011; Grinberg et al. 2020). Detailed modelling of the colour-colour diagram of wind-fed systems can give constraints on the physical properties of the wind. Earlier studies made use of empirical functions or simple (neutral) stellar wind models to describe the characteristic shape of colour-colour tracks (e.g. Hanke et al. 2008; Nowak et al. 2011; Hirsch et al. 2019). Recently, Grinberg et al. (2020) explored more complex wind models, concluding that a proper description of the observed colour-colour tracks requires accounting for different levels of ionisation in the wind. In this paper, building on the work of Grinberg et al. (2020), we model the colour-colour diagrams of Cyg X-1 in order to constrain the parameters of the wind as a function of the orbital phase. This information can be used to infer physical properties such as the stellar mass-loss rate and the mass of the clumps (e.g. El Mellah et al. 2020).

In this paper, we first describe the data reduction procedure (Sect. 2), then we test a number of stellar wind models (Sect. 3) to determine the one that best describes the colour-colour diagram of Cyg X-1 in its hard state (the most affected by the wind, Nowak et al. 2011). In Sect. 4, we investigate the temporal evolution of the column density and covering factor of the wind by fitting time-resolved colour-colour diagrams. We discuss our results in Sect. 5 and present our conclusions in Sect. 6.

2. Observation and data reduction

We focussed our analysis on a single XMM-Newton observation of Cyg X-1 (obsID 0745250201, hereafter 201) in its hard state, which is part of the ∼1.5 orbital period-long monitoring from the multi wavelength campaign “Cyg X-1 Hard state Observations of a Complete Binary Orbit in X-rays” (CHOCBOX). The full CHOCBOX coverage is shown in Fig. 1. The monitoring consists of four pointings for a total duration of ∼7 days, and catches two consecutive passages of the X-ray source at superior conjunction.

Fig. 1. Orbital phase coverage of the XMM-Newton CHOCBOX monitoring of Cyg X-1 in the hard state. The four arc labels correspond to the ObsIDs of each pointing (see Table 1 in Lai et al. 2022 for the full ObsID). The orbital phase, ϕorb = 0, indicates the passage at superior conjunction. Observation 201 (in purple) is the focus of the analyses presented in this paper. The starting orbital phase of the monitoring (corresponding to ϕorb = 0.82) is also indicated. HDE 226868 and Cyg X-1 are only schematically represented, no Roche lobe was graphically considered in the scheme.

In Lai et al. (2022), we investigated the changes driven by the stellar wind in the X-ray spectral-timing properties of the source throughout the entire monitoring. Among all the CHOCBOX observations, observation 201 turned out to be the most affected by wind absorption. It also showed the most prominent effects on the X-ray variability on short timescales (in the ∼0.1 − 10 s range) as compared to the other observations. Observation 201 covers orbital phases ϕ_orb = 0.82 − 0.06, corresponding to the first passage at superior conjunction (Lai et al. 2022). We focussed on observation 201 for the study of the stellar wind, although in Sect. 4 we also used observation ID 0745250501 (hereafter 501) to identify the region of the colour-colour diagram less affected by wind-absorption. Indeed, observation 501 covers the orbital phases ϕ_orb = 0.17 − 0.46. Being the closest to inferior conjunction with ϕ_orb = 0.5 it is also the least absorbed. For both observations we used EPIC-pn data (Strüder et al. 2001) in timing observing mode. We performed the data reduction using the XMM-NewtonScience Analysis Software (SAS, version 20.0.0) and calibration files (CCF) as of 2022 March.

Diez et al. (2023) showed that the extraction of EPIC-pn spectra in timing mode using the default Rate Dependent PHA (RDPHA) correction led to energies higher than expected for line features, and suggest a non-standard calibration to mitigate the problem. We followed the same approach proposed by Diez et al. (2023). We ran epproc on the Observation Data Files (ODFs), turning off the RDPHA correction (withrdpha = ‘N’) and applying the Rate Dependent CTI (RDCTI) correction using epfast (runepfast = ‘Y’). Calibrated event files were screened for the presence of strong background flaring events using light curves with 1s time resolution in the 10 − 15 keV energy band. No significant flares were observed. Using the SAS task epatplot, we found that a small fraction of pile-up is present in the data. We corrected for it by extracting source counts from a region that excludes the central pixels (RAWX = [36:39]); that is RAWX = [30:35; 40:46]. Background counts cannot be reliably estimated from timing mode observations of bright sources, as the entire CCD is illuminated by source photons. Since the PSF is energy-dependent, extracting the background from the outer columns is not recommended as this procedure modifies the intrinsic spectrum (Ng et al. 2010). Since the source is very bright, with ≳220 counts s⁻¹ in the 0.5 − 10 keV energy band, the background contribution is expected to be negligible. Therefore, we decided to not account for the background in our analysis. We made use of the SAS tasks arfgen and rmfgen to extract the Ancillary Response Files (ARF) and Redistribution Matrix Files (RMF). We produced two ARF files: one for the full region (RAWX = [30:46]) and one for the central excluded region (RAWX = [36:39]). Using the command addarf, we subtracted the latter from the former, generating the final ARF file (e.g. Wilkinson & Uttley 2009; Lai et al. 2022). The spectrum was rebinned to ensure a minimum of 20 counts per bin in order to apply chi-square statistics.

Spectral fits were done using XSPEC v12.10.1 (Arnaud 1996). We used the Interactive Spectral Interpretation System (ISIS) version 1.6.2 and a custom code, written in Python 3.10, for the calculation and the modelling of the colour-colour diagrams. Time series were extracted using stingray v1.1.2 (Huppenkothen et al. 2019a,b; Bachetti et al. 2023).

3. Models for the colour-colour diagram of Cyg X-1

The passage of absorbing clumps across the LOS manifests as dips in the X-ray light curve of length ranging between several seconds to hours (e.g. Hirsch et al. 2019). To explore the presence and duration of dips in observation 201, we show its light curve in Fig. 2. In the upper and middle panels, we report the soft (0.5 − 1.5 keV) and the hard (3 − 10 keV) band light curves, respectively (for the light curves of the entire modelling, see Lai et al. 2022¹).

Fig. 2. XMM–Newton EPIC-pn light curves of the observation 201 of Cyg X-1 with time bins of 10 s. The upper and the middle panels show, respectively, the 0.5 − 1.5 keV and the 3 − 10 keV light curves. The bottom panel reports the hardness, i.e. the ratio between count rates in the 3 − 10 keV and 0.5 − 1.5 keV energy bands. The vertical red line indicates the passage at superior conjunction (i.e. ϕorb = 0).

The presence of dips in the soft band light curve, which are approximately three times less intense in the hard band light curve, produces sharp increases in the hardness ratios (as is shown in the bottom panel of Fig. 2, where we choose to compare two non-contiguous energy bands to maximise the intensity of the dips). The single dips are quite short (≲1 ks). Therefore, given the swiftness of such events, their energy spectrum is difficult to study via time-resolved spectroscopy. An alternative approach is to investigate the broad band spectral variability caused by wind clumps using colour-colour diagrams (e.g. Nowak et al. 2011). These diagrams map the evolution of the ratio of the count rates in different energy bands. Specifically, a soft colour, defined as the ratio between count rates in a soft and an intermediate energy band, is plotted against a hard colour, defined as the ratio between count rates in the same intermediate band and a hard energy band. High values of both soft and hard colours characterise the least absorbed phases, and thus map the upper region of the diagram. As absorption along the LOS increases, the hard and soft colours change. In particular, for a partially covering absorber the resulting colour-colour tracks describe a ‘pointy’ or ‘nose-like’ shape as the source becomes more absorbed (i.e. in the lower region of the colour-colour diagram, see e.g. Hirsch et al. 2019). Following previous studies (Hanke et al. 2008; Nowak et al. 2011; Hirsch et al. 2019; Grinberg et al. 2020), we constructed the colour-colour diagram of observation 201 using light curve bins of 10 s, allowing us to detect also the spectral variability that could be caused by small clumps. The size of the clumps is directly related to the fly-by time across the LOS (El Mellah et al. 2020). With this time resolution, the minimum size of the clumps that can be tested is ∼2 × 10⁻⁴ R_*, assuming a wind terminal velocity of v_∞ = 2400 km s⁻¹ (Grinberg et al. 2015), where R_* is the radius of the companion star. The ratios were computed between the energy bands 0.5 − 1.5/1.5 − 3 keV (soft colour), and 1.5 − 3/3 − 10 keV (hard colour), see Fig. 3.

Fig. 3. Simulated tracks for a power law (with Γ fixed at 1.6, left panel, 1.7, middle panel, and 1.8, right panel) plus a neutral absorber, compared to the observed (grey points) colour-colour diagram tracks of observation 201 of Cyg X-1. The blue (red) tracks describe changes in the wind fc (NH, w) for a fixed value of NH, w (fc).

We then tested a number of physical scenarios in order to find a suitable stellar wind model able to describe the observed colour-colour diagram tracks in Cyg X-1. As our baseline model we chose a primary X-ray continuum modified by a partially covering absorber (Fürst et al. 2014; Fornasini et al. 2017; Grinberg et al. 2020). The adopted model is of the form:

$$\begin{array}{c}\hfill {\mathtt{abs}}_{\mathrm{ism}}\times \mathtt{continuum}\times (f_{\mathrm{c}}\times {\mathtt{abs}}_{\mathrm{wind}}+(1-f_{\mathrm{c}}))\end{array}$$(1)

where abs_ism is the absorption of the interstellar medium (ISM), continuum refers to the unabsorbed primary emission from the X-ray source and abs_wind is the absorption component associated with the stellar wind. This component blocks a fraction (f_c, or “covering fraction”) of the X-ray source photons, while only a percentage 1 − f_c reaches the observer modified only by the interstellar medium.

We made use of standard XSPEC models (Arnaud 1996) for the continuum and for the stellar wind to simulate spectra within ISIS (Houck & Denicola 2000; Houck 2002; Noble & Nowak 2008), using the response matrices of observation 201 (see Sect. 2). We let the parameters of interest of the model vary and, for each combination of parameters, we calculate the corresponding soft and hard colours. In order to speed up the process, for each set of simulations, only one parameter is left free to vary, while the others are kept fixed. This procedure allows us to build simulated colour-colour tracks, which are then compared to the data. We tested different models as is described in the following sections. Throughout this work ISM absorption was modelled with TBabs² using the wilm abundances (Wilms et al. 2000) and the vern cross-sections (Verner et al. 1996). The column density of the ISM was kept fixed at the tabulated value of ∼0.7 × 10²² cm⁻² (HI4PI Collaboration 2016). We note that in this paper we are interested in testing for variability of the stellar wind component; therefore, in our simulations, we do not investigate intrinsic spectral variations of the X-ray source.

3.1. Model 1: Single power law plus neutral stellar wind

We first simulated theoretical colour-colour tracks assuming a power law for the underlying continuum and a partially covering neutral absorber. For this simple model, the shape of the colour-colour tracks mainly depends on three parameters: the power law photon index Γ of the continuum, the wind covering fraction f_c and the wind column density N_H, w.

Fig. 3 shows the resulting simulated tracks overplotted on the colour-colour diagram of observation 201. For fixed values of Γ and N_H, w, a variable f_c describes the curvature of the blue-cyan tracks shown in Fig. 3. Within each track the f_c parameter samples the range of values f_c = 0.1 − 1 (at steps of 0.1). The different tracks from deep blue to light cyan are obtained by increasing the value of N_H, w from a minimum value of N_H, w = 0.1 × 10²² to a maximum value of N_H, w = 32 × 10²² cm⁻² (for higher values of N_H, w hard colours start to increase again, producing an upward tail that is not observed in the data). On the other hand, for fixed values of Γ and f_c, a variable N_H, w, within the range N_H, w = 0.1 − 32 × 10²²cm⁻², describes the curvature of the red-orange tracks shown in Fig. 3. The different tracks (from deep red to light orange) are obtained by increasing the value of f_c (from 0.1 to 1). We repeated the same simulations for Γ = 1.6, 1.7, and 1.8 (panels from left to right in Fig. 3), which match typical values for the hard state (e.g. Joinet et al. 2008; Motta et al. 2009; Gilfanov 2010; Basak et al. 2017; Zhou et al. 2022).

We confirm the results of previous studies (e.g. Hirsch et al. 2019; Grinberg et al. 2020), which showed that this simple model fails to describe the observed colour-colour diagram. In particular, a variable N_H, w, while approximating the curvature of the observed track, does not properly reproduce its pointy shape at low hard colours (giving the track a nose-like shape). We observe that a varying Γ shifts the tracks vertically. Therefore, in line with previous studies (e.g. Grinberg et al. 2020), we infer that the data require a more complex modelling of the stellar wind and/or the intrinsic X-ray continuum. Indeed, several X-ray spectroscopic studies highlight the need for a more complex primary continuum for Cyg X-1 (e.g. Basak et al. 2017; Tomsick et al. 2018) as well as a structured stellar wind made of differently ionised material (Hirsch et al. 2019).

3.2. Model 2: Structured continuum plus ionised stellar wind

We then changed our assumptions on the underlying primary continuum and tested a more complex model. To this aim, we assumed the model presented in Lai et al. (2022), which includes a soft and a hard Comptonisation component, as well as direct thermal and reflected emission from the disc. The employed XSPEC model is: TBabs × [diskbb + nthComp + relxillCp]. In order to properly assess the time-averaged parameters of the continuum we fit a subset of good time intervals (GTIs), selected so as to be the least affected by the intervening stellar wind. To this aim, we adopted the same criterion defined in Lai et al. (2022); that is, we selected data with hard colour ≥0.95 and soft colour ≥0.7, so as to single out GTIs from the upper, least absorbed region of the colour-colour diagram. Details about the continuum model used in these simulations are reported in Appendix A.

Fig. 4 shows the colour-colour tracks computed for a neutral gas partially absorbing the adopted structured continuum. As opposed to simulations shown in Fig. 3, the new tracks all start right at the centre of the observed distribution of data points (upper part of the diagram), meaning that the new continuum model describes the unabsorbed time-averaged spectrum much better than a simple power law. Therefore, hereafter we use this more complex model to describe the continuum in our simulations.

Fig. 4. Simulated tracks for the complex continuum absorbed by a neutral wind model, as compared to the colour-colour diagram of Cyg X-1.

We then considered the effects of an ionised wind on the simulated colour-colour tracks, substituting the abs_wind component in Eq. (1) with the warmabs (v2.31) analytic photoionisation model (Kallman et al. 2009). As in Grinberg et al. (2020), we used the standard population files delivered with warmabs, with densities of 10¹² cm⁻³, typical of stellar wind densities close to the compact object (Lomaeva et al. 2020). We first verified how different levels of ionisation modify the simulated tracks. To this aim, we considered the ionisation parameter ξ = L_X/nr² (Tarter et al. 1969), where L_X is the ionising luminosity above 13.6 eV, n is the absorbing gas density, and r its distance from the ionising X-ray source. Values of log ξ between −4 and 2 were considered. For each fixed value of log ξ, we let the N_H, w vary between 0.1 and 32 × 10²² cm⁻², leading to the tracks shown in Fig. 5. We reproduced the same simulations for three different values of covering factor (f_c = 0.8, 0.87, 0.95, Fig. 5 from left to right).

Fig. 5. Simulated tracks for a homogeneously ionised absorber (with log ξ values ranging between −4 and 2) compared to the colour-colour diagram of observation 201. The different panels show the simulated tracks computed for different values of fc: 0.8 (left), 0.87 (middle), 0.95 (right). The black curve in the middle panel represents the one that best describes the data.

By increasing the ionisation parameter, the simulated tracks become more pointy, thus more closely resemble the observed tracks. The effect of increasing the covering factor is mostly to increase the range of encompassed hard and soft colours. In particular, for sufficiently high N_H, w, a higher f_c produces more absorption at low-to-intermediate energies, causing the absorption dips to be more spectrally hard (lower values of hard and soft colours during the most absorbed stages). Among the simulated curves, we observe that the one better resembling the observed tracks corresponds to f_c = 0.87, log ξ = 1.6, and N_H, w varying between 0.1 and 32 × 10²² cm⁻², as is highlighted in black in the middle panel of Fig. 5, suggesting mild ionisation for the partially covering wind.

3.3. Model 3: Stellar wind with variable ionisation

The simulations shown in Sect. 3.2 assumed a constant ionisation parameter throughout the different stages of the absorption dips. Nonetheless, Hirsch et al. (2019) reported the appearance of less ionised species in the most absorbed stages corresponding to the lower part of the colour-colour diagram. This suggests the presence of structured clumps, made of differently ionised material. Such results highlight the need to consider changes in the ionisation parameter (Grinberg et al. 2020). Therefore, we now account for variations in the stellar wind ionisation parameter as a function of N_H, w.

Following Grinberg et al. (2020), we test two empirical functions to describe the dependence of the ionisation parameter on the column density of the absorbing material. The first function is defined as:

$$\begin{array}{c}\hfill log\xi =log(A/[N_{\mathrm{H},\mathrm{w}}/{10}^{22}{\mathrm{cm}}^{-2}]),\end{array}$$(2)

with A > 0. It is based on the definition of the ionisation parameter and assumes a linear decrease in ξ as a function of N_H, w. The second function is defined as:

$$\begin{array}{c}\hfill log\xi =log\frac{B+[N_{\mathrm{H},\mathrm{w}}/{10}^{22}{\mathrm{cm}}^{-2}]}{[N_{\mathrm{H},\mathrm{w}}/{10}^{22}{\mathrm{cm}}^{-2}]}+C,\end{array}$$(3)

with B, C > 0. It assumes that the ionisation parameter deviates from a linearly decreasing trend at high column densities, staying relatively high, as expected if an additional source of ionisation plays an important role in denser environments (e.g. Feldmeier et al. 1997). We arbitrarily chose the values of the constant parameters to be A = 100, B = 10, and C = 1. This choice of parameters is such that the two functions match at low column densities N_H, w/10²² cm⁻² < < 10, and start deviating significantly at higher densities, N_H, w/10²² cm⁻² ≳ 1. We let N_H, w vary between 0.01 and 32 × 10²² cm⁻². Sampling down to values of N_H, w = 0.01 × 10²² cm⁻² allows us to test values of ionisation parameter as high as the highest value allowed by the warmabs model. Therefore, the first function spans the range of log ξ ∼ 0.5 − 4, while the second function spans the range of log ξ ∼ 1.12 − 4. The two empirical functions are shown in Fig. 6.

Fig. 6. Empirical functions describing the dependence of the ionisation parameter on the column density of the absorbing gas. Left panel: Assumed dependencies of the ionisation parameter, log ξ, on the column density, NH, w. Middle and right panel: Simulated tracks for a warm absorbing gas, assuming log ξ = log(100/[NH, w/1022 cm−2]) (middle panel), and $log\xi =log\frac{10+[N_{\mathrm{H},\mathrm{w}}/{10}^{22}{\mathrm{cm}}^{-2}]}{[N_{\mathrm{H},\mathrm{w}}/{10}^{22}{\mathrm{cm}}^{-2}]}+1$ (right panel). The different shades correspond to different values of covering factor as reported in the labels.

We built models comprising the complex continuum defined in Sect. 3.2 plus each of the two empirical functions describing variable ionisation. Since the data in the colour-colour diagrams are not normally distributed, a simple χ² minimisation method cannot be used to select the best-fit model. In order to find the best-fit model without any prior assumption on the distribution of the data we employed a non-parametric method.

We used the kernel density estimation (KDE) method to infer the probability density function (PDF) of the data in the colour-colour diagram (e.g. Hill 1985; Cavecchi & Patruno 2022). The KDE is a non-parametric estimator, and its functional form is obtained by combining as many building blocks – kernels – as the number of data points. In the KDE a single type of kernel K is used. The choice of the kernel is arbitrary, but for large datasets this choice does not have any significant effect on the final output. Therefore, we used a simple kernel with a Gaussian shape. The kernel K_i centred on the data point x_i is defined as

$$\begin{array}{c}\hfill K_i(\mathit{x})\propto exp(-\frac{{(\mathit{x}-{\mathit{x}}_i)}^2}{2h^2})\end{array}$$(4)

where x is the point in the colour-colour diagram at which the kernel function is estimated. The width of the kernel, h, controls the smoothing. The normalised PDF is then computed as:

$$\begin{array}{c}\hfill \mathrm{PDF}(\mathit{x})=\frac{1}{n}{\displaystyle \sum _i^n}{K}_i(\mathit{x}),\end{array}$$(5)

which sums the contribution of the kernels centred on each data point, with n being the total number of data points. We implemented the method using the SciPy function gaussian_kde, with the optimal kernel width automatically determined by the function (default parameter bw_method = ‘scott’³). Using the KDE method we estimated the PDF of the data. Then, for each simulated model we computed the combined likelihood as the product of the values of the probability density at each point of the model. The best-fit model is the one characterised by the highest combined likelihood.

In Fig. 6, we show the best-fit track obtained using the functions log ξ = log(100/[N_H, w/10²² cm⁻²]) (middle panel) and $log\xi =log\frac{10+[N_{\mathrm{H},\mathrm{w}}/{10}^{22}{\mathrm{cm}}^{-2}]}{[N_{\mathrm{H},\mathrm{w}}/{10}^{22}{\mathrm{cm}}^{-2}]}+1$ (right panel) to model the ionised absorber. The best-fit track is obtained in both cases by fitting for the covering factor f_c (between 0.7 and 0.95, at steps of 0.01) parameter while N_H, w spans the entire range 0.01 − 32 × 10²² cm⁻² for each fit value of the covering factor. In Fig. 6, we also show the tracks corresponding to f_c = 0.7 and f_c = 0.95 and note that outside these limiting values the simulated curves do not intersect the data in the lower part of the diagram.

During the least absorbed stages (i.e. at hard colours ≳0.8), the simulated tracks are not very sensitive to the covering factor, and they all resemble fairly well the observed shape of the track in this part of the diagram. However, at lower values of hard colours, corresponding to the most absorbed phases, the dependence on the covering factor becomes significant. The two variable ionisation models considered here describe the data well, with a mild 5% difference between their inferred probabilities, which does not allow us to prefer either one of the two models. Notably, both functions resemble the pointy shape of the data much better than the homogeneously ionised absorber considered in Sect. 3.2, giving more support to the hypothesis of a structured absorber. Finally, we observe that the data all lay within the simulated tracks for f_c = 0.7 and f_c = 0.95, implying that the scatter of data points in this part of the diagram might be in part due to intrinsic variations in the covering factor between ∼0.7 − 0.95. We investigate this problem further in the following section.

4. Time-resolved colour-colour diagrams

To study the phase-dependent variability of the absorbing column density and covering fraction of the wind, we investigated the temporal evolution of the shape of the colour-colour diagram. We divided observation 201 in 10 segments of the same duration of ∼11 ks each, corresponding to segments of ∼0.024 in orbital phase, and produced colour-colour diagrams for each segment. The corresponding tracks are all shown in Fig. 7, colour-coded based on the range of orbital phases. The plot shows a strong temporal evolution of the tracks. The same data are shown in Fig. B.1, split in separate time-resolved colour-colour diagrams. Using the KDE method (see Sect. 3.3), we inferred the PDF of the data in each diagram. The resulting probability distribution maps are shown in Fig. 8.

Fig. 7. Time-resolved colour-colour diagrams of observation 201, extracted from segments with a duration of 11 ks (corresponding to ∼0.024 in orbital phase) and time resolution of 10 s. The resulting ten colour-colour diagrams are overplotted using different colours that correspond to different phase intervals, according to the colour map reported on the right.

The passage at superior conjunction happens at t ∼ 86.36 ks after the starting time of the observation (see Fig. 2). Therefore, it is contained in panel “h” corresponding to ϕ_orb = 0.988 − 0.012. The fact that the LOS passes through the densest regions of the stellar wind in this phase of the orbit is clearly seen from the shape of the corresponding colour-colour track, whereby ‘only’ the lower region of the diagram is populated (i.e. at hard colours ≲0.75). A similar distribution is observed both at earlier and later phases of the orbit from panel “f” to panel “i” due to strong dipping events occurring several hours before and after superior conjunction (see Fig. 2).

We built simulated tracks for each time-resolved colour-colour diagram in order to constrain the changes in the wind parameters, which can cause the observed changes in the colour-colour tracks. To this aim, we used the complex continuum model fitted to the spectrum of each time-resolved epoch described in Sect. 3.2 and the variable ionisation absorption component described by the empirical function $log\xi =log\frac{10+[N_{\mathrm{H},\mathrm{w}}/{10}^{22}{\mathrm{cm}}^{-2}]}{[N_{\mathrm{H},\mathrm{w}}/{10}^{22}{\mathrm{cm}}^{-2}]}+1$, introduced in Sect. 3.3 (red curve in the left panel of Fig. 6). Although this function cannot be statistically preferred to the other one tested in Sect. 3.3 (blue curve in the left panel of Fig. 6), our choice is dictated by the willingness to include the effects of additional sources of ionisation becoming important in denser environments (e.g. Feldmeier et al. 1997; Sundqvist & Puls 2018), such as the base of the wind, intercepted by the LOS at superior conjunction. We selected as best-fit tracks those with the highest combined likelihood, as described in Sect. 3.3. The best-fit tracks are shown in Fig. 8 (as well as in Fig. B.1), and the corresponding best-fit covering factors for each time-resolved diagram are listed in Table 1.

Fig. 8. Probability distribution maps of each time-resolved colour-colour diagram obtained using the KDE method. In blue, the best-fit simulated track. The grey-shaded closed curves represent the probability distribution map for orbital phases ϕorb = 0.43 − 0.46 at the 99.7% (in black), 95% (in grey), and 68% (in light grey) confidence levels.

We notice that the best-fit tracks do not always extend to high hard colours as much as the data do. This particularly happens in the less absorbed phases of the orbit when the source is out of superior conjunction (i.e. panels a, b, c, and j). This is because the model reaches the maximum allowed value for the ionisation parameter (in our model, the lowest values of N_H, w correspond to the highest values of ξ, Fig. 6, left panel). The upper portion of the diagram not covered by the wind model track is populated by the least absorbed data bins. Therefore, the scattering in the measured hard and soft colours within this region might be completely driven by the intrinsic spectral variability of the source. To verify this, we calculated a colour-colour diagram from observation 501, which is one of the least affected by the wind as it was carried out right before inferior conjunction (ϕ_orb = 0.5; see Lai et al. 2022). We note that the CHOCBOX campaign did not catch the passage at inferior conjunction, and therefore from this observation we selected the closest orbital phases to it (ϕ_orb = 0.43 − 0.46). Observation 501 is also consecutive to observation 201. This allows us to rely on the implicit assumption that the intrinsic continuum spectrum of the X-ray source has not changed significantly between the two observations, and that any scatter in the colour-colour diagram would only be due to fluctuations in the parameters of the same continuum model. Using the KDE method (Sect. 3.3), we calculated the probability distribution of soft and hard colours in the diagram. The colour-colour diagram of observation 501 and the corresponding contour plots are reported in Appendix C and Fig. C.1. As it can be seen, the data are quite scattered, suggesting intrinsic short-term spectral variability of the continuum. In Figs. 8 and B.1 we overplot the 99.7%, 95% and 68% confidence contours of observation 501, to mark the area in the diagram where the intrinsic emission from the X-ray source dominates. In other words, in this region the probability for the source to be unabsorbed by the stellar wind is the highest. The data points not reached by the wind model tracks are thus fully consistent with being free from wind absorption.

From the fit of the time-resolved colour-colour diagrams we infer a steady and significant increase in the covering factor, by a factor of ∼1.2 between ϕ_orb = 0.820 − 0.844 and ϕ_orb = 0.036 − 0.06 (see Table 1). The observed trend is plotted in Fig. 9 in the upper panel. The maximum value of the covering factor is registered at superior conjunction, but the data suggest this parameter remains high up to at least ϕ_orb = 0.06. From the PDFs of Fig. 8 and our best-fit tracks we extracted the average column density ${\overline{N}}_{\mathrm{H},\mathrm{w}}$ and its scatter δN_H, w, as this parameter is an indicator of the amount of short-term variability of N_H, w within each epoch. These values are reported in Table 1 and plotted in Fig. 9 (middle and lower panels) for each epoch. Our procedure is this: given the best-fit track, we extract the probability of each point from the KDE, which is thus a function of N_H, w. This probability is a proxy for the number of data points with a given N_H, w. Then, we consider the range of N_H, w that encloses 68% of the total (renormalised) probability for the corresponding track to be a proxy for the total number of data points. The width of the range gives us δN_H, w and its average value gives us ${\overline{N}}_{\mathrm{H},\mathrm{w}}$. We reckon this approach to be more informative than, for example, extracting the N_H, w with the highest probability, since the latter method would result in a loss of information, especially for the distributions that show more than one peak. For example, the fit of the colour-colour diagram of panel “a” of Fig. 8 would return a low N_H, w, thus missing the information associated with the second peak of the PDF. Our procedure also provides a straightforward way to obtain estimates of the confidence intervals for each pair of ${\overline{N}}_{\mathrm{H},\mathrm{w}}$ and δN_H, w values. We now consider the two tracks corresponding to the extrema of the 68% confidence level of f_c and repeat the procedure obtaining new confidence level of N_H, w. To be conservative, we take the overall maximal and minimal confidence level and from these obtain the intervals of ${\overline{N}}_{\mathrm{H},\mathrm{w}}$ and δN_H, w.

Fig. 9. Stellar wind parameters as obtained from the fit of the time-resolved colour-colour diagrams of observation 201, plotted as a function of the orbital phase (see also Table 1). The parameters are: the covering factor, fc (upper panel), the mean column density, ${\overline{N}}_{\mathrm{H},\mathrm{w}}$ (middle panel), and the scatter in column density parameter, δNH, w (lower panel). The shaded areas represent the 68% confidence contours. The yellow curves are the neutral column density profiles produced by the smooth wind model (El Mellah et al. 2020) at high (solid) and low (dashed) stellar mass loss rates.

We observe significant variations in both the ${\overline{N}}_{\mathrm{H},\mathrm{w}}$ and the δN_H, w parameters. While the two show the same trend of variations (when the average column density increases, its scatter also increases), their trend differs from that of f_c. In particular, both the ${\overline{N}}_{\mathrm{H},\mathrm{w}}$ and the δN_H, w show two peaks, the second one happening at superior conjunction, when also f_c reaches its maximum value. However, when f_c is still at its maximum, ${\overline{N}}_{\mathrm{H},\mathrm{w}}$ and δN_H, w drop by a factor of ∼2.5 and ∼2, respectively. This is readily visible in panel “j” of Figs. 8 and B.1, where the data populate again the region of high hard colours of the colour-colour diagram (as in panels “a–e”) while also stretching down to the lowest registered hard colours. Finally, it is worth noticing that the measured values of ${\overline{N}}_{\mathrm{H},\mathrm{w}}$ (between ∼6 × 10²² cm⁻² and ∼16 × 10²² cm⁻²) all fall in the regime where our model predicts additional ionisation in high density environments to become significant (see Fig. 6, left panel). In this regime, strong variations in ${\overline{N}}_{\mathrm{H},\mathrm{w}}$ along the dips are not mirrored by equally strong variations in the ionisation parameter, which remains constrained between log ξ ∼ 1.2 − 1.4.

4.1. On the soft tail of colour-colour tracks

Before discussing our results, we note that none of the models used to fit the colour-colour diagrams (Sects. 3 and 4) produces tracks that can explain the extended tail at high values of soft colours characterising the most absorbed stages of the dips (lower part of the colour-colour diagram, see e.g. Figs. 6 and B.1).

Such high-value soft colours indicate the presence of an additional soft component that becomes significant in the deeper parts of the dips. The most likely candidates for explaining this behaviour are the contribution from emission lines produced by the wind material, and/or from a dust-scattering halo (e.g. Nowak et al. 2011). The former possibility is supported by results presented in Hirsch et al. (2019), who revealed the presence of an emission line spectrum from the photoionised plasma around the BH, emerging only when the X-ray source is highly absorbed (down in the dips). Additional contribution from a dust-scattering halo cannot be excluded. Indeed, source photons intercepted by foreground interstellar dust on the LOS are scattered away, thus dimming the source, but interstellar dust grains outside the LOS can redirect photons back to the observer, producing a halo around the source. Given the energy dependence of the dust-scattering cross-section (Corrales et al. 2016) soft X-ray photons will have more chance to be scattered back, producing an excess of soft X-ray flux in the halo (Maeda et al. 1996).

In order to check for the presence of a dust-scattering halo component in the XMM-Newton EPIC-pn data of Cyg X-1, we verified whether the spectrum significantly softens when extracting photons from increasingly more external regions, as described in Jin et al. (2017). It is worth noticing that the timing mode does not allow us to select regions totally unaffected by a dust-scattering halo since all pixels in the same column lose their spatial information along RAWY. Nonetheless, if a scattering halo is present and the X-ray source is bright, the inner columns will still be harder than the external columns due to the X-ray source dominating the central columns and the halo dominating the external ones.

We extracted spectra in the 2 − 10 keV band during the time interval ∼82 − 94 ks that comprises the passage at superior conjunction (see Fig. 10). Indeed, the emission from the halo is expected to be more significant when the X-ray source is highly absorbed (e.g. Jin et al. 2018), and thus in the deepest stages of the dips. We selected three detector regions at increasingly higher distance from the central pixels (and discarding the central RAWX = [36:39] to avoid pile-up effects, Sect. 2): RAWX = [32:35]−[40:43] covering an angular diameter D = 8″ − 16″ from the central pixel, RAWX = [30:32]−[43:46] (D = 16″ − 32″) and RAWX = [27:30]−[46:49] (D = 32″ − 44″). A gradual softening of the Γ parameter is observed, up to ∼10% in the most external region. In particular, Γ = 1.42 ± 0.01 for RAWX = [32:35]−[40:43], Γ = 1.54 ± 0.01 for RAWX = [30:32]−[43:46], and Γ = 1.58 ± 0.01 for RAWX = [27:30]−[46:49]. This gradual softening hints at the presence of a dust-scattering halo component. We conclude that the dust-scattering halo can potentially contribute to the soft tail observed in the colour-colour tracks. Nonetheless, we notice that both a dust-scattering halo and the emission line component from the wind would contribute the most when the source is highly absorbed. Thus, we expect a soft tail to be particularly prominent when the uppermost part of the colour-colour diagram is not populated. However, all the observed time-resolved colour-colour diagrams follow these expectations (more visible in Fig. B.1), except the one at ϕ_orb = 0.012 − 0.036 (see panel “i” of Fig. B.1), which does not show a prominent soft tail despite being highly absorbed (${\overline{N}}_{\mathrm{H},\mathrm{w}}\sim 10.80\times {10}^{22}{\mathrm{cm}}^{-2}$ and f_c ∼ 0.92). We suggest this might be due to a delayed response of the illuminated dust halo or emission line region to a change of the irradiating flux from the X-ray source (see Fig. 2).

Fig. 10. Unfolded (against a constant model fixed at 1) 2 − 10 keV EPIC-pn spectra (upper panel) extracted in the time interval t ∼ 82 − 94 ks, catching the passage at superior conjunction. The chosen RAWX selections for each spectrum are reported in the labels. The softening of the spectra suggests that there is contribution from a dust-scattering halo, easily noticeable in the ratios of the unfolded spectra (bottom panel).

5. Discussion

Several studies have demonstrated the complexity of the stellar wind in Cyg X-1 (e.g. Grinberg et al. 2015; Miškovičová et al. 2016; Hirsch et al. 2019; Grinberg et al. 2020). The wind appears to be formed of over dense clumps (e.g. Owocki et al. 1988; Feldmeier et al. 1997; Oskinova et al. 2012; Sundqvist & Owocki 2013) crossing our LOS to the X-ray source, producing stochastic variability that typically manifests itself as absorption dips in the X-ray light curves, and adds up to the intrinsic variability produced in the inner accretion flow (Lai et al. 2022). Absorption dips occurrence and concurrent absorption variability are more prominent near superior conjunction (e.g. Bałucińska-Church et al. 2000; Lai et al. 2022), when the companion star is between the observer and the BH (ϕ_orb = 0), as in these phases the LOS crosses deeper wind layers.

In this work, we have presented a characterisation of the stellar wind properties in Cyg X-1, through the analysis and modelling of the colour-colour diagrams of the source during a passage at superior conjunction (between orbital phases 0.82 and 0.06, see Figs. 1 and 2), with the X-ray source in its hard state. The prominence of absorption dips in this specific state, rather than in softer states, are explained as due to the X-ray source not being strong enough to fully ionise the wind (e.g. Miškovičová et al. 2016; Grinberg et al. 2015). In addition to previous studies, we employed a non-parametric fitting method based on the kernel density estimation analysis (Hill 1985), which allowed us to select the model that best describes the observed colour-colour diagrams and constrain the physical parameters of the stellar wind (Sects. 3.3 and 4).

5.1. Modelling of colour-colour diagrams and evolution of stellar wind parameters

We confirmed (see also Grinberg et al. 2020 who analysed a Chandra 2004 observation of the source covering similar orbital phases) that the characteristic ‘pointy’ or ‘nose-like’ shape of the colour-colour tracks requires the absorbing gas to be partially ionised (Sect. 3.2), in agreement with independent results from high-resolution spectroscopy (e.g. Hanke et al. 2009; Miškovičová et al. 2016). In particular, Hirsch et al. (2019) revealed different levels of ionisation at different absorption stages. This suggests an (unknown) functional dependence of log ξ on N_H, w (Grinberg et al. 2020). We tested two empirical functions (Sect. 3.3), one assuming a simple linear scaling and the other allowing for additional sources of ionisation at high column densities (Fig. 6). Although the data do not allow us to statistically prefer either of the two models, we notice that simulations show that collisions in the ambient stellar wind occurring in denser regions (high N_H, w regime) are expected to produce significant X-ray emission (Feldmeier et al. 1997; Sundqvist & Puls 2018). This emission can contribute to further ionise the gas.

In order to characterise the evolution of the wind as a function of the orbital phase around superior conjunction, we extracted time-resolved colour-colour diagrams, and fit them with a continuum plus variable ionisation model (see Sect. 4). The fits show a steady increase in the covering factor, f_c, by a factor of ∼1.2 between ϕ_orb ∼ 0.8 and ϕ_orb ∼ 0, which remains high also after superior conjunction (up to at least ϕ_orb ∼ 0.06, Fig. 9 top panel and Table 1). On the other hand, the mean column density ${\overline{N}}_{\mathrm{H},\mathrm{w}}$ shows two peaks, one at ϕ_orb ∼ 0.9 and the other at superior conjunction. After superior conjunction, while f_c is still high, ${\overline{N}}_{\mathrm{H},\mathrm{w}}$ drops to a minimum (Fig. 9 middle and Table 1). We also measured the scatter parameter, δN_H, w, which quantifies the amount of variability of N_H, w within each time-resolved colour-colour diagram (Fig. 9 bottom panel and Table 1). This parameter follows the same trend displayed by ${\overline{N}}_{\mathrm{H},\mathrm{w}}$ (see also Fig. 11). Given the sampling adopted for the extraction of time-resolved colour-colour diagrams (10 segments of ∼11 ks), while the inferred ${\overline{N}}_{\mathrm{H},\mathrm{w}}$ parameter probes relatively long-term modulations, the δN_H, w parameter probes more rapid variability (between 10 s, the time bins used to build the colour-colour diagrams, and ∼11 ks) likely driven by the smallest-scale inhomogeneities in the wind. Our study hints at a one-to-one relation between the amount of absorption at a given phase and its associated rapid variability: the higher the ${\overline{N}}_{\mathrm{H},\mathrm{w}}$, the higher its scatter δN_H, w (see Fig. 11). This is reminiscent of the “rms-flux” relation characterising the stochastic flux variability in compact object systems (e.g. Uttley et al. 2005; Scaringi et al. 2012), which implies a log-normal distribution of the flux and requires variability on different timescales to combine multiplicatively (Uttley et al. 2005). Since the size of the clumps influences the observed column density, confirmation of such correlation might have implications on the distribution of clump sizes and on the way they combine to form bigger clumps.

Fig. 11. Relation between the measured ${\overline{N}}_{\mathrm{H},\mathrm{w}}$ and δNH, w parameters, respectively sampling long-term and rapid variations in the column density parameter of the stellar wind (the values are listed in Table 1).

5.2. Stellar wind mass loss rate and clump mass estimations

The variations in the absorbing column density stem from two aspects. First, we expect a periodic variation due to the orbital motion of the BH around its stellar companion. This smooth component can be described by the absorption from a smooth wind, but it falls short at capturing the stochastic variability on timescales much shorter than the orbital period. The latter comes from clumps intercepting the LOS between the observer and the BH (El Mellah et al. 2020). The average column density time series match the smooth wind profile of same mass loss rate, while the dispersion above and below tell us about the clumps’ properties.

The data considered in this paper cover only a small portion of the orbital period, so they cannot be used to constrain the periodic variations in column density. Instead, we plotted in the middle panel in Fig. 9, for two different stellar mass loss rates, the periodic component of the column density due to the orbital motion (yellow lines). It was computed based on the smooth and spherically symmetric stellar wind model of El Mellah et al. (2020): it accounts for all the material along the LOS to the BH, regardless of its ionisation state. Therefore, it is an estimate of the neutral column density, necessarily higher than the measurable ionised column density. If we neglect the orbital eccentricity and assume a β-law for the wind velocity profile, the shape of the neutral column density curve only depends on the orbital inclination i, on the ratio of the orbital separation d to the stellar radius R_* and on the β exponent of the velocity profile. We take i = 27° and d/R_* = 2.3, which correspond to the values given by Miller-Jones et al. (2021), and β = 1.5, a representative value for UV-line driven winds of O supergiants (Rubio-Díez et al. 2022). Owing to the low inclination of Cyg X-1, the peak-to-peak variation remains modest across the orbit and no significant changes are expected for realistic values of d/R_* and β. On the other hand, the scale of the neutral column density curve is set by ${\dot{M}}_{*}/(m_{\mathrm{p}}{v}_{\infty }{R}_{*})$, where Ṁ_* is the stellar mass loss rate and m_p is the proton mass. We use the values reported by Miller-Jones et al. (2021) to set v_∞ = 2100 km/s and R_* = 22 R_⊙. Through Hα diagnostics, Gies et al. (2003) reported a stellar mass loss rate ${\dot{M}}_{*}\sim 2.6\times {10}^{-6}{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}$ in the hard state. However, they do not take ionisation into account, and they rely on the stellar parameters used by Herrero et al. (1995), which do not correspond to the updated values. Therefore, there is a risk that this mass-loss rate is underestimated. And indeed, when plotted over our data points, the neutral column density profile with ${\dot{M}}_{*}\sim 2.6\times {10}^{-6}{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}$ (dashed yellow line) lies below the data points. However, since the profiles, in Fig. 9, produced by the smooth wind model of El Mellah et al. (2020) represent the neutral column density while the data points stand for the absorbing column density (i.e. only the material at ionisation levels low enough to contribute to absorption), the smooth profile should be an upper limit. Therefore, we think a higher stellar mass loss rate, for instance ${\dot{M}}_{*}\sim 7\times {10}^{-6}{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}$ (solid yellow line), is more compatible with the column density we measure at superior conjunction. Measurements of ${\overline{N}}_{\mathrm{H},\mathrm{w}}$ over multiple consecutive orbits are needed to obtain better constraints.

Moreover, both the ${\overline{N}}_{\mathrm{H},\mathrm{w}}$ and the δN_H, w parameters show significant short-term variability (Fig. 9). El Mellah et al. (2020) presented a model that explores the impact of the wind clumpiness on the rapid time-variability of its column density, and describes how this observed variability can be linked to the clumps’ physical properties, such as the size and the mass. In the clumpy wind model of El Mellah et al. (2020), the clumps’ properties are directly linked to the scatter parameter, δN_H, w. This parameter is expected to be maximum around superior conjunction, as indeed observed (Table 1 and Fig. 9, lower panel). The model also predicts the scatter to be an excellent tracer of the ratio $\sqrt{m_{\mathrm{cl}}}/R_{\mathrm{cl}}$, where m_cl and R_cl are the mass and radius of the clumps. In Lai et al. (2022) we measured the minimum timescales at which the passage of clumps causes significant excess X-ray variability in the power spectra of the source at superior conjunction. These timescales correspond to a clump radial size of ≳10⁻⁴ R_*. Using this as the smallest clumps radius, the measured maximum value of δN_H, w (Table 1), and the above estimate for the mass loss rate, and plugging them into Equation (18) of El Mellah et al. (2020) we infer m_cl ∼ 10¹⁷ g as an estimate for the characteristic mass of the clumps (for v_∞ = 2100 km/s and R_* = 22 R_⊙). This estimate is in agreement with values inferred from radiative hydrodynamical simulations (Sundqvist et al. 2018), as well as from spectral analysis (Härer et al. 2023).

5.3. The soft colour tail

The observed colour-colour tracks also reveal the presence of a soft emission component in the most absorbed stages, producing an extended soft-colour tail. We verified (Sect. 4) that the analysed XMM-Newton data are consistent with the presence of significant contribution from a dust-scattering halo (e.g. Jin et al. 2017). However, additional contribution from the wind itself in the form of emission lines from the diffused photoionised gas around the BH is also expected (Hirsch et al. 2019). Indeed, both contributions from a dust-scattering halo and the photoionised emitting plasma are expected to be prominent when the primary radiation from the X-ray source is significantly blocked by the absorbing material; that is, in the deepest stages of the dips. A more complex modelling (which is beyond the scope of this paper) would allow the two contributions to be disentangled, and might in turn be useful to put stronger constraints on wind models, as well as probing the interstellar dust towards Cyg X-1.

As a final remark, it is worth noting that additional spectral changes extrinsic to the wind – intrinsic changes of the hard X-ray source properties – have been neglected in our fits of time-resolved colour-colour diagrams. For example, such changes might be due to time variations in the spectral index, Γ, as induced by variations in either the temperature or the optical depth of the Comptonising gas (or both). In Sect. 4, we verified that the scattering in hard and soft colours in the uppermost, unabsorbed part of the diagrams are most likely driven by spectral variability. Mastroserio et al. (2021) showed that Γ variability of a few percents is plausible in hard state sources with ∼10% fractional rms variability, typical of the hard state. Skipper et al. (2013) found a larger Γ variability (∼10%) for the same source and spectral state but how these spectral changes influence the most absorbed parts of the diagram is yet to be tested.

6. Conclusions

In this paper, we have presented an analysis of the colour-colour diagrams of Cyg X-1 during a passage at superior conjunction, with the aim of constraining the properties of the stellar wind. We employed the KDE method to select wind models that best fit the colour-colour diagrams. This allowed us to overcome the problem of dealing with data that are not normally distributed. The main results are the following:

• We found that the model that best describes the characteristic ‘pointy’ or ‘nose-like’ shape of the diagrams implies the wind to be partially ionised.

• We revealed a strong temporal evolution of the colour-colour diagrams around superior conjunction. Our fits suggest this evolution to be strongly influenced by concurrent variations in the column density and covering factor of the wind.

• Both the column density and the covering factor are maximum at superior conjunction (since the LOS crosses deeper wind layers), but their overall variations follow different trends.

• We report a one-to-one scaling between long-term (> 11 ks) and rapid (between 10 s and 11 ks) variations in the column density (reminiscent of the ‘rms-flux’ relation characterising stochastic variability in compact objects). The existence of such a correlation might have implications for the distribution of the wind clumps and for the way they combine into bigger clumps.

• Using the clumpy wind model proposed in El Mellah et al. (2020), we have estimated a wind mass loss rate of ${\dot{M}}_{*}\sim 7\times {10}^{-6}{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}$ and a characteristic clump mass of m_cl ∼ 10¹⁷ g.

Future applications of this analysis approach will require the sampling of longer stretches of the binary orbit and the application of more complex models (e.g. including contribution from scattered wind emission and from the dust-scattering halo).

Acknowledgments

This work is based on observations obtained with XMM-Newton, an ESA science mission instrument, and contributions directly funded by ESA Member States and NASA. EVL thanks Alex Markowitz, Piotr Życki and Maura Pilia for useful discussions. The research leading to these results has received funding from the European Union’s Horizon 2020 Programme under the AHEAD2020 project (grant agreement n. 871158). EVL and AR were partially supported by the Polish National Science Centre, grant no. 2021/41/B/ST9/04110. EVL and MB are supported by the Italian Research Center on High Performance Computing Big Data and Quantum Computing (ICSC), project funded by European Union – NextGenerationEU – and National Recovery and Resilience Plan (NRRP) – Mission 4 Component 2 within the activities of Spoke 3 (Astrophysics and Cosmos Observations). BDM acknowledges support via Ramón y Cajal Fellowship (RYC2018-025950-I) and the Spanish MINECO grant PID2022-136828NB-C44. YC acknowledges support from the grant RYC2021-032718-I, financed by MCIN/AEI/10.13039/501100011033 and the European Union NextGenerationEU/PRTR. BDM, YC, GS, and JJ thank the Spanish MINECO grant PID2020-117252GB-I00 and the AGAUR/Generalitat de Catalunya grant SGR-386/2021. MC acknowledges support from the “Universitas Copernicana Thoruniensis In Futuro” project nr POWR.03.05.00-00-Z302/17 and from the MNiSW grant DIR/WK/2018/12. AAZ acknowledges support from the Polish National Science Center grants 2019/35/B/ST9/03944 and 2023/48/Q/ST9/00138, and from the Copernicus Academy grant CBMK/01/24. MB was funded in part by PRIN TEC INAF 2019 “SpecTemPolar! – Timing analysis in the era of high-throughput photon detectors”. This research has made use of NASA’s Astrophysics Data System Bibliographic Service (ADS) and of ISIS functions (isisscripts http://www.sternwarte.uni-erlangen.de/isis/) provided by ECAP/Remeis observatory and MIT. The authors thank Mirjam Oertel, developer of the first versions of some of the scripts calculating the colour-colour tracks. They also thank John E. Davis for the development of the slxfig (http://www.jedsoft.org/fun/slxfig/) module used to prepare some of the figures in this work. Colour schemes used in the first part of this paper were based on Paul Tol’s color-blindness-friendly palettes and templates (https://personal.sron.nl/~pault/).

References

Appendix

The model used to describe the broad band primary continuum for simulations reported in Sects. 3.2 and 3.3 includes a disc black body component, a soft Comptonisation component, and a hard Comptonisation component and its associated reflection, that is, TBabs × [diskbb + nthComp + relxillCp] in XSPEC. We froze the BH spin at the maximum value a = 0.998 (which allows the inner disc truncation radius to span the largest range of values), the inclination of the reflector i at 27° (Miller-Jones et al. 2021, see also recent results by Krawczynski et al. 2022 and Poutanen et al. 2023) and the high-energy cut-off of the hard Comptonisation component at 100 keV (Basak et al. 2017). The nthComp component is included in the model to account for the observed “soft excess”, as done in e.g. Basak et al. (2017). The seed photon temperature of the soft Comptonisation component is tied to the best-fit inner disc temperature of the diskbb component.

This model was fit to the spectrum of the XMM-Newton observation 201 of Cyg X-1 after filtering out strong dips (the NWA dataset selected in Lai et al. 2022). After fitting this model, the data still show narrow residuals, mainly between ∼0.5 − 2.5 keV. Since these residuals appear close to the absorption edges of the response matrix, i.e. at E ≈ 0.528 keV, and between E ≈ 1.83 − 1.87 keV, they are likely due to calibration problems. In order to account for these residuals, we modified the response matrix using the gain function in XSPEC. This command permits to change the slope or the intercept of the effective area curve, shifting the energies at which the response matrix is defined. Therefore, we introduced a linear gain shift (intercept parameter in XSPEC) of 0.01 keV (slope parameter of gain fixed to 0). As a result, the best-fit strongly improves (from χ²/dof = 4404.17/1891 to χ²/dof = 3923.7/1891, see Fig. A.1).

Fig. A.1. Comparison of the data to model ratios before (upper panel) and after (bottom panel) using the gain correction (intercept parameter set to 0.01 keV and slope parameter fixed to 0).

Fig. A.2. The best-fit model of the NWA time-averaged spectrum of observation 201. The complex best-fit model is overplotted in black, while the single components are shown in different colours: in magenta the disc blackbody, in blue the soft excess, in green the hard Comptonisation component and its relative reflection and in dotted grey the additional gaussians. Ratios of the data to the best-fitting model are shown in the bottom panel.

Nonetheless, some residuals were still present, including a strong narrow excess in the Fe Kα region. To account for these features, we added narrow (σ < 0.1 keV) emission and absorption gaussian components (see Tab. A.1). We did not investigate the nature of these features, but we point out that they might be due to incorrect calibration or residual wind absorption after the selection process of the NWA GTIs (Lai et al. 2022). On the other hand, the narrow feature observed in the Fe Kα region might indicate the need to include a second reflection component (such as from the outer disc) to obtain a better description of the broad band continuum. Nonetheless, we stress that our analysis, which focusses on testing the effects of variable absorption on the colour-colour tracks, is not strongly influenced by the details of the continuum model as long as it provides a good description of the broad-band spectrum.

The final best-fit model yields χ²/dof = 2263.76/1873. Fig. A.2 shows the spectrum of the NWA time-averaged spectrum of observation 201 and its best-fit model. The most relevant parameters are listed in Tab. A.1.

Appendix B: Time-resolved colour-colour diagrams: best-fit models

In Fig. B.1 we show the time-resolved colour-colour diagrams of observation 201 (colour-coded as in Fig. 7), and the best-fit models obtained from the fits discussed in Sect. 4.

Fig. B.1. Time-resolved colour-colour diagrams of observation 201. Overplotted (black tracks) are the best-fit models obtained assuming $log\xi =log\frac{10+[N_{\mathrm{H},\mathrm{w}}/{10}^{22}{\mathrm{cm}}^{-2}]}{[N_{\mathrm{H},\mathrm{w}}/{10}^{22}{\mathrm{cm}}^{-2}]}+1$. The contour plots represent the 99.7% (in black), 95% (in grey) and 68% (in light grey) confidence levels of the data distribution during the least absorbed stages of the orbit (ϕorb = 0.43 − 0.46, see Appendix C).

Appendix C: Probability distribution close to inferior conjunction

We show the colour-colour diagram calculated from orbital phases ϕ_orb = 0.43 − 0.46, covered by observation 501 (Lai et al. 2022). These phases are selected to be the closest to inferior conjunction (i.e. ϕ_orb ∼ 0.5), and thus the least affected by wind absorption. Using the KDE method, we computed the probability distribution for this dataset (Sect. 4). Fig. C.1 shows the resulting 99.7%, 95% and 68% confidence contours overplotted to the data.

Fig. C.1. Colour-colour diagram and probability distribution map for orbital phases ϕorb = 0.43 − 0.46. The curves correspond to the 99.7% (in black), 95% (in grey) and 68% (in light grey) confidence regions.

All Tables

All Figures

Fig. 1. Orbital phase coverage of the XMM-Newton CHOCBOX monitoring of Cyg X-1 in the hard state. The four arc labels correspond to the ObsIDs of each pointing (see Table 1 in Lai et al. 2022 for the full ObsID). The orbital phase, ϕorb = 0, indicates the passage at superior conjunction. Observation 201 (in purple) is the focus of the analyses presented in this paper. The starting orbital phase of the monitoring (corresponding to ϕorb = 0.82) is also indicated. HDE 226868 and Cyg X-1 are only schematically represented, no Roche lobe was graphically considered in the scheme.

In the text

Fig. 2. XMM–Newton EPIC-pn light curves of the observation 201 of Cyg X-1 with time bins of 10 s. The upper and the middle panels show, respectively, the 0.5 − 1.5 keV and the 3 − 10 keV light curves. The bottom panel reports the hardness, i.e. the ratio between count rates in the 3 − 10 keV and 0.5 − 1.5 keV energy bands. The vertical red line indicates the passage at superior conjunction (i.e. ϕorb = 0).

In the text

	Fig. 3. Simulated tracks for a power law (with Γ fixed at 1.6, left panel, 1.7, middle panel, and 1.8, right panel) plus a neutral absorber, compared to the observed (grey points) colour-colour diagram tracks of observation 201 of Cyg X-1. The blue (red) tracks describe changes in the wind fc (NH, w) for a fixed value of NH, w (fc).
In the text

	Fig. 4. Simulated tracks for the complex continuum absorbed by a neutral wind model, as compared to the colour-colour diagram of Cyg X-1.
In the text

	Fig. 5. Simulated tracks for a homogeneously ionised absorber (with log ξ values ranging between −4 and 2) compared to the colour-colour diagram of observation 201. The different panels show the simulated tracks computed for different values of fc: 0.8 (left), 0.87 (middle), 0.95 (right). The black curve in the middle panel represents the one that best describes the data.
In the text

Fig. 6. Empirical functions describing the dependence of the ionisation parameter on the column density of the absorbing gas. Left panel: Assumed dependencies of the ionisation parameter, log ξ, on the column density, NH, w. Middle and right panel: Simulated tracks for a warm absorbing gas, assuming log ξ = log(100/[NH, w/1022 cm−2]) (middle panel), and $log\xi =log\frac{10+[N_{\mathrm{H},\mathrm{w}}/{10}^{22}{\mathrm{cm}}^{-2}]}{[N_{\mathrm{H},\mathrm{w}}/{10}^{22}{\mathrm{cm}}^{-2}]}+1$ (right panel). The different shades correspond to different values of covering factor as reported in the labels.

In the text

	Fig. 7. Time-resolved colour-colour diagrams of observation 201, extracted from segments with a duration of 11 ks (corresponding to ∼0.024 in orbital phase) and time resolution of 10 s. The resulting ten colour-colour diagrams are overplotted using different colours that correspond to different phase intervals, according to the colour map reported on the right.
In the text

	Fig. 8. Probability distribution maps of each time-resolved colour-colour diagram obtained using the KDE method. In blue, the best-fit simulated track. The grey-shaded closed curves represent the probability distribution map for orbital phases ϕorb = 0.43 − 0.46 at the 99.7% (in black), 95% (in grey), and 68% (in light grey) confidence levels.
In the text

Fig. 9. Stellar wind parameters as obtained from the fit of the time-resolved colour-colour diagrams of observation 201, plotted as a function of the orbital phase (see also Table 1). The parameters are: the covering factor, fc (upper panel), the mean column density, ${\overline{N}}_{\mathrm{H},\mathrm{w}}$ (middle panel), and the scatter in column density parameter, δNH, w (lower panel). The shaded areas represent the 68% confidence contours. The yellow curves are the neutral column density profiles produced by the smooth wind model (El Mellah et al. 2020) at high (solid) and low (dashed) stellar mass loss rates.

In the text

Fig. 10. Unfolded (against a constant model fixed at 1) 2 − 10 keV EPIC-pn spectra (upper panel) extracted in the time interval t ∼ 82 − 94 ks, catching the passage at superior conjunction. The chosen RAWX selections for each spectrum are reported in the labels. The softening of the spectra suggests that there is contribution from a dust-scattering halo, easily noticeable in the ratios of the unfolded spectra (bottom panel).

In the text

	Fig. 11. Relation between the measured ${\overline{N}}_{\mathrm{H},\mathrm{w}}$ and δNH, w parameters, respectively sampling long-term and rapid variations in the column density parameter of the stellar wind (the values are listed in Table 1).
In the text

	Fig. A.1. Comparison of the data to model ratios before (upper panel) and after (bottom panel) using the gain correction (intercept parameter set to 0.01 keV and slope parameter fixed to 0).
In the text

Fig. A.2. The best-fit model of the NWA time-averaged spectrum of observation 201. The complex best-fit model is overplotted in black, while the single components are shown in different colours: in magenta the disc blackbody, in blue the soft excess, in green the hard Comptonisation component and its relative reflection and in dotted grey the additional gaussians. Ratios of the data to the best-fitting model are shown in the bottom panel.

In the text

Fig. B.1. Time-resolved colour-colour diagrams of observation 201. Overplotted (black tracks) are the best-fit models obtained assuming $log\xi =log\frac{10+[N_{\mathrm{H},\mathrm{w}}/{10}^{22}{\mathrm{cm}}^{-2}]}{[N_{\mathrm{H},\mathrm{w}}/{10}^{22}{\mathrm{cm}}^{-2}]}+1$. The contour plots represent the 99.7% (in black), 95% (in grey) and 68% (in light grey) confidence levels of the data distribution during the least absorbed stages of the orbit (ϕorb = 0.43 − 0.46, see Appendix C).

In the text

	Fig. C.1. Colour-colour diagram and probability distribution map for orbital phases ϕorb = 0.43 − 0.46. The curves correspond to the 99.7% (in black), 95% (in grey) and 68% (in light grey) confidence regions.
In the text