© ESO 2020

1 Introduction

Massive stars live a forceful life during which they shape galaxies with their mechanical, radiative, and chemical feedback. Evolutionary pathways of massive stars are largely controlled by their line-driven winds, with typical mass-loss rates of 10⁻⁶ M_⊙ yr⁻¹ and terminal wind speeds on the order of 1000 km s⁻¹ (Puls et al. 2008). The flow is provided with net outward momentum thanks to the resonant line absorption of ultraviolet (UV) stellar photons by partly ionized metal ions present in the outer envelope of the star. As the wind accelerates, it keeps tapping higher energy photons thanks to the Doppler shift (Lucy & Solomon 1970; Castor et al. 1975).

Originally, these line-driven winds were treated as smooth, steady-state outflows without any internal shocks or overdensity. However, it has been know for a long time now that they are, in fact, highly structured and variable on a broad range of spatial and temporal scales (for overviews, see Hamann et al. 2008; Sundqvist et al. 2012a; Puls et al. 2015). Indeed, the dependence of the acceleration on the velocity profile itself causes the emergence of a self-growing feedback loop, the line-deshadowing instability, which is responsible for internal shocks and the formation of overdense regions or “clumps” (Owocki & Rybicki 1984; Owocki et al. 1988; Feldmeier 1995). This clumping can then severely affect the radiative transfer computations used to derive synthetic observables; for example, neglecting wind clumping leads to empirically inferred mass-loss rates that may differ by an order of magnitude for the same star, depending on which observational diagnostic is used to estimate the mass loss (Fullerton et al. 2006). Typically, mass-loss rates derived from diagnostics with opacities depending on the square of the density (e.g., the Hα line and thermal radio emission) are overestimated if clumping is neglected, whereas those derived from diagnostics depending linearly on density (usually UV resonance lines) can be underestimated if the light-leakage effects of a wind that is porous in velocity-space are not properly taken into account (e.g., Sundqvist & Puls 2018). X-ray emission line shapes from single stars offer an alternative mass-loss rate diagnostic based on continuum photoelectric opacity which is linear in density and insensitive to velocity field inhomogeneity (Owocki & Cohen 2001; Cohen et al. 2014). This diagnostic might still be compromised by geometrical porosity (Feldmeier et al. 2003; Sundqvist et al. 2012b), although observational constraints disfavor the extreme degree of wind structure that would be needed for this (Leutenegger et al. 2013). Due to the faintness of X-ray emission from single stars and the modest size of existing X-ray telescopes, this diagnostic is currently limited to a handful of nearby objects.

High-mass X-ray binaries (HMXBs) can shed a unique light upon the mass loss and micro-structure of line-driven winds. In Supergiant X-ray binaries (SgXBs), a subclass of HMXBs, a compact object, generally a neutron star (NS), sometimes a black hole (BH), orbits a blue supergiant and accretes part of its winds (Walter et al. 2015; Martínez-Núñez et al. 2017). The material captured heats up as it falls into the gravitational potential well and emits copious amounts of X-rays from the immediate vicinity of the compact object. Two methods have been proposed in order to use the X-ray point source as a probe and bring constrains on the wind clumpiness independent from the ones obtained with diagnostics for isolated massive star (Martínez-Núñez et al. 2017).

In the first method, if the episodic enhancements of the intrinsic X-ray luminosity observed in SgXBs are due to the directcapture of a clump, the amount of mass in a clump and the number of clumps per unit volume could be estimated (in ’t Zand 2005). This idea was applied both in persistent SgXBs (Ducci et al. 2009, also called classic SgXBs, with low intrinsic variability,) and in supergiant fast X-ray transients (SFXTs, Bozzo et al. 2017a). In a 1D framework, Oskinova et al. (2012) have calculated synthetic X-ray lightcurves that would result from this approach and compared the expected variability to the observed. This method presupposes that the variability induced by the wind micro-structure remains unaltered from the orbital scale all the way down to the X-ray emitting region (either the NS surface or the inner rim of the accretion disk surrounding the BH). Under this assumption, Oskinova et al. (2012) obtained a peak-to-peak amplitude of the X-ray luminosity several orders of magnitude higher than what observed in SgXBs. However, additional structures form, such as a hydrodynamics bow shock around the accretor, which have been shown to lower the variability compared to a purely ballistic model of the clump capture (El Mellah et al. 2018). The amountof material involved in an observed flare is thus an upper limit on the mass of a clump: a flare might be triggered by the capture of a clump, as noticed by Bozzo et al. (2015), but likely involves additional material which has been previously piling up at the outer rim of the NS magnetosphere due to the magneto-centrifugal gating mechanism (or propeller effect, Illarionov & Sunyaev 1975; Grebenev & Sunyaev 2007; Bozzo et al. 2008), in a quasi-spherical hot shell on top of the NS magnetosphere (Shakura et al. 2013) or in a disk-like structure (El Mellah et al. 2019; Karino et al. 2019). Also, Hainich et al. (2020) showed that there was no systematic difference between the stellar properties in classic SgXBs and in SFXTs, hinting at similar wind micro-structures (Driessen et al. 2019). Finally, Pradhan et al. (2018) concluded that the higher variability of the emission in SFXTs compared to SgXBs was likely due to mechanisms inhibiting accretion near the accretor, rather than direct clump capture, although flares in SFXTs have recently been found by Ferrigno et al. (2020) to be associated with massive structures approaching the accretor. To summarize, X-ray variability is, at best, an indirect tracer of the clumpiness of the wind because of intermediate regions where the clumps are mixed and where other instabilities can kick in, in particular those associated to the different accretion regimes.

In a second method, one can make use of the highly variable absorption local to the HMXBs. If the variations in the absorbing column density are due to unaccreted clumps passing through the line-of-sight (LOS), the accretor could play the role of a steady backlight whose shimmering betrays the stochastic properties of the wind it is embedded into (Oskinova et al. 2012; Grinberg et al. 2015). Time resolved spectroscopy enables us to monitor the absorbing column density directly (e.g., Haberl et al. 1989; Martínez-Núñez et al. 2014; Grinberg et al. 2015). Shorter time scales may be indirectly accessible via changes in X-ray color (e.g., Grinberg et al. 2017; Bozzo et al. 2017a), especially when using multiple energy bands and working with color–color diagrams (Nowak et al. 2011; Hirsch et al. 2019).

Grinberg et al. (2015) laid the foundations of a model of the clumps to relate the observed dips in the X-ray light curve of the HMXB Cygnus X-1 to the micro-structure of the wind. They assumed spherical clumps emitted by the donor star on radial trajectories. However, their exploration was driven by the measurements available using the particular set of observations discussed in their work and limited to the well known parameters of one particular system whose observable properties their model aimed to reproduce. In parallel, new insights have recently been provided on the clump properties thanks to multidimensional radiative hydrodynamics simulations (Sundqvist et al. 2018) and on the velocity profiles thanks to 1D radiative transfer calculations in a steadily expanding atmosphere in nonlocal thermodynamical equilibrium (Sander et al. 2017a). These results drove us into expanding on the first ansatz introduced by Grinberg et al. (2015). We aim at designing a minimal representation of a clumpy wind suitable to grasp the complexity of the flow while still including analytically derivable predictions to prove its robustness. We want to determine whether, in realistic conditions, we can use the median value, the amplitude variations and the coherence time scales of the column density as a function of the orbital phase in order to draw conclusions on the wind micro-structure, in spite of the large number of clumps and of the diversity of clump properties along the LOS. To do so, we carry out a parametric study based on physically-motivated ranges for parameters controlling the orbit, the wind velocity profile and the mass and radius of the clumps. We discuss the observational consequences to be expected according to this theoretical framework but we postpone to later a detailed conversation about the instrumental capabilities and limitations of the current and upcoming X-ray satellites. Also, we focus on the impact of the clumps and discard other sources of column density variability such as larger scale structures produced by the orbital motion for instance.

In Sect. 2, we describe the ingredients of the model we developed and the numerical procedure we follow to compute the column density produced by clumps in number high enough to match the expected number of clumps around massive stars in HMXBs. In Sect. 3, we present the analytic derivation of the column density variance that we confront to the numerical results in Sect. 4, where we also report on the characteristic time scales of absorption episodes and the dependences with respect to the parameters of the model. In Sect. 5, we present how these results can be used to constrain the wind micro-structure from observations, while the limitations of our framework are acknowledged in Sect. 6. Finally, we summarize our findings in Sect. 7.

2 Model

2.1 Wind velocity profile

Hot stars provide their outer envelopes with plenty of high energy photons able to drive a wind. Most of the predictions from the seminal papers on the theory of radiation-driven winds by Lucy & Solomon (1970) and Castor et al. (1975) have been qualitatively confirmed by observations of isolated hot stars, on the main sequence but also evolved (Puls et al. 2006). In the Sobolev approximation, the radial optical depth can be written as a local quantity modulated by the radial velocity gradient (Sobolev 1960). The wind velocity profile can then be determined analytically and it has been common since then to prescribe, for the dynamics of the wind, a radial velocity profile which obeys the β-law, a generalization of the first analytic formula derived (Puls et al. 2008). Its two parameters are the terminal wind speed v_∞, highly supersonic, and a β exponent which quantifies how quickly the terminal speed is reached: $$v(r)=v_{\infty }{\left(1-R_{\star }/r\right)}^{\beta },$$(1)

where r is the distance to the center of the star and R_⋆ is the stellar radius. The larger β, the more gradual the acceleration.

In Vela X-1, observations of UV spectral lines led Gimenez-Garcia et al. (2016) to propose, for the wind of the donor star HD 77581 of spectral type B0.5Iae (van Kerkwijk et al. 1995), a β-law with v_∞~ 700 km s⁻¹ and β = 1. Later on, Sander et al. (2017b) computed a hydrodynamic steady-state structure for the wind stratification of the donor star in Vela X-1 based on a nonlocal thermodynamical equilibrium computation of the population on involved energy levels. Between the photosphere and the compact object and for our current purpose, the velocity profile predicted by Sander et al. (2017b) can be reasonably well fit by a β-law with v_∞~ 500 km s⁻¹ and β ~ 2.3. In our simulations, we consider two different values for β: 0.5, representative of a quickly accelerating wind, and 2, for a slowly accelerating wind. The terminal wind speed in Vela X-1 lies on the lower end of what is usually found in HMXBs (see Table 1 for a few other systems), which leads to a wind focused toward the accretor and susceptible to form transient wind-captured disks (El Mellah et al. 2019), a prediction recently corroborated by observations (Liao et al. 2020). Given how high the terminal wind speed usually is compared to the orbital speed in SgXBs, where the orbital period is typically of at least a few days, the wind ballistic streamlines are not significantly curved (El Mellah & Casse 2017). In several HMXBs, the gravitational influence of the compact object has little impact at the orbital scale and the wind should be essentially radial and isotropic around the donor star.

2.2 Clumps

In the two-zones model we work with, all the wind mass is contained in spherical clumps, each being uniformly dense, with an effectively void inter-clump environment. This approximation is justified given the observational findings that most of the mass of the wind material in the two prototypical HMXBs Cyg X-1 and Vela X-1 is concentrated in a small fraction of the volume (Sako et al. 1999; Rahoui et al. 2011). The clumps are carried with the wind, on purely radial orbits. They do not merge with each other, which enables us to use the continuity equation to deduce n, the number of clumps per unit volume at a certain distance r from the stellar center (not to be confused with the particle number density of the smooth wind): $$n=\frac{\stackrel{.}{\mathcal{N}}}{4\pi r^{2}v},$$(2)

where $\stackrel{.}{\mathcal{N}}$ is the rate at which clumps are emitted from the stellar photosphere, given by Ṁ ∕m. Ṁ is the stellar mass-loss rate and m is the mass of one clump, assumed independent of the distance to the star. All the clumps have the same mass m and the same radius R_cl,2 at two stellar radii, which are set as input parameters of the model. The evolution of the clump radius with the distance to the star can be set based on observational or physical insights concerning how clumpy the wind is. The latter is quantified through the clumping factor f or the volume filling factor f_vol: $$\{\begin{array}{l}f=\langle {\rho }^2\rangle /{\langle \rho \rangle }^2\hfill \\ f_{\text{vol}}=\rho /{\rho }_{\text{cl}},\hfill \end{array}$$(3)

where ρ is the smooth wind density and ρ_cl is the mass density of a clump. When the inter-clump environment is empty, ⟨ρ⟩ = ρ_clf_vol but $\langle {\rho }^2\rangle ={\rho }_{\text{cl}}^2{f}_{\text{vol}}$, such that $\langle {\rho }^2\rangle /{\langle \rho \rangle }^2=1/f_{\text{vol}}=f$.

In the case of a uniform clumping factor through the wind, the evolution of the clump radius R_cl with the distance r to the stellar center is best represented by: $$R_{\text{cl}}\propto {(v{r}^2)}^{1/3}.$$(4)

It is also theprofile which ensures that the mass density of constant mass clumps evolves as the density of a smooth wind (Ducci et al. 2009). Alternatively, it has been proposed based on the observed and simulated radial stratification of f in O-stars (Sundqvist & Owocki 2013) that in a simple model such as here the clumps might expand according to (see the Appendix A in Sundqvist et al. 2012b): $$R_{\text{cl}}\propto r.$$(5)

Although the exact clump radius expansion law is unknown, radiative hydrodynamics simulations indicate that the radius of the clumps tend to increase as they flow away from the star (Sundqvist et al. 2018). Hereafter, these two expansion laws are referred to as smoothly expanding clumps and as linearly expanding clumps respectively.

The absorption produced by clumps can be quantified through the porosity length, the effective clump column density divided by the mean density (Owocki & Cohen 2006). The porosity length h can be interpreted as the distance in a smooth wind along which the column density is the same as the column density of one clump averaged over all possible impact parameters p: $${\rho }_{\text{cl}}\left(\frac{1}{2R_{\text{cl}}}{\displaystyle {\int }_{-R_{\text{cl}}}^{R_{\text{cl}}}2}\sqrt{R_{\text{cl}}^2-p^2}\text{d}p\right)={\rho }_{\text{cl}}\frac{\pi }{2}{R}_{\text{cl}}=\langle \rho \rangle h,$$(6)

which means that the porosity length profile can be written as a function of the clumping factor f: $$h=\frac{\pi }{2}f{R}_{\text{cl}}.$$(7)

The profiles for the clump radius, the porosity length, and the clumping factors are shown in Fig. 1 for expansion laws in Eqs. (4) and (5) (in blue and red, respectively). Linearly expanding clumps lead to a clumping factor peaking at r = (β + 1)R_* while in a wind with smoothly expanding clumps, the clumping factor is uniform. Additionally included in the bottom panel in Fig. 1 are shaded regions which represent typical clumping factor profiles obtained from 1D line-driven instability simulationsof two O supergiants and two B supergiants (Driessen et al. 2019). Observations of isolated O stars yield clumping factors which qualitatively agree with the ones computed in these simulations (Puls et al. 2006; Najarro et al. 2011; Sundqvist & Owocki 2013). We point out that we are not aiming here to obtain absolute clumping factors but rather focus on the comparison with the formalism developed in this paper. Indeed, variations in clumping factor are to be expected between one-dimensional and multidimensional line-driven instability simulations (Sundqvist et al. 2018). The latter shows that clumping factors are reduced with respect to one-dimensional simulations because in 1D, the flow is unable to circumvent clumps formed in the simulation (really shells in one dimension). Accelerated interclump matter will flow into dense, clumpy material which artificially increases the clump density ρ_cl and, hence, the clumping factor. However, in Fig. 1, the comparison between the profiles computed in numerical simulations (shaded regions) and those computed based on the analytic description of the clump radius expansion models (solid lines) indicates that the two simple expansion laws described above are still realistic for the wind of O and B supergiants and that the corresponding clumping factor profiles follow reasonably well the trends derived by the simulations.

Simulations and observations further suggest that clumps masses (resp. radii) lie in-between a few 10¹⁷ and 10¹⁹ g (resp. 0.5 and 8% of the stellar radius), which, after considering realistic normalization in SgXBs, motivates our choice of dimensionless parameters (see Tables 2 and 3).

On purpose, we compute the variability of the total column density without mentioning the X-ray luminosity itself, the optical depth or the extinction coefficient since the model we develop does not include energy-dependent opacity. In follow up work, we plan to address the question of the ionization structure of the wind, with a proper treatment of the X-ray ionizing feedback from the X-ray source. The present work concerns only the total column density, that is to say the amount of material integrated along the LOS, not the absorbing column density deduced from fits of the continuum in X-ray spectra. Changes of the ionization parameter and thus, of the opacity, along the LOS could induce significant variations in the absorbing column density while leaving the total column density unchanged (see also discussion in Sect. 6.1.2).

Fig. 1 Top panel: clump radius (solid lines) and porosity length (dashed lines) in units of R* as a function of the distance to the stellar center. In blue (resp. in red) are represented the profiles for the clump radius expansion law (4) (resp. Eq. (5)). The dimensionless clump radius at r = 2, indicated with a black dot, was set to 0.005, and the dimensionless clump mass to 4 × 10−7 (see Sect. 2.4 for normalization units). Bottom panel: clumping factor profiles for both expansion laws. The red (resp. blue) shaded region represents the shape of the clumping factor profile obtained for O (resp. B) supergiants by Driessen et al. (2019). Spectroscopically, these OB supergiants can be considered representing typical early- and middle-type stars.

2.3 Orbiting X-ray source

In HMXBs, most of the X-rays we observe originate from close to the NS magnetic poles or from close to the innermost regions of the accretion disk surrounding the BH, depending on the nature of the accretor. At the orbital scale, these regions of a few tens to hundreds km are so small that the accretor can be represented as an orbiting X-ray point source in our model.

In classic SgXBs, the eccentricity of the orbit of the compact object is low (Walter et al. 2015) and so is the orbital separation, typically within ~1.5–3 stellar radii,motivating our choice of orbital separations. The present model could be extended to include any type of predefined orbit but for the sake of simplicity, we only treat circular orbits (see Sect. 6.1.3 for a discussion on eccentric orbits). The orbit is entirely determined by two parameters: the orbital separation a and the inclination angle i between the line of sight and the orbital angular momentum axis. In Fig. 2 we show two classic HMXBs, Vela X-1 and Cygnus X-1 as seen from Earth. The NS-hosting Vela X-1 is essentially edge-on as it is an eclipsing HMXB, that is it has a high inclination angle (Joss & Rappaport 1984). The BH-hosting Cygnus X-1 is much more face-on, with a low inclination angle of ~27° (Orosz et al. 2011).

2.4 Physical parameters of the model

Our model includes a very limited number of six degrees of freedom to evaluate the impact of the wind micro-structure on the time variability of the column density. The orbit is determined by the orbital separation a, the orbital inclination i and the orbital period P. The clump properties are the mass and the radius at 2R_*. The dynamics of the wind is set by the β-exponent. The terminal wind speed v_∞ is used as the normalization quantity for speeds and thus, does not affect the shape of the numerical solutions. The stellar radius R_* is the normalization of the lengths and with thestellar mass-loss rate Ṁ as an additional normalization unit, all our results can be scaled to physical units. The values for the six shape parameters we explore are given in Table 2 and correspond to 600 different configurations. Among these 600 configurations, a few correspond approximately to the values in the two classic HMXBs, Vela X-1 and Cygnus X-1, for which we provide, for illustrative purposes, the corresponding scales in Table 3. However, our aim is not to limit our study to these two systems specifically but rather to explore a set of realistic configurations which are thought to cover most of the known SgXBs, including Vela X-1 and Cygnus X-1. From now on, unless specified otherwise, the equations are given in their dimensionless form.

Fig. 2 Representation of the orbit (in green) of the compact object around its stellar companion (in blue) for Vela X-1 (top panel) and Cygnus X-1 (bottom panel) as seen from Earth, where the eccentricity of the orbit has been neglected. Only a fraction of the clumps susceptible to contribute to the column density along the orbit have been represented (in semi-opaquegrey). The black cross indicates the point of inferior conjunction on the orbit.

2.5 Numerics

2.5.1 Clumps

We work in the observer frame centered on the donor star and neglect stellar rotation and wobbling around the center of mass. The latter assumptions are expected to have little impact on the column density variability since the wind essentially flows away from the star, with terminal wind speeds typically larger than the orbital speed by a factor of a few. The position of the clumps is advanced from a lower spherical boundary of radius r_m = 1.2 to an upper spherical boundary of radius r_M. Between r = 1 and r = r_m, the wind is assumed to be smooth, to represent the finite distance from the photosphere beyond which significant clumping is expected (see Fig. 1 for theoretical estimates, and for example Puls et al. 2006, for observational constraints). Clump emission is isotropic in the sense that clumps are injected at the inner boundary r = r_m with an equal probability per solid angle. At each time step, clumps beyond r = r_M are deleted while clumps are created at the basis at a rate $\stackrel{.}{\mathcal{N}}=1/m_{cl}$. We can estimate a priori the number $\mathcal{N}$ of clumps in the simulation space by writing: $$\mathcal{N}={\displaystyle {\int }_{r_m}^{r_M}n}4\pi r^2\text{d}r=\frac{1}{m_{\text{cl}}}{\displaystyle {\int }_{r_m}^{r_M}\frac{\text{d}r}{v}}.$$(8)

A preliminary integration phase of duration $\mathcal{N}{m}_{\text{cl}}$ guarantees that when we start computing column densities, the whole simulation space is populated with clumps, from the lower to the upper radial boundaries. The position of each clump is advanced using a 4th order Runge-Kutta integration method, with a time step Δt equal to P∕(3 × 10⁴). This time step can be compared to the minimum of the duration for a clump to pass in front of a static point source: $$\Delta t=\text{min}\left[2R_{\text{cl}}(r)/v(r)\right],$$(9)

which is found, for both expansion laws, at r = β + 1. Even in the least favorable configuration, this characteristic time scale is still hundreds of times larger than Δt, which ensures that we resolve the impact of all clumps on the column density, even the smallest ones with grazing impact parameters with respect to the LOS. For realistic parameters and r_M = 30a (see next section), the total number of clumps can reach up to 10⁹ which need to be integrated over 3 × 10⁵ of time steps in order to get the evolution of the column density over 10 orbital periods, our full integration time. In order to reduce the number of clumps and speed up the computation, we retain only the clumps susceptible to intercept the LOS and discard the others. As visible in Fig. 2, it means that the number of clumps to consider is much smaller, in particular when the orbit is close to edge-on.

2.5.2 Maximum distance of integration

Every time step, we cast a ray from the X-ray point source to the observer and sum up the contribution of each segment which intercepts a clump to obtain the total column density at this time step. To determine at which distance r_M we should stop the integration of the column density along the LOS, we quantified the contribution of each logarithmic radial bin between the compact object and large distances. We worked at inferior conjunction (ϕ = 0.5) such as these results are a fortiori valid when the X-ray source is buried deeper into the wind. In Fig. 3 is represented the cumulated column density in blue, and the contribution of each bin in red. As we move away from the accretor (located at r = 2 in this example), the column density per logarithmic radial bin decreases quickly: 10 orbital separations away from the stellar center, the column density integrated over a certain interval is ~50 times lower than the one integrated over a range 10 times shorter but around r = a. This common trend for all the sets of parameters we considered drove us into setting the outer edge of the simulation domain at r_M = 30a, a distance at which the column density has essentially converged.

3 Analytical predictions

3.1 Smooth wind scenario

In the case of a smooth (i.e., homogeneous, without clumps) and isotropic wind with a radial velocity profile which follows a β-law, the normalized column density profile along the orbit of the accretor can be written as: $$N_H(\varphi )=\frac{1}{4\pi }{\displaystyle {\int }_{z(\varphi )}^{\infty }\frac{\text{d}z}{r^2{\left(1-1/r\right)}^{\beta }}},$$(10)

where r (function of z and of the orbital elements) and z are respectively the distance to the stellar center and the coordinate along the LOS joining the X-ray point source to the observer, both normalized with the stellar radius R_*, and ϕ is the orbital phase (from 0 to 1). This profile depends on three shape parameters and one scale. The three shape parameters are the exponent β of the velocity profile and the orbital elements. For a circular orbit, the latter are the inclination angle i of the orbital angular momentum with respect to the LOS and a, the orbital separation divided by the stellar radius. They both play a role to determine the inferior bound of the integral, z(ϕ), which locates the compact object along the LOS, and to convert r into z. The integral can then be performed numerically to compute any column density profile in a smooth and isotropic wind. In eclipsing HMXBs, the median N_H profile at egress is the most constraining due to its sharp decrease. The fit of the measured column density profile with Eq. (10) also grants access to the mass-loss rate via the column density scale N_H,0, provided we can estimate the terminal wind speed and the stellar radius (see Table 3 for scalings).

To illustrate how a, i, β and $N_{H_0}$ can be obtained based on this fitting procedure, we focus on an actual HMXB. The classic SgXB 2S0114+650 (also known as LSI+65 010) is an excellent candidate to apply the model presented in this paper, thanks to its grazing orbit. Sanjurjo-Ferrrín et al. (2017) provide a ratio a∕R_* of 1.34–1.65. Given the high mass ratio (Reig et al. 1996), the absence of accretion disk (Koenigsberger et al. 2003) and the moderate X-ray luminosity (1.4 × 10³⁶ erg s⁻¹, Sanjurjo-Ferrrín et al. 2017), we can safely assume that the star does not fill its Roche lobe. It brings us to the conclusion that the upper limit of this range is to be preferred for the ratio a∕R_*. Pradhan et al. (2015) showed that it was not eclipsing by measuring a finite column density at superior conjunction. However, the sharp decrease of the column density following superior conjunction indicates that the LOS probes deep layers of the wind, near the photosphere. Consequently, the orbital inclination has to be between 40° and 50° approximately, which is also consistent with results obtained from radial velocity measures of the donor star (Crampton et al. 1985; Farrell et al. 2008; Koenigsberger et al. 2006). We fit the N_H data between orbital phases 0 and 0.1 reported in Pradhan et al. (2015) with the smooth wind profile given by Eq. (10). We could find a reasonable match for the following parameters: a = 1.6, i = 50°, β = 1 and N_H,0 = 1.2 × 10²³ cm⁻². Reig et al. (1996) estimated the terminal wind speed to be 1200 km s⁻¹, such that if the stellar radius is R_* = 25R_⊙, we deduce from N_H,0 a stellar mass-loss rate Ṁ ~ 7 × 10⁻⁷ M_⊙ yr⁻¹, very similar to the inferred mass-loss rate of the donor star in Vela X-1 (Gimenez-Garcia et al. 2016). Spectral characterization in the optical and UV waveranges of the stellar atmosphere with non local thermodynamical equilibrium analysis can bring additional constrains to break the degeneracy and disentangle between the different parameters which come into play in the value of N_H,0 (Hainich et al. 2020).

Due to the conservation of mass, the averaged column density should remain the same whether the mass is gathered in clumps or spread in a smooth wind: the properties of a smooth wind are expected to match the time-averaged properties of a clumpy wind where clumps are realistically small and numerous, a conjecture which will be validated in Sect. 4.2. The micro-structure of the wind will only manifest through the time variability of the column density.

Fig. 3 Column density, per logarithmic radial bin (red, logarithmic axis on the left) and cumulative (blue, linear axis on the right), as a function of the distance to the stellar center. The position of the compact object is indicated with a dashed black line. The total column density up to 30 orbital separations is ~97% of the column density computed for a smooth wind integrated up to thousands of stellar radii. The dimensionless parameters are mcl = 10−5, Rcl = 0.08, a = 2, P = 80 and i = 45°.

3.2 Standard deviation

Following the arguments by Owocki & Sundqvist (2018), we may write the variance of the column density between a source located at z₀ along a LOS and an observer at z = ∞ as, in physical units and for a mean molecular weight of one proton mass m_p : $$\delta N_H^2={\displaystyle {\int }_{z(\varphi )}^{\infty }{\left(\frac{\rho }{m_p}\right)}^2}h\text{d}z,$$(11)

where h is again the porosity length and ρ = ⟨ρ⟩ is the mean mass density of the isotropic wind. Using the continuity equation in combination with Eq. (6), we obtain the dimensionless standard deviation δN_H in column density around the value set by the smooth and isotropic wind as a function of the orbital phase: $$\delta N_H=\frac{3}{32\pi }{\left(m_{\text{cl}}{\displaystyle {\int }_{z(\varphi )}^{\infty }\frac{\text{d}z}{R_{\text{cl}}^2{r}^{2}v}}\right)}^{1/2}.$$(12)

The expansion law for the clumps can then be reinjected in the formula above. For example, for linearly expanding clumps, we then have, if the velocity profile is given by the β-law (1): $$\delta N_H=\frac{3}{8\pi }{\left(\frac{m_{\text{cl}}}{R_{\text{cl},\text{2}}^2}{\displaystyle {\int }_{z(\varphi )}^{\infty }\frac{\text{d}z}{r^4{\left(1-1/r\right)}^{\beta }}}\right)}^{1/2}.$$(13)

This shows explicitly how, for a given orbit and wind acceleration, the standard deviation of the column density depends on the ratio $\sqrt{m_{\text{cl}}}/R_{\text{cl},\text{2}}$ for any power-law clump radius expansion.

We performed simulations with both clump radius expansion laws for the 300 simulations with high-mass clumps (m₀ = 10⁻⁵), less computationally-demanding due to the lower total number of clumps in the simulation space. The results were similar enough that all the conclusions we draw in the next sections are equally valid for both laws. It was to be expected because in addition to the expression (12) above, Fig. 1 indicates that clump radii and porosity lengths within 3R_*, in the region where the clumps contribute the most to the final column density, are alike. Hereafter, results are presented for the smoothly expanding law (Eq. (4)).

4 Results

In this section, we focus on the results of our simulation, in particular the time-evolution of column density (Sect. 4.1), peak to peak variability (Sect. 4.2), and coherence time scales (Sect. 4.3). We compare to theoretical predictions and relate the so derived properties of absorption variability to the underlying properties of the clumpy winds. The observational implications and in particular the challenges and prospects are discussed later, in Sect. 6.

4.1 Time evolution of the column density

The main output of each simulation is a time series of the column density. This data is represented for three simulations in Fig. 4 over a few orbital periods. As the compact object orbits the donor star, the amount of material integrated along the LOS changes because of two effects: the motion of the X-ray point source around the star and the clumps crossing the LOS. The latter induces a variability on short time scales whose statistical properties relate to the micro-structure of the wind and are discussed in more detail in Sects. 4.2 and 4.3. The orbital motion, on the other hand, causes a smooth periodic modulation of the column density, more contrasted for higher inclination and orbital separation, and which vanishes for face-on configurations where the distance between the accretor and the observer remains constant (see bottom panel in Fig. 4). At high inclination, eclipses can occur during which the column density can not be measured (see gaps in the top panel in Fig. 4, for a 90° inclination). The duration d of an eclipse with respect to the orbital period P relates to the ratio a∕R_* via: $$\frac{d}{P}=\frac{1}{\pi }\arcsin \left(\frac{R_{*}\sin i}{a}\right)~\frac{1}{\pi }\frac{R_{*}\sin i}{a},$$(14)

where the approximation is valid within 10% for any inclination provided a > 1.6 R_*, a condition usually met in HMXBs (e.g., Falanga et al. 2015). In Fig. 4, we represented the smooth wind solution obtained with Eq. (10) with a dashed green line. The measured column densities oscillate around this profile. Because we only cover the case of circular orbits in this paper, the smooth wind as a function of the orbital phase ϕ is symmetric with respect to the axis ϕ = 0 and ϕ = 0.5 (superior and inferior conjunction respectively). We highlighted the short term variability of the column density with inserts showing zoomed-in versions of the time series. Each insert contains approximately 200 data points which are visible and indicate that the impact of the clumps is resolved up to very short time scales. The total column density consists of cumulative contribution of many clumps with different radii and densities, which indicates that accurate knowledge on the clumpiness of the stellar wind can be gained from the time variability of the column density.

4.2 Column density peak-to-peak variability

4.2.1 Variability amplitude as a function of the orbital phase

In Fig. 4, the spread of the column density around its smooth wind value changes from a simulation to another. In this section,we determine the main parameters which set the measured standard deviation but beforehand, we represent the computed column density as a function of the orbital phase in Fig. 5 for a fiducial set of parameters. We computed the column density during each simulation at 50 locations equally spaced on the orbit simultaneously every 100 time steps such as each orbital phase bin in Fig. 5 contains 6000 values of N_H (after folding the orbit from ϕ = 0 to ϕ = 0.5). For each bin, we compute the median (black steps) and the standard deviation. The black shaded region around the median shows the values between the 20th and the 80th percentiles while the grey shaded region is delimited by the 5th and the 95th percentiles. Simultaneous computation of N_H along the orbit enables us to get a much better statistics to compute the standard deviation at each orbital phase than using the time series presented in Sect. 4.1. It does not introduce biases due to time coherence of the signal associated to clumps if the simulation duration (ten orbital periods) is high compared to the coherence time scale, an assumption we validate a posteriori (see Sect. 4.3). These results indicate that except for very large (R_cl = 0.08) and sparse clumps (m_cl = 10⁻⁵), the median N_H profile is representative of the smooth wind profile with the same parameters but where clumps are replaced with a smooth wind of same mass-loss rate. It means that when a column density is measured from observations around a certain orbital phase (e.g., during a 10 ks observation, which would correspond to an orbital phase width on the order of 1∕50 in Cygnus X-1), it can be compared to Fig. 5 for a set of parameters believed to match the observed system. If the column density at this very same orbital phase is repeatedly measured over different orbits, the histogram of values should reproduce a vertical slice in Fig. 5, with the median corresponding to the value for a smooth wind of same mass-loss rate.

Fig. 4 Three different evolutions of the column density as a function of time over several orbits for i = 90° (top panel), i = 45° (middle panel) and i = 0° (bottom panel) – and different clump radii, orbital separations, periods and β exponent.The dashed green line is the corresponding smooth wind column density profile. The top profile displays eclipses due to the high inclination. The insets show zoomed in portions of the curves where the numerical time sampling of the column density, 3 × 104 times per orbital period, is visible. The origin of time is set at inferior conjunction (i.e., t = 0 at ϕ = 0.5).

Fig. 5 Column density from superior conjunction (ϕ = 0) to inferior conjunction (ϕ = 0.5) for an orbital inclination i = 45°. In black is represented the median over 10 orbits and the red dashed line is the smooth wind solution from Eq. (10). The dark (resp. light) grey shaded region is bounded by the 20th to 80th percentiles (resp. the 5th to 95th percentiles). The parameters are mcl = 10−5, Rcl = 0.01, β = 0.5 and a = 2.

4.2.2 Comparison to analytical predictions

For each simulation, we compared the profiles of the standard deviation as a function of the orbital phase to the predicted profiles derived in Sect. 3.2 (see Fig. 6). Apart from near superior conjunction in a few configurations, they match within a few percent. The discrepancies sometimes observed near superior conjunctions can be fully attributed to two effects. First, when the system is eclipsing, the orbital bin which contains ingress or egress underestimates the variability because it averages the standard deviation over a finite orbital phase interval which partly includes the eclipse, during which there is no variability. Second and more importantly, when the LOS intercepts the region r ∈ [1;1.2] around the donor star where we assumed that the wind is smooth, the standard deviation is also lowered with respect to the prediction of the analytical model where the wind is assumed to be clumpy all the way down to the stellar photosphere.

Observational hints in favor of clumpiness near the star exist (Puls et al. 2006; Cohen et al. 2011; Torrejón et al. 2015) but our model provides a way to quantify this: if the standard deviation at superior conjunction measured in HMXBs where the orbit is grazing is lower than what predicted from fitting the standard deviation profile far from superior conjunction with the formula (11), then it would support the scenario of a smooth wind base. Monitoring the δN_H profile in highly inclined systems might thus enable us to constrain the onset of the clump-forming region. For instance, the bottom panel in Fig. 6 corresponds to a configuration where, although the system is not eclipsing (since at superior conjunction, when ϕ = 0, the projected distance to the star is a sin i > 1), the LOS to the X-ray source probes the deep smooth layers near the photosphere (since at superior conjunction, a sin i < 1.2). We see thatthe standard deviation is lower than predicted by the analytic model (which assumes that the wind is clumped all the way down to the stellar surface) up to an orbital phase of ~0.07, when the LOS no longer intersects the smooth wind region. In SgXBs where the orbital inclination and separation are such that the LOS probes the smooth wind region near the star, the location of the maximum in the δN_H profile is thus a direct measure of the distance to the stellar photosphere where clumping occurs.

Fig. 6 Standard deviation of the column density as a function of the orbital phase computed numerically (in blue) and analytically (black dots) for simulations with mcl = 10−5, Rcl = 0.02. Top panel; β = 2, a = 2 and i = 90° and bottom panel: β = 0.5, a = 1.6 and i = 45°. The discrepancy at low orbital phase is partly indicative of the region near the stellar photosphere where the wind in the numerical model is assumed to not be clumpy yet.

Fig. 7 From top to bottom row: orbital inclination is 65°, 25° and 0°. Left column: standard deviation of the column density at inferior (resp. superior) conjunction in red (resp.blue) as a function of the ratio $\sqrt{m_{\text{cl}}}/R_{\text{cl},\text{2}}$. Each couple (a, β) leads to a different proportionality constant between δNH and $\sqrt{m_{\text{cl}}}/R_{\text{cl},\text{2}}$, hence the large dispersion. The dashed line in the top left panel indicates the y = x line. Right column: same as before but Rcl,2 has been replaced by the clump radius at the point along the LOS between the X-ray source and the observer the closest from the donor star. The dispersion is reduced, especially at inferior conjunction.

4.2.3 Dependence on parameters

Finally, we compare the standard deviation of the column density at inferior and superior conjunctions, δN_H,inf and δN_H,sup, for different parameters. The dependency of δN_H on $\sqrt{m_{\text{cl}}}/R_{\text{cl}}$ in Eq. (11) invited us torepresent δN_H,inf and δN_H,sup as a function of $\sqrt{m_{\text{cl}}}/R_{\text{cl},\text{2}}$ in the left column in Fig. 7 for an orbital inclination of 65°, 25° and 0° (top, middle and bottom rows). For eclipsing HMXBs, the standard deviation at superior conjunction is not defined and thus not represented. We vary the orbital separation (different markers) and differentiate quickly and slowly accelerating winds (filled and empty markers respectively). In blue (resp. red) are shown the values of δN_H,sup (resp. δN_H,inf). These results do not depend on the orbital period, which was expected since the orbital period is so much longer than the time scale at which N_H varies. All other parameters being fixed, we see that the relation $\delta N_H\propto \sqrt{m_{\text{cl}}}/R_{\text{cl},\text{2}}$ whose slope is indicated in the first panel with a black dashed line is verified for each set of ten points (five different values of R_cl,2 and two values of m_cl). The ratio between the standard deviation at superior and inferior conjunctions is hardly sensitive to the orbital separation. As the orbital separation increases, the difference between slowly and quickly accelerating winds vanish at inferior conjunction: while the wind speed at r = 1.6 is more than four times higher for β = 0.5 than for β = 2, this ratio drops to less than 2 at r = 3. The difference between β = 0.5 and β = 2 remains significant at superior conjunction though, even for r = 3, since the LOS probes regions r < a enclosed within the orbit. A higher β leads to a higher standard deviation due to the higher number of clumps per unit volume at a given distance from the star (see Eq. (2)). The standard deviation at superior conjunction is always higher than at inferior conjunction, except for face-on orbits where they coincide: at superior conjunction, the LOS intercepts a region closer from the star than at inferior conjunction, where the clumps are small and packed, which increases the standard deviation.

Given the dominating contribution of the densest clumps pinpointed in Fig. 3, we changed the clump radius used in the x-axis of the plots in Fig. 7 so as it corresponds to the smallest clump radius encountered along the LOS, R_cl,min, that is the clump radius at the point along the LOS the closest from the donor star (right column in Fig. 7). At inferior conjunction, this point is simply the position of the compact object but at superior conjunction, it is a point located deep in the wind if the inclination is significantly high, such as the relevant clump radius becomes, for smoothly expanding clumps: $$R_{\text{cl},\text{min}}(\varphi =0)=R_{\text{cl},\text{2}}{\left(\frac{a\cos i}{2}\right)}^{2/3}{\left(2-\frac{2}{a\cos i}\right)}^{\beta /3}.$$(15)

This correction enables us to get a much lower dispersion of the points along the $\delta N_H\propto \sqrt{m_{\text{cl}}}/R_{\text{cl},\text{min}}$ axis and proves that the standard deviation is an excellent tracer of the ratio $\sqrt{m_{\text{cl}}}/R_{\text{cl},\text{min}}$. Once the parameters a, i, β and N_H,0 have been constrained thanks to the median column density profile as a function of the orbital phase, a precise measure of the standard deviation grants us access to the ratio $\sqrt{m_{\text{cl}}}/R_{\text{cl},\text{min}}$. The accuracy of this obtained ratio is maximal if δN_H is measured at inferior conjunction, where the dispersion is the lowest. The latter effect can be explained by the fact that at inferior conjunction, the LOS is best aligned with the radial direction from the donor star (with exact alignment for an edge-on system) so the column density is mostly localized near the accretor and the contributions quickly drop as we go toward the observer. Another advantage of working at inferior conjunction is to obtain an accurate median profile $N_H\left(\varphi \right)$ with fewer observations. The spread which subsists in the right column in Fig. 7 is due to the minor but non negligible contribution of clumps with radii R_cl > R_cl,min on the standard deviation of the column density.

These results illustrate how the standard deviation on column density can be exploited to constrain the ratio $\sqrt{m_{\text{cl}}}/R_{\text{cl},\text{min}}$, or alternatively the porosity length since $\sqrt{m_{\text{cl}}}/R_{\text{cl}}\propto \sqrt{h}$. In order to disentangle the clump mass and radius, it is necessary to analyze another output of our simulations, the coherence time scale.

4.3 Coherence time scale

4.3.1 Autocorrelation function

In the previous section, we extracted information on the characteristic standard deviation of the column density at each orbital phase for each set of parameter, but we did not exploit the information on the time evolution of the column density. The time series contain additional information, the coherence time scale, which is the characteristic duration over which the changes in column density are smooth and, on average, predictable. If observations of HMXBs are carried out over intervals lasting a significant fraction of the orbital period, the column density as a function of time can, in theory, be extracted on time scales shorter than this coherence time scale (see Sect. 6.2 for example estimate of necessary time resolution).

In order to determine the coherence time scales in each simulation, we first compute the relative difference between the time series presented in Sect. 4.1 and the smooth wind solution (10), and subtract the median of the obtained signal afterwhile. Since we expect the coherence time scale to vary with the orbital phase, we must treat separately the different orbital phases to avoid smearing out the coherence time scale. We cut the time series in twenty equal orbital phase intervals. In spite of the motion of the X-ray source along its orbit, we can still use the symmetry with respect to ϕ = 0.5 and gather the data at ϕ = 0.25 and ϕ = 0.75 for instance, since both are expected to have the same coherence time scale. For each simulation, we are left with ten orbital phase intervals of width 0.05, between ϕ = 0 and ϕ = 0.5, each containing twenty time series of the column density of duration P∕20. The autocorrelation function of each of this individual time series is then computed and the twenty autocorrelation functions of each orbital phase bin are averaged to obtain a representative autocorrelation function for each of the ten intervals¹. With this method, we can capture coherence time scales ranging from P∕(3 × 10⁴) (the time resolution) up to P∕40. As an illustration, the ten autocorrelation functions of a fiducial simulation are represented in Fig. 8, from inferior conjunction (dark purple) to superior conjunction (yellow).

The autocorrelation function is maximal in zero since the time series is fully correlated to its zero time shifted counterpart. As the shift increases, the correlation decreases until the time shift is so large that the two functions are no longer correlated at all. The coherence time scale corresponds to the time beyond which the signal changed enough to consider that it lost the memory of its previous values. We measure it by taking the time shift when the autocorrelation has been halved, τ. We evaluate the uncertainty on the obtainedvalue of τ by looking at the face-on configurations, where τ should not depend on the orbital phase: it remains indeed constant, and varies by at much 5%.

The coherence time scale measured at inferior conjunction is typically higher than the one at superior conjunction. The LOS probes deeper wind layers at superior conjunction where clumps are smaller but also slower, such that a shorter coherence time scale is not obvious without accounting for the orbital speed, which turns out to be larger than the speed of the clumps near the stellar photosphere.

4.3.2 Dependence on parameters

In this section, we focus on the coherence time scales computed in bins centered on phases ϕ = 0, ϕ = 0.25 and ϕ = 0.5 and of width 0.1 each. To study the evolution of the coherence time scale with the parameters, we want to compare the coherence time scale τ to the amount of time a single clump intercepts the LOS to the accreting compact object that is the duration for the X-ray source to cross a clump projected along the LOS. We call it the self-crossing time and write it τ_SC. Due to the short orbital periods and the high total masses in HMXBs, the orbital speed can be as high as a few hundreds of kilometers per second. This speed is not negligible compared to the speed of the wind, especially if the acceleration is very gradual (i.e., high β). We compute the local relative speed Δv between the orbiting X-ray point source and the densest clumps along the LOS: $$\begin{array}{ll}\Delta v=\hfill & [v_{\text{orb}}^2\left(1-{\sin }^2\Phi {\sin }^{2}i\right)+v_{\text{cl}}^2\left(1-{\cos }^2\Phi {\sin }^{2}i\right)\hfill \\ \hfill & \mathrm{...}-2v_{\text{orb}}{v}_{\text{cl}}\cos \Phi \sin \Phi {\sin }^{2}i{]}^{1/2}\hfill \end{array}$$(16)

where Φ = 2πϕ, v_orb = 2πa∕P and v_cl = v(r = a). From this, we deduce the self-crossing time: $${\tau }_{\text{SC}}=2R_{\text{cl}}/\Delta v.$$(17)

To compute the clump radius R_cl which enters this equation, we focus on the clumps which contribute the most to N_H, like we did in Sect. 4.2.3: at inferior conjunction and ϕ = 0.25, we compute the clump radius at r = a while at superior conjunction, we compute the clump radius at the point along the LOS which is the closest from the donor star.

Figure 9 provides an overview of the ratios τ∕τ_SC for all simulations. Except for a handful of simulations where clumps are the largest and the scarcest, the coherence time scale matches the self-crossing time within 50%. This result means that a fairly accurate estimate of the clump self-crossing time and, by then, of the characteristic dimension of the clumps, can be obtained by measuring the full width at half maximum of the autocorrelation function of a column density time serie. Within the ranges of parameters we explored, τ_SC depends on all the 6 dimensionless parameters summarized in Table 2 except on the clump mass. It provides a convenient way to distinguish the size and the mass of the clumps, which were intertwined when we analyzed the standard deviation oncolumn density in Sect. 4.2.

Fig. 8 Autocorrelation function of the column density as a function of the time shift. From purple to green and yellow, we represent the autocorrelation functions of the data at different orbital phases, from inferior conjunction (dark purple) to superior conjunction (yellow). Simulation with mcl = 4 × 10−7, Rcl = 0.04, β = 2, a = 3, i = 65° and P = 85.

Fig. 9 Stacked histograms of the ratios of the measured coherence time scale τ compared to the clump self-crossing time τSC for all 600 simulations at inferior conjunction (ϕ = 0.5, purple), ϕ = 0.25 (green) and superior conjunction (ϕ = 0, yellow).

5 Observational diagnostics

5.1 Procedure

The procedure we propose to determine clump properties from the variability of the column density follows the 3 successive steps below.

First, we need to have constrained values for the orbital inclination i and separation a, the shape of the velocity profile (i.e., the β-exponent) and the column density scale $N_{H_0}$. This can be done based on the median column density profile and if possible, additional insights (e.g., radial velocity profiles or the extent of the stellar Roche lobe).

Then we can make use of the N_H time series provided by observations, for instance through time-resolved X-ray spectroscopy (e.g., Martínez-Núñez et al. 2014) or color–color diagrams (Nowak et al. 2011; Grinberg et al. 2020). The clump radius manifests as a coherence time scale imprinted in the signal which can be extracted from the width of the autocorrelation function at half its maximum. It typically ranges from a few 10⁻⁴ to a few 10⁻² times the orbital periods, which means that a time resolution of ~100 s is enough to determine the clump radius if a coherence time scale is detected, and to provide constraining upper limits on the size of the clumps otherwise. The relative speed projected on the plane of the sky of the clumps with respect to the X-ray point source is lower near inferior conjunction. Consequently, for systems where the inclination is low enough that we can afford it (i < 65°), measures at inferior conjunction should be privileged since they yield longer coherence time scales that are easier to observationally constrain. Otherwise, measures at ϕ ~ 0.2−0.3 are to be preferred.

Finally, we can use the standard deviation of the column density, δN_H, to estimate the clump mass based on Eq. (12). At inferior conjunction, δN_H∕N_H,0 provides a direct measure of the ratio $\sqrt{m_{\text{cl}}}/R_{\text{cl}}$ whose accuracy is within a factor of 2 in the worst case, when no other parameter (a, i, P and β) is known (see Fig. 7). Since the clump radius R_cl entering this ratio is the same as the one deduced from the coherence time scale, the characteristic mass of a single clump can be derived.

The diagram in Fig. 10 summarizes the main ideas underlying this procedure for observational diagnostics of the wind micro-structure, showing the relationship between the mass and radius of an individual clump, and the standard deviation and coherence time scale of the column density δN_H, here calculated for a specific case. For illustrative purposes and to provide masses in physical units, we considered arbitrary scaling units and a simulation with intermediate parameters. The coherence time scales τ are given inorbital periods and directly map for clump radii. The additional information on δN_H enables to constrain the range of possible clump masses using: $$\begin{array}{ll}{m}_{\text{cl}}~\hfill & 3\times {10}^{22}\text{g}{\left(\frac{R_{\text{cl},m}}{R_{*}}\right)}^2\hfill \\ \hfill & \times \mathrm{...}{\left(\frac{R_{*}}{20R_{\odot }}\right)}^3\left(\frac{v_{\infty }}{1000\text{km}\times {\text{s}}^{-1}}\right)\hfill \\ \hfill & \times \mathrm{...}{\left(\frac{\dot{M}}{{10}^{-6}{M}_{\odot }\times {\text{yr}}^{-1}}\right)}^{-1}{\left(\frac{\delta N_H}{5\times {10}^{21}{\text{cm}}^{-2}}\right)}^2.\hfill \end{array}$$(18)

The corresponding porosity length is, for the parameters and mass unit used in Fig. 10: $$h~2\%R_{*}\left(\frac{m_{\text{cl}}}{4\times {10}^{17}\text{g}}\right){\left(\frac{R_{\text{cl},\text{2}}}{0.01R_{*}}\right)}^{-2}.$$(19)

This small value of the porosity length is typical of an environment with many low mass clumps susceptible to intercept the LOS and capture photons.

Fig. 10 Synthetic diagram summarizing the dependence of the coherence time scale τ and of the standard deviation on the column density δNH on the size of the clumps and their mass, for intermediate parameters (a = 2, i = 45°, P = 24 and β = 2). We used as normalization units m0 = 1024g for masses and NH,0 = 3 × 1023cm−2 for column densities.

5.2 Comparison to previous results for clumps in HMXBs

For clump sizes and masses chosen based on radiative hydrodynamics simulations (R_cl ~ 0.01 R_* and m_cl~ a few 10¹⁷ g, Sundqvist et al. 2018), the dispersion of the N_H values around the smooth wind profile is limited. The clumps with masses m_cl~ 10²²⁻²³ g and radii R_cl ~ a few R_*∕10 deduced from hard X-ray flares in SFXTs (Walter & Zurita-Heras 2007), which trace the intrinsic emission and thus, the mass accretion rate, are more likely not representative of the capture of a single clump by the compact object. Instead, the flares could be due either to large scale structures crossed by the accretor (e.g., corotating interaction regions, Bozzo et al. 2017b) or either to instabilities and gating mechanisms at the outer rim of the NS magnetosphere such as the propeller effect (Bozzo et al. 2008; Parfrey & Tchekhovskoy 2017), or within the accretion disk surrounding the black hole. These instabilities could be triggered bya clump, especially in SFXTs (Ferrigno et al. 2020), but involve much larger amounts of mass. For the edge-on HMXB VelaX-1, Grinberg et al. (2017) carried out a time-resolved analysis of the X-ray spectra and identified two coherent episodes of enhanced absorption of approximately an hour each, which would correspond to clumps of size R_cl ~ 5% R_*. However, the amplitude of the N_H increase, of several 10²² cm⁻², is too high to be explained by clump-induced stochastic variations. In edge-on systems, these observationscould be due to variability in the accretion wake or, when observations are performed at ϕ ~ 0.25 like in Grinberg et al. (2017), to the formation of transient high density axisymmetric structures around the NS magnetosphere such as an accretion disk (El Mellah et al. 2019; Liao et al. 2020).

6 Limitations of the model

6.1 Model-intrinsic limitations

6.1.1 Caveats for edge-on systems

Since the first simulations of wind accretion by Blondin et al. (1990), the flow in HMXBs is believed to form, in the orbital plane, a hydrodynamics bow shock which trails the accreting compact object due to the orbital motion. Manousakis & Walter (2015) showed that this dense structure was highly variable due to the X-ray ionizing feedback from the accretor onto the feeding material upstream, which alters the line acceleration as first described by Hatchett & McCray (1977) (see Sect. 6.1.2). We intentionally left large scale density modulations out of the model introduced in the paper so as to isolate the effect of the clumps themselves on the variability of the column density. In order to apply the present model in HMXBs with high orbital inclination (eclipsing and/or i > 65°), observations from superior conjunction (ϕ = 0), or egress if eclipsing, up to ϕ ~ 0.2−0.3 are relevant to determine the wind micro-structure. In this interval, the accretion wake points away from the observer and does not intercept the LOS. Elsewhere, the column density is significantly impacted by the contribution of structures not associated to the clumps. Instead, for intermediate inclinations (25° < i < 65°), data up to inferior conjunction can safely be used, especially if the wind speed is large compared to the orbital speed, leading to a low opening angle of the bow shock around the accretor. If the inclination is low (i < 25°), X-ray data at any orbital phase can be used but the weak modulation of photometric and spectroscopic signals with the orbital motion makes quantities such as the orbital period difficult to narrow down. Notice that a low inclination is much less likely than an edge-on configuration since the probability to observe a system with orbital inclination i goes as sin i.

6.1.2 X-ray photoionization

Wind accretion being less efficient at transferring stellar material to the compact object than Roche lobe overflow, the mass accretion rate in wind-fed HMXBs is generally moderate, leading to typical X-ray luminosities on the order of 10³⁶ erg × s⁻¹. Nevertheless, the influence of the X-rays emitted from the immediate vicinity of the compact object on the ionization state of the upstream wind can often not be neglected (see e.g., Boroson et al. 2003; Miškovičová et al. 2016; Grinberg et al. 2020, for observational results). In addition of altering the wind dynamics (e.g., Čechura et al. 2015; Sander et al. 2017b), the X-rays ionize the material surrounding the compact object and significantly decreases its opacity (Hatchett & McCray 1977; Ho & Arons 1987). In particular in the soft X-ray range used for absorption measurements in HMXBs, the effects of changing ionization lead to complex, nonlinear changes in observed spectral shape. In general, segments of the LOS near the accretor will pass through material that is less absorbent; the apparent column density inferred from X-ray spectroscopy is then, if a neutral absorber is assumed, smaller than the actual total column density of material and thus also than the column density computed in this paper.

A first correction can be made to the smooth wind profile (Eq. (10)) by accounting only for the material outside a region where the ionization parameter² ξ is lower than a critical value ξ_c. However, the exact value of ξ_c will depend on details of the observations used to obtain the absorption measurements from X-ray spectra, for instance the energy ranges considered and underlying spectral degeneracies, and on detailed (micro-) structure and composition of the wind. Only a proper radiative transfer computation can provide the insights we need on the wind ionization structure to correct the total column density N_H (Watanabe et al. 2006; Mauche et al. 2008). However, such simulations are not yet available for complex, dynamic clumpy winds. Interestingly, in SFXTs, the lower X-ray luminosity than in classic SgXBs should limit the impact of photoionization on the column density (Pradhan et al. 2018). Also, the SgXB 2S0114+650 discussed in Sect. 3.1 may, due to its moderate X-ray luminosity, be a good candidate to constrain the total column density. Investigations of the stochastic properties of the column density along the grazing orbit would be a precious element to gain knowledge on the size, mass and onset region of the clumps in the wind of the B1 supergiant donor star in 2S0114+650.

An additional problem may be posed by inhomogeneous and variable X-ray irradiation, such as when considering time scales smaller than the pulsar rotation period in HMXBs with accreting pulsars as compact object. In such systems, changes in X-ray flux can reach a factor of a few within 10–100 s, resulting in a dynamic illumination of the wind material. Around active galactic nuclei, the resulting effects of the X-rays on the clumps have been studied thanks to radiative transfer computations by Waters et al. (2017). A similar dedicated investigation in HMXBs is still required to evaluate how impacted the wind ionization structure can be by the spinning pulsar.

6.1.3 Eccentric orbits

For the sake of simplicity, we kept the number of degrees of freedom of our model limited to a minimum and discarded eccentric orbits. However, the present model can straightforwardly be extended to study a specific system where the orbit is suspected to be eccentric, especially when it comes to the fitting of the median N_H profile by the smooth wind solution (Eq. (10)). Although classic SgXBs such as Vela X-1 are largely circularized, SFXTs exhibit more eccentric orbits and would benefit from such an extension. For instance, the N_H profile obtained by García et al. (2018) in the SFXT IGR J16320-4751 displays a very sharp increase which the authors interpret as due to a grazing orbit. In addition, the flat N_H profile they measure away from this spike can only be reproduced either with a low inclination, incompatible with a grazing orbit, or with an eccentric orbit. They derived an eccentricity of 0.2, possibly indicative of the conditions of formation of the compact companion (Brandt & Podsiadlowski 1995). Even if they are more rare, eccentric orbits can occur in classic SgXBs too: in 2S0114+650 discussed in Sect. 3.1, some authors found an eccentricity close to 0.17 (Crampton et al. 1985; Grundstrom et al. 2007) although Koenigsberger et al. (2006) argued in favor of a zero eccentricity, which could account for the super-orbital period measured in this system.

6.1.4 Varying clump properties

A final limitation to be acknowledged is the unique clump mass and radius at r = 2 we consider. In reality, numerical simulations indicate various masses and radii for the clumps (Sundqvist et al. 2018). Ducci et al. (2009) considered the impact of a distribution of sizes and masses of the clumps to study the time variability of the mass accretion rate. We avoided introducing new parameters in the present study but such a distribution could be included in our simulations to seed the clumps. Similarly, it would be possible to study in more detail the impact of nonspherical clumps on the column density (Oskinova et al. 2012).

6.2 Observational challenges and prospects

In the Galaxy, we currently know ~20 wind-accreting SgXB and ~10 SFXTs withNS compact objects (Martínez-Núñez et al. 2017) as well as one wind-accreting SgXB with a BH, which constitutes a reasonable sample for observational studies. However, next to the model-intrinsic caveats listed above, our theoretical approach and suggested procedure do not yet take into account observational limitations and biases for measurements of N_H on the required time scales. An in detail discussion of these issues is out of scope of this paper, as it would require simulations for different instruments, observational set ups, source properties, and analysis approaches. We thus limit ourselves to giving some first estimates for observational challenges and prospects that could include the following.

• In regards to the possible cadence of X-ray measurements, X-ray satellites in low Earth orbit may not be able to provide uninterrupted measurements long enough to obtain meaningful measurements of the coherence time scale. As an example, we assume comparatively small clumps (R_c = 0.005 R_* and thus τ∕P ~ 8.6 × 10⁻⁴, cf. Fig. 10) in the black hole SgXB Cygnus X-1 (P = 5.6 d), resulting in a coherence time scale of ~400 s. This has to be compared to the typical length of an uninterrupted observation with a low Earth satellite such as RXTE of on average ~1.7 ks and a maximum of ~3.3 ks, that is only 4 and 8 full coherence cycles, respectively (Grinberg et al. 2015). Using such data for measurements would thus be challenging. On the other hand, instruments on elliptical orbits can potentially continuously observe a source for longer periods, up to ~100 ks (e.g., Haberl et al. 1989; Odaka et al. 2013; Miškovičová et al. 2016; Grinberg et al. 2017).

• It might be challenging to achieve a time resolution high enough for reliable absorption measurements. In the above example, we would need a time resolution of ~40 s to cover the coherence time scale with 10 measurements. The necessary time resolution would increase for larger clumps, toward ~600 s for R_c = 0.08 R_*, cf. Fig. 10, and be longer for systems with larger orbital periods. Such measurements are challenging. In the eclipsing HMXB Vela X-1, absorption has been directly modeled at different orbital phases on time scales of ~1 ks using Suzaku (Odaka et al. 2013) and XMM-Newton’s EPIC-pn (Martínez-Núñez et al. 2014) and on time scales of individual pulsar pulses, that is ~280 s, using NuSTAR (Fürst et al. 2014). The variability in 4U 1700−37 has been studied by Haberl et al. (1989) on time scales of ~1 ks. Typical formal errors on the measurements listed here will depend on the absolute value of N_H, but measurements of δN_H ~ 0.5…1 × 10²² cm⁻² are realistic (Odaka et al. 2013). However, direct absorption modeling from spectral fits on shorter time scales seems unrealistic for today’s instruments. With current instruments, the variability of absorption may be accessible on even shorter time scales using indirect methods such color-color diagrams (Nowak et al. 2011; Grinberg et al. 2020), although exact measurements are yet lacking for this approach. High throughput X-ray missions such as XRISM (XRISM Science Team 2020), eXTP (in ’t Zand et al. 2019), and Athena (Nandra et al. 2013) will enable even better measurements on short time scales and make a high number of fainter sources accessible for short-term absorption variability studies. For Athena high quality spectra of HMXBs can be obtained on time scales as short as 300 s (e.g., Lomaeva et al. 2020) and it is worth noting that estimates for the minimum time scale for absorption measurements with any of the listed future instruments have not yet been made.

• Absorption models used in X-ray astronomy and the assumed underlying continuum emission introduce systematic uncertainties in the measurements of the absorbing column density. This includes effects of cold vs. (lightly) ionized absorbers and different cross-sections and abundances between different models employed (i.e., the chosen opacity of the wind), but also known systematic effects such as the degeneracy between the spectral slope of the continuum emission and obtained N_H values for CCD-resolution instruments (e.g., Suchy et al. 2008). To our knowledge, there are no systematic studies of such effects for HMXBs yet that go beyond discussion of individual measurements and observations. An estimate of the overall reliability is thus not yet possible and in particular no comparison of such systematic uncertainties to the intrinsic variability of column density as predictedby our model.

Overall,the observational test of the presented model appears challenging, but possible.

7 Summary

Motivated by the pressing need for a better understanding of stellar outflows from hot stars, we attempted to bridge the gap between theory and observations of line-driven winds from massive stars in orbit with a compact object. We designed a model both complex enough to fairly account for results derived by numerical simulations of line-driven winds and simple enough to be of some use to interpret the observational data from HMXBs. We highlighted how the wind micro-structure could be constrained based on the variability of the column density. We provided the following standard procedure to determine the radius R_cl and mass m_cl of the overdense clumps formed in the wind of the massive donor star.

• We identified the fundamental parameters of the problem and illustrated how the stellar mass-loss rate could be estimated based on the scale of the orbital N_H profile averaged over several orbits.

• We showed that the standard deviation δN_H of the column density could be computed analytically based on the porosity length (Owocki & Cohen 2006), characterized here as the ratio between the clump column mass and the wind mean density. The analytical prediction is in agreement within a few percent with the measured values in the simulations based on the model we designed.

• In systems where the orbit is grazing, the maximum of δN_H is expected to occur at an orbital phase when the LOS is tangent to the region above the photosphere where the clumps start to form: the larger this orbital phase, the more extended the smooth wind region is. The measurements of the maximum of δN_H and the possible deviation of the measured δN_H curve from analytical predictions, thus provides a precious observational diagnostics to determine the location of this region.

• The standard deviation, δN_H, was found to be a good measure of the ratio $\sqrt{m_{\text{cl}}}/R_{\text{cl}}$. We also found that the coherence time scale of N_H time series was an excellent tracer of the clump self-crossing time at the point on the LOS the closest from the donor star, provided the relative speed between the clumps and the orbiting X-ray point source was accounted for. The coherence time scale is independent of the number of clumps per unit volume, which enables us to estimate the clump radius. Using the standard deviation δN_H, we can subsequently obtain an estimate on the mass of individual clumps without resorting on arguments using the time variability of the intrinsic X-ray emission.

We further discussed the current limitation of our model, both in terms of intrinsic limitations of the current simulation set-up and possible observational challenges. In practical terms, the suggested procedure will likely requires multiple X-ray exposures of a few 10 ks each of HMXBs with moderate luminosities and intermediate orbital inclinations, around orbital phases ϕ = 0.25 and at inferior conjunction (ϕ = 0.5) with sensitive X-ray instruments, that allow ≳100 individual measurements of absorbing column density in each exposure.

The heuristic model we designed is aimed at systems where the column density is essentially due to unaccreted clumps passing by the LOS, but it falls short of explaining several observations in HMXBs. Noticeably, the amount of mass accreted during X-ray flares is 10¹⁹⁻²³g, orders of magnitude larger than the individual clump mass obtained by radiative hydrodynamics simulations of the line-deshadowing instability. Also, the clump sizes deduced from the duration of enhanced absorption episodes are a few times larger than expected from first principles computations. This tension between observation and theory deserves to be addressed with additional investigations of the column density in HMXBs where the intermediate to low orbital inclination guarantees that the variability of N_H is likely due tothe clumps and not to large scale structures produced in the orbital plane by tidal forces.

With the advent of a new generation of X-ray satellites in the decade to come, time-resolved X-ray high resolution spectroscopy will grant us access to exquisite data. Valuable information on the morphology of the stellar wind in HMXBs, both at small andlarge scales, will be within our scope. Other domains will benefit from the upcoming modernization of the observing fleet. In the Active Galactic Nuclei community, several models of parametrized small scale obscurers have been designed to interpret the X-ray observations (Siebenmorgen et al. 2015; Hönig & Kishimoto 2017; Buchner et al. 2019). The origin of these structures remains largely unknown though, contrary to their clump counterparts in the wind of massive stars where they naturally form due to the line-deshadowing instability. Combined efforts to understand the propagation of X-rays in these multiphase environments will shed unprecedented light upon the structures surrounding accreting compact objects but also on the accretion process itself, strongly influenced by the ambient medium.

Acknowledgements

The authors warmly thank the referee for their particularly insightful report that improved both the scientific content and the readability of the papers, in particular for bringing up the possibility to constrain the basis of the clump-forming region from the δN_H profiles. I.E.M. wishes to thank Tomer Shenar for information concerning the system 2S0114+650, Antonios Manousakis for fruitful discussions about the underlying computational aspects, Lida Oskinova for insights concerning the modeling part and the members of the X-wind collaboration for the observational consequences of the present analysis. IEM has received funding from the ERC consolidator grant 646758 AEROSOL and from the Research Foundation Flanders (FWO) and the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 665501. V.G. is supported through the Margarete von Wrangell fellowship by the ESF and the Ministry of Science, Research and the Arts Baden-Württemberg. I.E.M. and J.O.S. are grateful for the hospitality of the International Space Science Institute (ISSI), Bern, Switzerland which sponsored a team meeting initiating a tighter collaboration between massive stars wind and X-ray binaries communities. The simulations were conducted on the Tier-1 VSC (Flemish Supercomputer Center funded by Hercules foundation and Flemish government). J.O.S. and F.A.D. acknowledge support by the Belgian Research Foundation Flanders (FWO) Odysseus program under grant number G0H9218N. M.A.L. acknowledges support from NASA’s Astrophysics Program.

References

All Tables

All Figures

Fig. 1 Top panel: clump radius (solid lines) and porosity length (dashed lines) in units of R* as a function of the distance to the stellar center. In blue (resp. in red) are represented the profiles for the clump radius expansion law (4) (resp. Eq. (5)). The dimensionless clump radius at r = 2, indicated with a black dot, was set to 0.005, and the dimensionless clump mass to 4 × 10−7 (see Sect. 2.4 for normalization units). Bottom panel: clumping factor profiles for both expansion laws. The red (resp. blue) shaded region represents the shape of the clumping factor profile obtained for O (resp. B) supergiants by Driessen et al. (2019). Spectroscopically, these OB supergiants can be considered representing typical early- and middle-type stars.

In the text

Fig. 2 Representation of the orbit (in green) of the compact object around its stellar companion (in blue) for Vela X-1 (top panel) and Cygnus X-1 (bottom panel) as seen from Earth, where the eccentricity of the orbit has been neglected. Only a fraction of the clumps susceptible to contribute to the column density along the orbit have been represented (in semi-opaquegrey). The black cross indicates the point of inferior conjunction on the orbit.

In the text

Fig. 3 Column density, per logarithmic radial bin (red, logarithmic axis on the left) and cumulative (blue, linear axis on the right), as a function of the distance to the stellar center. The position of the compact object is indicated with a dashed black line. The total column density up to 30 orbital separations is ~97% of the column density computed for a smooth wind integrated up to thousands of stellar radii. The dimensionless parameters are mcl = 10−5, Rcl = 0.08, a = 2, P = 80 and i = 45°.

In the text

Fig. 4 Three different evolutions of the column density as a function of time over several orbits for i = 90° (top panel), i = 45° (middle panel) and i = 0° (bottom panel) – and different clump radii, orbital separations, periods and β exponent.The dashed green line is the corresponding smooth wind column density profile. The top profile displays eclipses due to the high inclination. The insets show zoomed in portions of the curves where the numerical time sampling of the column density, 3 × 104 times per orbital period, is visible. The origin of time is set at inferior conjunction (i.e., t = 0 at ϕ = 0.5).

In the text

Fig. 5 Column density from superior conjunction (ϕ = 0) to inferior conjunction (ϕ = 0.5) for an orbital inclination i = 45°. In black is represented the median over 10 orbits and the red dashed line is the smooth wind solution from Eq. (10). The dark (resp. light) grey shaded region is bounded by the 20th to 80th percentiles (resp. the 5th to 95th percentiles). The parameters are mcl = 10−5, Rcl = 0.01, β = 0.5 and a = 2.

In the text

Fig. 6 Standard deviation of the column density as a function of the orbital phase computed numerically (in blue) and analytically (black dots) for simulations with mcl = 10−5, Rcl = 0.02. Top panel; β = 2, a = 2 and i = 90° and bottom panel: β = 0.5, a = 1.6 and i = 45°. The discrepancy at low orbital phase is partly indicative of the region near the stellar photosphere where the wind in the numerical model is assumed to not be clumpy yet.

In the text

Fig. 7 From top to bottom row: orbital inclination is 65°, 25° and 0°. Left column: standard deviation of the column density at inferior (resp. superior) conjunction in red (resp.blue) as a function of the ratio $\sqrt{m_{\text{cl}}}/R_{\text{cl},\text{2}}$. Each couple (a, β) leads to a different proportionality constant between δNH and $\sqrt{m_{\text{cl}}}/R_{\text{cl},\text{2}}$, hence the large dispersion. The dashed line in the top left panel indicates the y = x line. Right column: same as before but Rcl,2 has been replaced by the clump radius at the point along the LOS between the X-ray source and the observer the closest from the donor star. The dispersion is reduced, especially at inferior conjunction.

In the text

	Fig. 8 Autocorrelation function of the column density as a function of the time shift. From purple to green and yellow, we represent the autocorrelation functions of the data at different orbital phases, from inferior conjunction (dark purple) to superior conjunction (yellow). Simulation with mcl = 4 × 10−7, Rcl = 0.04, β = 2, a = 3, i = 65° and P = 85.
In the text

	Fig. 9 Stacked histograms of the ratios of the measured coherence time scale τ compared to the clump self-crossing time τSC for all 600 simulations at inferior conjunction (ϕ = 0.5, purple), ϕ = 0.25 (green) and superior conjunction (ϕ = 0, yellow).
In the text

	Fig. 10 Synthetic diagram summarizing the dependence of the coherence time scale τ and of the standard deviation on the column density δNH on the size of the clumps and their mass, for intermediate parameters (a = 2, i = 45°, P = 24 and β = 2). We used as normalization units m0 = 1024g for masses and NH,0 = 3 × 1023cm−2 for column densities.
In the text