Accretion from a clumpy massivestar wind
in Supergiant Xray binaries
Abstract
Supergiant Xray Binaries (sgxb) host a compact object, often a neutron star (NS), orbiting an evolved O/B star. Mass transfer proceeds through the intense linedriven wind of the stellar donor, a fraction of which is captured by the gravitational field of the NS. The subsequent accretion process onto the NS is responsible for the abundant Xray emission from sgxb. They also display peaktopeak variability of the Xray flux by a factor of a few 10 to 100, along with changes in the hardness ratios possibly due to varying absorption along the lineofsight. We use recent radiationhydrodynamic simulations of inhomogeneities (aka clumps) in the nonstationary wind of massive hot stars to evaluate their impact on the timevariable accretion process. For this, we run 3D hydrodynamic simulations of the wind in the vicinity of the accretor to investigate the formation of the bow shock and follow the inhomogeneous flow over several spatial orders of magnitude, down to the NS magnetosphere. In particular, we show that the impact of the wind clumps on the timevariability of the intrinsic mass accretion rate is severely tempered by the crossing of the shock, compared to the purely ballistic BondiHoyleLyttleton estimation. We also account for the variable absorption due to clumps passing by the lineofsight and estimate the final effective variability of the column density and mass accretion rate for different orbital separations. Finally, we compare our results to the most recent analysis of the Xray flux and the hardness ratio in Vela X1.
keywords:
accretion, accretion discs – Xrays: binaries – stars: neutron, supergiants, winds, outflows – methods: numerical1 Introduction
Massive, isolated supergiant O/B stars (Sg) have been thoroughly investigated in the last years under several complementary points of view : using asteroseismology (Aerts et al., 2017), investigating rotation properties and chemical abundance patterns (e.g. Martins et al., 2016), and not the least by analyses of their powerful linedriven wind outflows (see review by Puls et al., 2008). But while such studies have brought many significant insights regarding the evolution of massive stars, the recent detection of gravitational waves from coalescing compact objects (CO, Abbott et al., 2016a, b) have urged even more the community to evaluate the impact of binarity on the evolutionary tracks.
In particular, we expect to find good candidates as progenitors of compact binaries among the High Mass Xray Binaries (hmxb), systems where a stellarmass black hole or a neutron star (NS), accretes a fraction of the wind emitted by the stellar companion. This mass transfer mechanism called wind accretion is much less constrained than the one in Low Mass Xray Binaries, where it essentially proceeds through Roche Lobe Overflow (Paczynski, 1971). Until now, the BondiHoyleLyttleton (BHL) sketch (see the review by Edgar, 2004) applied to a smooth wind is used to estimate the mass accretion rate onto the compact object (Davidson & Ostriker, 1973). However, the underlying idealized representation, both of the wind and of the accreted flow, might jeopardize any attempt to accurately account for wind mass transfer in stellar evolutionary codes.
In the subclass of classic Supergiant Xray Binaries (sgxb), the evolved Sg companion undergoes dramatic mass losses through a fast, dense and isotropic linedriven wind : because of resonant line absorption of ultraviolet photons by partly ionized metal ions, the outer stellar layers are provided net outwards linear momentum (Lucy & Solomon, 1970; Castor et al., 1975). Such linedriven winds are known to be the siege of an intrinsic and remarkably strong instability (the linedeshadowing instability, LDI, Lucy & White, 1980; Owocki & Rybicki, 1984) whose development leads to internal shocks within the wind. The nonlinear growth of this instability has been computed first in onedimensional simulations (see e.g. Owocki et al., 1988; Feldmeier et al., 1997) but it is only in the 2000’s that the first multidimensional simulations were performed (Dessart & Owocki, 2005). These first multidimensional simulations witnessed a transverse fragmentation of the shells into higher density regions called "clumps", but could not properly resolve the lateral extension of the clumps. Using a twodimensional pseudoplanar grid sampling a restricted angular region, Sundqvist et al. (2017) (hereafter S17) recently managed to capture the dimensions of the clumps in both directions and to follow them up to two stellar radii. Their results agree with expectations from linear analysis (Rybicki et al., 1990) and are consistent with observations (Oskinova et al., 2007; Fürst et al., 2010). Investigating the way this structured wind can be beamed and captured by the gravitational field of the orbiting NS in a sgxb is the object of the present paper.
Previous numerical models of wind of a massive star with accretion onto an orbiting compact object have been made in Blondin et al. (1990) and Manousakis & Walter (2015). The authors used a global modeling approach of the wind and the compact object, accounting for the Xray ionizing feedback on the flow, which fully ionizes the metal ions and subsequently inhibits the acceleration (Hatchett & McCray, 1977). They also included the orbital effects, twisting and shearing the wind, but considered a smooth wind (suppressing the LDI by using a lineforce based on the socalled Sobolev approximation). On the contrary, Oskinova et al. (2012) and Bozzo et al. (2016) modeled the clumps but they relied on the LDI simulations by Feldmeier et al. (1997) which were unidimensional and thus did not resolve the transverse structure of the clumps ; the authors warned that the simple BHL analytic recipe for accretion then yields excessive variability, unconsistent with the observations of sgxb. More recently, CruzOsorio et al. (2017) worked out the influence of rigid and randomly distributed static spheres on the accretion of a planar supersonic homogeneous BHL flow in a twodimensional spherical slice where invariance around the polar axis is assumed. Ducci et al. (2009) adopted a statistic representation of spherical clumps in the wind to evaluate their impact on the timevariability of the mass accretion rate. However, they did not solve the dynamical problem and resorted on the BHL prescription to deduce the mass accretion rate as clumps intersect the transverse section of the accretion cylinder.
In the present paper, we consider a physicallymotivated threedimensional wind derived from full radiationhydrodynamic simulations. Its dynamics is not altered by the presence of a neutron star (NS) companion up to a few accretion radii around the accretor. Then, we solve the equations of hydrodynamics to compute the evolution of the inhomogeneous flow as it is beamed by the gravitational field of the NS within a few accretion radii and down to a few hundred times the size of the NS (i.e. at the outer edge of its magnetosphere). For different orbital separations and mass of the accretor, we compute, as a function of time, the mass accretion rate , and the column density in any direction. We use those quantities to estimate the impact of a realistic wind clumpiness on the time variability of the absorption and Xray luminosity.
In Section 2, we describe how we used the most recent simulations on the nonlinear growth of clumps. We then use this wind as outer upstream boundary condition for a numerical setup centered on the accretor and described in more detail in Section 3 and the appendices. The obtained results on and are discussed in Section 4 and conclusions are summarized in Section 5.
2 Clumpy wind representation
Numerical radiationhydrodynamic simulations of a linedriven wind generally show that the nonlinear growth of the LDI leads to highspeed rarefactions that steepen into strong reverse shocks, whereby the wind plasma becomes compressed into spatially narrow "clumps" separated by large regions of rarefied gas (Owocki et al., 1988; Sundqvist & Owocki, 2013; Feldmeier & Thomas, 2017). In previous 1D models, such clumps really are radial shells ; S17 recently performed 2D radiationhydrodynamic simulations within 2 stellar radii () which demonstrate how these characteristic LDIshells are broken up by the influence of oblique radiation rays and basic hydrodynamic instabilities. The width of the two dimensional simulation box, , is large enough to capture the transverse extension of the clumps, while still small enough to work in a pseudoplanar approach. After a few wind dynamical flow timescales, the multidimensional structures develop into a complex but statistically quite steady flow, characterized by localized highdensity clumps of small spatial extent embedded in larger regions of much lower density. The full numerical setup is described in detail in S17.
These simulations displayed a transition from an initially smooth analytic profile to a statistically steady state. As described above, the transition is due to the LDI (Owocki & Rybicki, 1984) whose growth first forms shells due to the internal shocks before lateral fragmentation occurs under the action of the lateral component of the lineforce. The obtained overdensities range from being sometimes roughly bowshaped to sometimes more isotropic, have typical sizes of a few at and are denser than the interclump environment by a factor up to a few 100 (see Figure 1). Their total mass is g, available to enhance the permanent mass accretion rate if such a clump would happen to pass close enough to the accretor. Those dimensions and masses are consistent with the observations (Fürst et al., 2010; Leutenegger et al., 2013; Grinberg et al., 2015) and theoretical estimates (Ducci et al., 2009), and they favor the lower end of the previously claimed values.
S17 considered an O star with parameters given in Table 1. The main difference between the star they considered and the typical Supergiant O/B host star in a classical sgxb (such as HD 77581 in Vela X1) is mostly the effective liberation speed at the photosphere, 20% lower in the current case. It has been shown that the predicted structure of the wind should be little altered by reasonable changes in those parameters.
In this paper, we use the results from S17 as outer boundary conditions to our simulation space centered on the accretor. We extend the results to 3D while still retaining the stochastic properties of the flow (histograms and correlation map) following the methods presented in more detail in Appendix A.
3 Numerical model
The numerical mesh we work with is a three dimensional extension of the one introduced in El Mellah & Casse (2015), a spherical stretched AMR mesh centered on the accretor whose outer spatial extension is described in Section 3.1 below. Within this space, we solve the conservative equations of hydrodynamics using the MPIAMRVAC finite volume code (Porth et al., 2014). To inject the highly supersonic clumpy wind described in the previous section in the simulation space, we incorporate timevariable outer boundary conditions (Section 3.3). For a more thorough description of the numerical setup itself, we refer to Appendix B.
3.1 Self defined accretion radius
To avoid any spurious reflection at the outer boundary downstream, once the flow has passed the accretor, we must make sure that the flow is supersonic once it reaches the outer edge (Blondin & Pope, 2009). The incoming wind is largely supersonic but gravitational beaming of streamlines with an impact parameter smaller than the accretion radius, leads to the formation of a detached bow shock (El Mellah & Casse, 2015) with given by :
(1) 
where is the gravitational constant, the mass of the accretor and the velocity of the flow before it is significantly altered by the gravitational field of the accretor^{1}^{1}1Not to be confused here with the wind terminal speed . (see eg Edgar, 2004). In the wake of the accretor (a.k.a. accretion trail), the flow which manages to escape accretion remains subsonic for a few accretion radii before being supersonic again (Blondin & Pope, 2009). A spherical simulation space with an outer radius of is enough to witness the formation of the shock without numerical artifacts due to the outer boundary conditions downstream which can simply be taken as continuous.
We compute in a low and high orbital separation configuration (respectively the close and wide configurations in Table 1). The effective velocity at the orbital separation, , is computed using the values of radial speed from S17 : once averaged transversally and in time around the orbital separation, , we obtain a speed representative of the amount of specific kinetic energy provided to the wind by the lineforce process before it starts to experience the gravitational grasp of the accretor. It serves as the normalization speed. We validated a posteriori that this estimate is valid and does give a suitable outer radial extension to the simulation space. We tune the NS mass so as to work with a similar physical accretion radius in both configurations. For realistic values of the mass of the accretor, of the wind speed at the orbital separation and of the stellar radius, the outer radius of the simulation space is a few times larger than the width of the simulation stripes in our LDI wind models, hence the need for the transverse prolongation described in Appendix A.1. Furthermore, this sphere of radius centered on the accretor coincides approximately with the Roche lobe of the compact object, which justifies to discard the gravitational influence of the compact object outside this region and to rely on the wind simulation from S17.
Parameters  Close config.  Wide config. 

Stellar mass  50M  
Stellar radius  20R  
Mass loss rate  1.310Myr  
Orbital separation  1.6R  2R 
NS mass  1.3M  2M 
Effective velocity at NS  925 kms  1140 kms 
Density at orb. sep.  7.210gcm  3.610gcm 
Accretion radius  cm R 
3.2 Ionization surface
What about the radiative influence of the Xrays emitted in the immediate vicinity of the accretor on the inflowing wind? The wind launching mechanism derived by Lucy & Solomon (1970) and Castor et al. (1975) for hot stars relies fully on the capacity of partly ionized metal ions to resonantly absorb ultraviolet photons. If those ions get to be fully ionized, no radiative momentum can be tapped by absorption and the acceleration of the wind is suddenly inhibited, an effect first reported by Hatchett & McCray (1977). The ionization parameter defines locally the state of ionization of the flow for an optically thin gas in local ionization and thermal balance :
(2) 
where is the Xray luminosity emitted by the accreting source, the distance to this source and the local atomic number density. The exact value of the critical ionization parameter above which the lineforce is significantly quenched has been a matter of debate, varying from ergcms (Hatchett & McCray, 1977; Ho & Arons, 1987) to ergcms (Blondin et al., 1990; Ducci et al., 2010; Manousakis & Walter, 2015) and even lower values (Stevens, 1991; Krticka et al., 2015). Oskinova et al. (2012) also showed that accounting for clumpiness lowers the ionization parameter compared to a smooth wind, by a factor of for clumping factors (aka wind inhomogeneity parameter) such as those obtained in our LDI simulations.
To evaluate the extension of this region in the upstream hemisphere of the region centered on the compact object, we compute the ionization parameter along a set of rays emitted from the accretor. The computation stops when reaches , which enables us to locate the ionization front for different ratios (Figure 2). Assuming optically thin conditions (Oskinova et al., 2012; Grinberg et al., 2015; MartínezNúñez et al., 2017), describes entirely the local ionization state of the flow. Figure 2 is plotted for ergs, consistent with the lowluminosity solution of Ho & Arons (1987), and . Those values are representative of sgxb. Depending on the directions and the time snapshot considered, the ionization front (red line in Figure 2) approximately matches the outer rim of the simulation space (black dashed arc in Figure 2). As a first approximation, we assume that they exactly coincide at any time, which means that the wind acceleration is fully inhibited within the simulation space and that the flow is unaltered downstream this region. Consequently, we can inject the clumps at the outer rim of the upstream hemisphere of the simulation space by using the timevarying values of mass density, speed and pressure from the LDI clumpy wind results of S17. In the present paper, we do not account for the timevarying Xray ionizing feedback from the immediate vicinity of the NS on the extension of this ionized region.
3.3 Outer boundary conditions
Once we set the orbital separation and the mass of the accretor, we can place the compact object within the wind and estimate the radial extension of the simulation space. Based on the three dimensional reconstruction described in Appendix A, we can determine the time varying outer boundary conditions downstream for the mass density, the velocity and the pressure. To do so, we sample the spherical mesh embedded in the Cartesian three dimensional data and assign values of the closest pixel. Such a straightforward interpolation method guarantees that we do not lower the contrasts determined by the LDI wind simulations. To be sure to resolve the clumps, in particular those with a small impact parameter, likely to be accreted, we realize this sampling with different mesh resolutions corresponding to the different Adaptive Mesh Refinement (AMR) levels described in Appendix B. An animation of the clumps entering the simulation space from the upstream hemisphere can be found at http://homes.esat.kuleuven.be/ileyk.
3.4 Preliminary relaxation
We first work with a planar and steady inflow to reach a numerically relaxed state and guarantee that the time variability measured in the simulations is fully induced by the inhomogeneities in the flow. The main computational cost of the simulations comes from (i) the large number of cells due to the full threedimensional geometry and (ii) the large number of time iterations required due to the large contrast between the inner boundary and the outer boundary, a few hundred times larger. To alleviate this difficulty and reach a numerically relaxed state at an affordable cost, we proceed in the following way. First, we work on a twodimensional radially stretched spherical grid similar to the one used in El Mellah & Casse (2015) with an inner boundary of radius (here and henceforth, length are in units of ). It is a larger inner boundary than the final one of , to reduce the number of time steps required to relax, but small enough to not significantly alter the detached bow shock formation. Then, the obtained axiallysymmetric relaxed state is resampled on a threedimensional spherical grid, rotated such that the polar axis is now orthogonal to the direction of the inflowing wind. We let this configuration relax in 3D and then, reduce the inner boundary radius down to its final value, . After a last stage of relaxation, we obtain a steady flow where the inner mass accretion rate varies at a level of less than 5%. We can now set the inhomogeneous and timevarying upstream outer boundary conditions described in Section 3.3 and let the simulations reach a statistically steady state.
4 Results & discussion
4.1 Structure of the flow
The structure of the flow follows qualitatively what has been computed with a planar homogeneous inflow (El Mellah & Casse, 2015). As the supersonic flow approaches the accretor, it is deflected towards its wake and forms a bow shock with a distance between the shock front and the accretor of approximately one accretion radius (see left panel in Figure 4). The inhomogeneities perturb the wake (see Figure 3). Also, the front shock moves back and forth but never reaches the inner boundary : the shock always remains detached.
Contrary to what has been reported from twodimensional simulations, we do not see any significant flipflop instability in the wake i.e. a wobbling motion of the instantaneous direction of the trailing flow. Historically, this instability has first been witnessed by Matsuda et al. (1987) in twodimensional numerical simulations and later confirmed and explored in more detail for different parameters (see e.g. Blondin & Pope, 2009). However, threedimensional simulations proved way more stable when precautions concerning the geometry of the grid and the size of the inner boundary were taken (Blondin & Raymer, 2012). Since we do not consider any orbital effect in the current paper, we cannot exclude the possibility that this would make a flipflop instability in the orbital plane possible. Within one accretion radius from the accretor, we do witness transient disclike structures (see Figure 6) but the alignment of their angular momentum with the polar axis suggests that they are essentially due to the numerical setup ; indeed, the planar BHL problem does not admit any preferential transverse direction but our present numerical setup introduces one, both because of the polar axis and because of the 3D reconstruction method (see Appendix A).
Concerning the wake itself, we measure a larger systematic opening angle of the bow shock compared to the steady homogeneous inflow configuration at high Mach number (see Figure 3). We see overdensities in the wake flowing away from the accretor. Their density contrast is much lower than in the upstream wind. Those overdensities form when a clump with an impact parameter too large (or, equivalently, with too much angular momentum) to be accreted reaches the shock, and propagates along the hull of the wake.
4.2 Mass accretion rate variability
We compute the mass accretion rate at the inner border (2 bottom panels in Figure 5), the mass flow rate entering through the outer upstream hemisphere (top panel) and the total net mass rate through the whole simulation outer border (second panel). The last serves as a proxy to estimate the moment when the simulation has reached a statistically steady state : it happens once the wind has crossed the whole simulation space, after approximately one hour. Once the numerically relaxed state for a homogeneous inflow has been reached (see Section 3.4), the mass accretion rate is constant to a 5% level or less.
Since clumps are typically 20 times smaller than the diameter of the simulation space, the inflowing mass rate through the outer upstream hemisphere remains fairly constant as a function of time (top panel in Figure 5). At a larger distance from the stellar surface, the clumps are larger and since the two configurations we consider have approximately the same accretion radius, it means that the clumps are comparatively larger for an accretor farther away from the star (wideconfiguration, red line in Figure 5). Hence a larger relative variability of the mass inflow rate, up to 20%. As a comparison, we compute the BHL mass accretion rates for a smooth wind :
(3) 
where is the density at the orbital separation (for an isotropic wind) and the 80% efficiency factor is measured during the preliminary relaxation steps (see Section 3.4) and in 2D simulations (El Mellah & Casse, 2015). We evaluate (dashed lines in bottom two panels in Figure 5) using the parameters from Table 1 and they approximately match the order of magnitude of the timeaveraged measured. However, such smoothwind BHL estimates fail to capture the time evolution of the mass rate through the inner edge of the simulation space, once the flow has been shocked. Let us now discuss this timevariability.
First, we examine the behavior of the flow at the shock. Sometimes, under the joint action of several clumps or of a particularly massive and cold one, the shock recedes down to the vicinity of the inner boundary. In this case, a clump can happen to get close to the inner boundary without undergoing much deformation. In Figure 6 for instance, although the front shock, delimited by the Mach1 contour (solid white line) stands firmly along the inflowing axis , it got breached on its side, at a distance of approximately one accretion radius from the accretor, by a stretched but unshocked clump. However, this clump has a slightly too high impact parameter and the angular momentum it carries is so important that it is not fully accreted ; it overshoots the accretor. Furthermore, at the same time, the other side of the shock is fairly unperturbed, with little inflow, and cannot feed angular momentum of opposite sign to facilitate accretion. This is a keyfeature of inhomogeneous BHL accretion : since the axisymmetry is broken, accretion is much less efficient, even within the accretion cylinder, because a laterally inflowing clump with no opposite axisymmetric clump to compensate for its angular momentum will not be accreted. It also explains why, when the first clumps reach the accretors, we saw a sharp drop in in all numerical setups. Matter piles up at the shock and in the wake but accretion itself is triggered only when the inflowing angular momentum changes sign. The most important flares are triggered by the conjunction between a series of clumps of low impact parameters. Then, the front of the shock recedes, the overdensities penetrate deep and contrary to the case above of high impact parameter, they are directly accreted, leading to a flaring by a factor of within a few minutes. As a consequence, the timevariability is enhanced with respect to the one of the mass inflowing rate from the upstream hemisphere at 8, but not as much as expected from the BHL argument where dissipation is idealized.
An important qualitative conclusion to draw from the simulations is that an incoming large clump does not straightforwardly yield a flare in . Conversely, it means that much caution should be taken when we try to trace back the origin of a flare to a clump mass. Between the ballistic clumps and the Xray emission regions lie many different regions which dilute the clumps, store matter and trigger their own instabilities and gating effects : time variability can be either enhanced or lowered depending on the physical conditions encountered by the flow as it gets accreted. There is no direct onetoone relation between the time variability at the smaller scales and the one induced by the inhomogeneous wind. Matter can pile up for quite long in intermediate regions (e.g. the shock or the corotation radius) before an instability which makes accretion possible is triggered.
4.3 Column density
The column density is the integrated number density along the lineofsight (LOS) :
(4) 
with the number density of particles. Throughout this paper, we relate to the mass density in the simulations by assuming a mean molecular weight of 1. Once the opacity is known, can be used to evaluate the absorption along the LOS at a given wavelength.
4.3.1 Preliminary computation
We evaluate the timevarying for an edgeon system using the following process. We put a constant Xray point source at 1.8 from the stellar center, embedded in the clumpy wind described above (see Figure 7). We work at orbital phases between and (when the NS moves towards us). Then, we differentiate two regions :

the simulation space, where we can compute onthego in 256 directions in the equatorial plane, synchronously with the mass rates in Figure 5.

the region centered on the stellar center and extending up to , deprived of the simulation space (henceforth "the outer region").
To reconstruct the wind in the outer region from the stripe displayed in Figure 1, we convert the 2D pseudoplanar coordinates into proper polar ones. Then, we stick together, next to each other, stripes corresponding to different time snapshots, so as to laterally extend the signal. If we assume a selfsimilar extension of the clumps beyond , we can now extrapolate the results up to by applying a homogeneous radial stretching factor to each stripe, computed such as the LOS of the observer (see semitransparent straight line in Figure 7) intersects the stripe. To scale the mass density so as to guarantee a conserved mass loss rate at all radii, we fit the velocity profile computed in the LDI wind simulations (between and ) with a law (Castor et al., 1975) :
(5) 
and assume that this profile is valid up to . is the distance to the stellar center, the terminal speed of the wind and and are the parameters we fit for. Since we work around , the LOS does not sample regions too distorted by this stretching. We computed along the LOS in this outer region as the NS orbits from to and as the clumps move outwards. Beyond the expected systematic decrease in column density as the compact object moves towards us, we see a peaktopeak spread of approximately cm and a correlation time of at most 1 kilosecond. It is longer than the selfcrossing time associated to a single clump, and results in a larger enhancement, because several independent clumps along the LOS can happen to act jointly. From now on, we use this value of cm as an estimation of the variability of in the outer region at any . As long as the LOS does not intercept neither the trailing flow in the wake of the NS nor the deep layers of the wind close to the photosphere, this estimate should hold since the microstructure of the wind is not significantly altered beyond a few accretion radii from the NS (see Section 3.2).
4.3.2 Column density within the simulation space
Concerning the column densities computed within the simulation space, the equatorial plane of our mesh is representative of any plane which contains the main axis of the shock (i.e. the direction of the inflow, in Figure 6). Indeed, the present model is axisymmetric in the sense that, contrary to the approach of El Mellah & Casse (2016), the source terms do not break the axisymmetry (e.g. no Coriolis force and thus, no bending of the stellar wind). This procedure enables us to capture, with a very good time resolution, the behavior of between a hundredth of accretion radius (i.e. the outer rim of the NS magnetosphere) and 8 accretion radii (of the order of a tenth of the orbital separation) around the compact object. The aim of this procedure is twofold : (i) evaluating the respective impact on the column density variability of ballistic clumps on one hand, and the inner structures within the shocked region on the other hand, and (ii) checking whether phases of high accretion rates are associated to phases of high .
As discussed in Section 2, the clumps injected into our accretion simulation space have quite low sizes and mass. Although their extension is not uniquely defined and their shape is not spherical, the overdense regions tend to occupy a zone of a fifth of the accretion radius and to weigh a few g. In terms of , when an unperturbed clump happens to pass along the LOS between the observer and the Xray point source, it translates into an enhanced of the order of cm during approximately 100 seconds, in agreement with what we measure in our simulations.
Clumps can, on one hand, feed the accretor and on the other hand, raise the column density, whether they are accreted or not : a clump passing by the LOS but with an impact parameter much larger than the accretion radius will induce an enhanced absorption but not an Xray flare. In Figure 8, we plot, as a function of the inner mass accretion rate and for both configurations, the column density averaged in all directions (top panels), seen from the side (middle panels, at or i.e. same viewing angle as Figure 4) or from the wake (bottom panels, i.e. inferior conjunction). A flare would be associated to an excursion to the right, sometimes associated to a rise in column density if the clumps responsible for this flare happened to pass by the LOS just before being accreted (see e.g. the emerging loop in the upper right corner of the right middle panel). The upper panels display a slight positive correlation between the episodes of enhanced accretion and of larger quantities of mass integrated around the inner boundary, with Pearson correlation coefficients ranging from 0.3 to 0.7. However, for given viewing angles (i.e. exposure times small compared to the orbital period), those results show that we should not observe a significant correlation between the column density and the Xray luminosity, in agreement with the latest observations by Bozzo et al. (2017).
4.3.3 Orbital column density of a smooth wind
Our simulations last a few kiloseconds at most, which represent only a few percent of the total orbital period for a system like Vela X1. If we suppose that the time variability of measured within the simulation space is representative of its evolution at any , we can simply add the median value for each angle to the corresponding as a function of phase computed from a smooth wind in the outer region. This smooth wind is the isotropic arithmetic and angular average in time of the clumpy wind model used throughout this paper. For the spreading, we use the quadratic sum of cm obtained in Section 4.3.1 and the standard deviation computed in each direction within the simulation space. It yields an estimate of the variability of for each , along with a median value.
The results are displayed in Figure 9 for the close and wide configurations. Since we do not account for any non axisymmetric effect in the present paper, the results between to are the same as the ones for to . On the left of the diagram is the eclipse region, at upper conjunction. Except in this region, the results for an edgeon inclination and a grazing inclination depart little from each other and thus, only the former is plotted. Due to the conservation of mass, the absolute values of the column densities scale approximately linearly with the stellar mass loss rate. Here, it is set to Myr in agreement with the value measured by Sako et al. (1999), Watanabe et al. (2006) and GimenezGarcia et al. (2016). Without surprise, we see that for larger distance of the accretor from the photosphere, the median at any orbital phase is lower. Concerning the time variability relative to the median , represented by the colored area, it depends essentially on the size of the clumps with respect to the accretion radius : a larger orbital separation yields larger clumps and a smaller accretion radius (since the velocity of the wind is larger) while a smaller mass of the accretor gives a smaller accretion radius (see Equation 1). Both effects would increase the spreading around the median trend. As long as it remains small compared to the stellar mass, the mass of the accretor, , has little influence on the average profile. Within the statistically axisymmetric wake, within and outside of the simulation space are comparable.
The black rectangle in Figure 9 represents recent results on in Vela X1 at . During this interval, Grinberg et al. (2017) observed two periods of one and two hours where the hardness ratio increased dramatically. A spectral analysis of the continuum led them to attribute this change to a sudden rise of from cm to cm (the upper and lower end of the black rectangle in Figure 9, along with the uncertainties). Both the close and wide configurations yield column densities lying within this range but they fail to explain this dramatic change in , which led Grinberg et al. (2017) to discard clumps happening to cross the LOS as possible explainations. Indeed, leaving aside the uncertainty on the median value (which can be shifted by tuning, for instance, the stellar mass loss rate), the present model is unable to reproduce the observed contrast in column density at these phases. It might be due to underlying structures within the magnetosphere that we could not study in detail with the present numerical setup.
4.4 Effective variability
The upper panel of Figure 10 shows, with solid lines, the histogram of at the inner boundary compared to its maximal measured value. The time fraction represents the time spent by the system at this level of accretion, compared to the physical simulation duration (8h). We discard the eclipses in this section. To guide the reader and discuss the results, we also show the observed luminosity diagrams for 4 classical sgxb in lower panel of Figure 10 : Vela X1, 4U 1700+37, OAO 1657415 and 4U 1907+09. These results are extracted from Figure 8 of Walter et al. (2015) where the eclipses and the offstates have been clipped off and the systems have been observed with Swift/BAT, in an energy range where little absorption is expected.
The first comment to make is that without absorption (solid lines in upper panel of Figure 10), we retrieve approximately lognormal distributions (Fürst et al., 2008, 2010; Walter et al., 2015), although cut at high activity. A serious discrepancy with observations concerns the restricted timevariability we see without absorption : for the closeconfiguration, varies by a factor of approximately 10 while for the wideconfiguration, this factor is . The timevariability of the hard Xrays emission in those systems is somewhat higher, between a few 10 and 100. It might be partly due to the limited duration of the simulations compared to the orbital period which might not have let the time to more extreme events (at the high and low ends of the activity level) to manifest. This can however not fully explain this discrepancy because the simulations of S17 did reach a statistically steady state, with a characteristic selfcoherent time scale of the order of a few hundred of seconds at most. Consequently, we should already capture a representative variability of the inflow over the few hours of our simulations and should only miss marginally extreme events. The small scale clumps considered here are thus unlikely to fully explain the high timevariability observed in these systems, even on longer duration than 8 hours. Note however that the 2D LDI simulations used here are quite likely to somewhat underestimate the quantitative level of clumping, and moreover that the model has been computed for a generic OBsupergiant in the Galaxy and not been specifically tuned to the parameters of Vela X1 (see discussion below). Other potential mechanisms for enhancing the timevariability include varying Xray ionizing feedback at the orbital scale (Manousakis & Walter, 2015), magneticcentrifugal gating at the outer rim of the NS magnetosphere (Illarionov & Sunyaev, 1975; Bozzo et al., 2008) or quasispherical subsonic settling accretion (Shakura & Postnov, 2017). Those mechanisms would possibly enable larger quantities of matter to be placed in "stand by" before being suddenly accreted once the physical trigger is matched. On the other hand, given their dimensions and their interaction with the bow shock surrounding the accretor, the clumps are not able to be stored for so long and flushed so suddenly. As a consequence, the wind inhomogeneities must be larger than those considered here in order to produce larger timevariability than a factor of 20 in at approximately 1,000 NS radii from the accretor.
But is this relatively low contrast between high and low mass accretion regimes in our simulations due to a lack of high or low activity levels? The absolute values of tend to favor the former option. Indeed, in the case of the closeconfiguration, the maximum accretion luminosity (with a maximum accretion efficiency of 1) obtained in Figure 5 amounts to ergs. For the wideconfiguration, it hardly reaches ergs. Since the conversion of into Xray emission is generally considered to have an efficiency of a few 10% (Bozzo et al., 2008) and that the peak luminosity in a system like Vela X1 is closer to a couple of ergs, we likely underestimate the maximum by an order of magnitude. In the closeconfiguration, a 2 solar masses accretor and a 30% slower wind (Sander et al., 2017) could explain this discrepancy. Since the timevariability depends essentially on the size of the clumps comparatively to the accretion radius, clumps 4 times larger would be required to display similar levels of timevariability. In this respect, it is important to point out here that the LDI wind models considered in this paper deliberately stabilize the wind base at the stellar photosphere (see discussion in S17). While this provides a perfect, controlled environment for the present pilotstudy of clumpy wind accretion, it also somewhat underestimates the level of clumping. Indeed, in the preliminary rotating 2D LDImodels by S17, the smallscale clumps are embedded in largerscale structures emerging directly from the stellar surface. Future work will examine if considering these structures will bring up the the timevariability to the observed levels. Finally, more realistic lower wind speeds for systems like Vela X1 (GimenezGarcia et al., 2016) would increase the accretion radius but the altered size of the clumps it would be associated to remains to be studied to make any conclusive statement about the relative size of the clumps compared to the accretion radius and, by then, about the amplitude of the timevariability of the mass accretion rate.
In any case, the frequency of high activity levels is indeed underestimated in the present simulations. On the other hand, the low activity levels are captured accurately and match the ones of Vela X1, OAO 1657415 and 4U 1907+09, once we assume maximum an order of magnitude higher (which would shift the histograms to the bottom left in the upper panel of Figure 10). Overall, this suggests that the timevariability observed at low luminosity in sgxb is essentially due to clumps generated in the linedriven wind of the donor star, whereas further studies are required to determine if this could be true also for the highluminosity events.
Even if we do not expect significant absorption in the photoenergy ranges observed in the bottom panel in Figure 10, we still evaluated the impact of absorption on these results, for illustrative purposes. In the soft Xrays (keV) though, the emission is strongly absorbed and the column densities previously measured could not be discarded. Absorption would naturally account for the asymmetries in the aforementioned lognormal distributions, with an enrichment of the low luminosity end. To compute a timevariability from the Xray emission which accounts for absorption, we rely on the time and phasedependent within the simulation space and on the dependent out of the simulation space, computed from a smooth wind (see Section 4.3). Notice that due to the size of the inner border of our simulation space, we do not account for the column density between one and a few hundreds times the accretor radius. Also, the column densities we computed are for a stellar mass loss rate of Myr and the uncertainty reported in the literature around this value is at least a factor of 2 (Sako et al., 1999; Watanabe et al., 2006; GimenezGarcia et al., 2016). would scale as the mass loss rate, provided the velocity profile is unaltered. So as to span the full half orbital period (from inferior to superior conjunction), typically a few days, we repeated by concatenating the function of time displayed in Figure 5, with random initial phases. We also work with an edgeon inclination. To visualize the impact of those uncertainties, we consider low and high absorption configurations with a photoionization crosssection of respectively cm and cm (Wilms et al., 2000). We also compare those results to a zeroabsorption case, representative of the hard Xray emission (above 10keV). Once corrected for absorption, we obtain effective mass accretion rates or "pseudoluminosities" and plot their histogram in upper panel of Figure 10 (dashed and dotted lines). We see that when absorption is accounted for, the low activity part of the diagram is indeed enriched, while the high activity region is essentially unaltered.
5 Summary & outlook
We ran 3D simulations of the accreted flow, centered on the accretor and with a mesh suitable to resolve the inhomogeneities in the upstream wind and follow them over almost 3 orders of magnitude in space. As the flow enters the sphere of influence of the NS, it is beamed toward it and forms a detached bow shock. The successive clumps cross this shock and, provided they loose enough energy and angular momentum, get accreted. As expected from the BHL approach of clump accretion (Ducci et al., 2009), the accretion proceeds mostly through the accretion cylinder. However, the present hydrodynamical simulations show that the shock lowers the timevariability with respect to the BHL estimation. Also, it introduces a time lag and a phase mixing since the shocked material associated to a clump is not straightforwardly accreted : it might be stored in a transient disclike structure before accretion of matter with opposite net angular momentum triggers effective accretion.
Comparing to the Xray luminosity diagrams at high energy, we showed that the wind microstructure computed by the clumpy wind models considered in this paper is not sufficient, per se, to retrieve the timevariability levels observed in classical sgxb such as Vela X1, 4U 170037, OAO 1657415 or 4U 1907+09 (Walter et al., 2015). The behavior at low luminosity matches the observations but the largest luminosity events require the possibility to quickly tap larger amounts of matter the clumps we considered here are not able to provide, even accounting for the intermediate shocked region where the flow can pile up. Other storage stages can appear once the NS magnetosphere and radiative cooling is accounted for (Illarionov & Sunyaev, 1975; Bozzo et al., 2008; Shakura & Postnov, 2017) or at the orbital scale (e.g. the Xray ionizing feedback in Manousakis & Walter, 2015). However, as discussed above, while the LDI simulations used here provide a perfect testbed for this first study, such simulations do indeed yield larger clumps at a given orbital separation for a larger star and/or when clumping is enabled earlier on, near the stellar surface (where we know that clumps are already present, Cohen et al., 2011; Sundqvist & Owocki, 2013; Torrejón et al., 2015). Building on this first, generic study of 3D wind accretion from a linedriven clumpy wind outflow, future work will thus use LDI simulations more specifically tuned to Vela X1 to examine in more detail how well the clumpy accretion can explain also the recent observations by Grinberg et al. (2017).
Acknowledgments
The authors are indebted to the anonymous referee who brought up several insightful questions and helped to improve this paper. IEM has received funding from the Research Foundation Flanders (FWO) and the European Union’s Horizon 2020 research and innovation programme under the Marie SkłodowskaCurie grant agreement No 665501. IEM and JOS 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 Xray binaries communities. IEM also thanks Peter Kretschmar, Victoria Grinberg and Felix Fürst for the fruitful discussions and the relevant comments they made on the present work. This research was supported by FWO and by KU Leuven Project No. GOA/2015014 and by the Interuniversity Attraction Poles Programme by the Belgian Science Policy Office (IAP P7/08 CHARM). The simulations were conducted on the Tier1 VSC (Flemish Supercomputer Center funded by Hercules foundation and Flemish government).
References
 (1)
 Abbott et al. (2016a) Abbott, B. P., Abbott, R., Abbott, T. D. & Al, E. (2016a), Phys. Rev. Lett. 116(6), 061102.
 Abbott et al. (2016b) Abbott, B. P., Abbott, R., Abbott, T. D. & Al, E. (2016b), Phys. Rev. Lett. 116(24), 241103.
 Aerts et al. (2017) Aerts, C., SimonDiaz, S., Bloemen, S., Debosscher, J., Papics, P. I., Bryson, S., Still, M., Moravveji, E., Williamson, M. H., Grundahl, F., Andersen, M. F., Antoci, V., Palle, P. L., ChristensenDalsgaard, J. & Rogers, T. M. (2017), Astron. Astrophys. 602, 1–14.
 Blondin et al. (1990) Blondin, J. M., Kallman, T., Fryxell, B. & Taam, R. E. (1990), Astrophys. J. 356, 591–608.
 Blondin & Pope (2009) Blondin, J. M. & Pope, T. C. (2009), Astrophys. J. 700(1), 95–102.
 Blondin & Raymer (2012) Blondin, J. M. & Raymer, E. (2012), Astrophys. J. 752(1), 30.
 Bozzo et al. (2017) Bozzo, E., Bernardini, F., Ferrigno, C., Falanga, M., Romano, P. & Oskinova, L. (2017), eprint arXiv:1709.00716 .
 Bozzo et al. (2008) Bozzo, E., Falanga, M. & Stella, L. (2008), Astrophys. J. 683(2), 1031–1044.
 Bozzo et al. (2016) Bozzo, E., Oskinova, L., Feldmeier, A. & Falanga, M. (2016), Astron. Astrophys. Vol. 589, id.A102, 11 pp. 589.
 Castor et al. (1975) Castor, J. I., Abbott, D. C. & Klein, R. I. (1975), Astrophys. J. 195, 157.
 Cohen et al. (2011) Cohen, D. H., Gagné, M., Leutenegger, M. A., MacArthur, J. P., Wollman, E. E., Sundqvist, J. O., Fullerton, A. W. & Owocki, S. P. (2011), Mon. Not. R. Astron. Soc. 415(4), 3354–3364.
 CruzOsorio et al. (2017) CruzOsorio, A., SanchezSalcedo, F. J. & LoraClavijo, F. D. (2017), eprint arXiv:1707.05548 .
 Davidson & Ostriker (1973) Davidson, K. & Ostriker, J. P. (1973), Astrophys. J. 179, 585.
 Dessart & Owocki (2005) Dessart, L. & Owocki, S. P. (2005), Astron. Astrophys. 437(2), 657–666.
 Ducci et al. (2009) Ducci, L., Sidoli, L., Mereghetti, S., Paizis, A. & Romano, P. (2009), Mon. Not. R. Astron. Soc. 398(4), 2152–2165.
 Ducci et al. (2010) Ducci, L., Sidoli, L. & Paizis, A. (2010), Mon. Not. R. Astron. Soc. 408(3), 1540–1550.
 Edgar (2004) Edgar, R. G. (2004), New Astron. Rev. 48(10), 843–859.
 El Mellah & Casse (2015) El Mellah, I. & Casse, F. (2015), Mon. Not. R. Astron. Soc. 454(3), 2657–2667.
 El Mellah & Casse (2016) El Mellah, I. & Casse, F. (2016), Mon. Not. R. Astron. Soc. 467(3), 2585–2593.
 Feldmeier et al. (1997) Feldmeier, A., Puls, J. & Pauldrach, A. W. A. (1997), Astron. Astrophys. 322, 878–895.
 Feldmeier & Thomas (2017) Feldmeier, A. & Thomas, T. (2017), Mon. Not. R. Astron. Soc. 469(3), 3102–3107.
 Fürst et al. (2010) Fürst, F., Kreykenbohm, I., Pottschmidt, K., Wilms, J., Hanke, M., Rothschild, R. E., Kretschmar, P. & Schulz, N. S. (2010), 37, 1–8.
 Fürst et al. (2008) Fürst, F., Kreykenbohm, I., Wilms, J., Kretschmar, P., Klochkov, D., Santangelo, A. & Staubert, R. (2008), "Proceedings 7th Integr. Work. 8  11 Sept. 2008 Copenhagen, Denmark."" pp. 1–9.
 GimenezGarcia et al. (2016) GimenezGarcia, A., Shenar, T., Torrejon, J. M., Oskinova, L., MartinezNunez, S., Hamann, W.R., RodesRoca, J. J., GonzalezGalan, A., AlonsoSantiago, J., GonzalezFernandez, C., Bernabeu, G. & Sander, A. (2016), Astron. Astrophys. 591(A26), 25.
 Grinberg et al. (2017) Grinberg, V., Hell, N., El Mellah, I., Neilsen, J., Sander, A. A. C., Leutenegger, M. A., Fürst, F., Huenemoerder, D. P., Kretschmar, P., Kühnel, M., MartinezNunez, S., Niu, S., Pottschmidt, K., Schulz, N. S., Wilms, J. & Nowak, M. A. (2017), submitted .
 Grinberg et al. (2015) Grinberg, V., Leutenegger, M. A., Hell, N., Pottschmidt, K., Böck, M., García, J. A., Hanke, M., Nowak, M. A., Sundqvist, J. O., Townsend, R. H. D. & Wilms, J. (2015), Astron. Astrophys. 576.
 Harten et al. (1983) Harten, A., Lax, P. D. & van Leer, B. (1983), SIAM Rev. pp. 35–61.
 Hatchett & McCray (1977) Hatchett, S. & McCray, R. (1977), Astrophys. J. 211, 552.
 Ho & Arons (1987) Ho, C. & Arons, J. (1987), Astrophys. J. 316, 283.
 Illarionov & Sunyaev (1975) Illarionov, A. F. & Sunyaev, R. A. (1975), Astron. Astrophys. 39.
 Koren (1993) Koren, B. (1993), Notes Numer. Fluid Mech. 45, 117–138.
 Krticka et al. (2015) Krticka, J., Kubat, J. & Krtickova, I. (2015), Astron. Astrophys. 579, 16.
 Leutenegger et al. (2013) Leutenegger, M. A., Cohen, D. H., Sundqvist, J. O. & Owocki, S. P. (2013), Astrophys. Journal, Vol. 770, Issue 1, Artic. id. 80, 17 pp. (2013). 770.
 Lucy & Solomon (1970) Lucy, L. B. & Solomon, P. M. (1970), Astrophys. J. 159, 879.
 Lucy & White (1980) Lucy, L. B. & White, R. L. (1980), Astrophys. J. 241, 300.
 Manousakis & Walter (2015) Manousakis, A. & Walter, R. (2015), Astron. Astrophys. 58, A58.
 MartínezNúñez et al. (2017) MartínezNúñez, S., Kretschmar, P., Bozzo, E., Oskinova, L. M., Puls, J., Sidoli, L., Sundqvist, J. O., Blay, P., Falanga, M., Fürst, F., GímenezGarcía, Á., Kreykenbohm, I., Kühnel, M., Sander, A., Torrejón, J. M. & Wilms, J. (2017), Sp. Sci Rev .
 Martins et al. (2016) Martins, F., SimonDiaz, S., Barba, R. H., Gamen, R. C., Ekstroem, S., SimónDíaz, S., Barbá, R. H., Gamen, R. C. & Ekström, S. (2016), Astron. Astrophys. Vol. 599, id.A30, 14 pp. 30.
 Matsuda et al. (1987) Matsuda, T., Inoue, M., Sawada, K., Matsuda, T., Inoue, M. & Sawada, K. (1987), Mon. Not. R. Astron. Soc. (ISSN 00358711) 226, 785–811.
 Oskinova et al. (2012) Oskinova, L. M., Feldmeier, A. & Kretschmar, P. (2012), Mon. Not. R. Astron. Soc. Vol. 421, Issue 4, pp. 28202831. 421, 2820–2831.
 Oskinova et al. (2007) Oskinova, L. M., Hamann, W.R. & Feldmeier, A. (2007), Astron. Astrophys. 476(3), 1331–1340.
 Owocki et al. (1988) Owocki, S. P., Castor, J. I. & Rybicki, G. B. (1988), Astrophys. J. 335, 914.
 Owocki & Rybicki (1984) Owocki, S. P. & Rybicki, G. B. (1984), Astrophys. J. 284, 337.
 Paczynski (1971) Paczynski, B. (1971), Annu. Rev. Astron. Astrophys. 9(1), 183–208.
 Porth et al. (2014) Porth, O., Xia, C., Hendrix, T., Moschou, S. P. & Keppens, R. (2014), Astrophys. J. Suppl. Ser. 214(1), 4.
 Puls et al. (2008) Puls, J., Vink, J. S. & Najarro, F. (2008), Astron. Astrophys. Rev. 16(34), 209–325.
 Rybicki et al. (1990) Rybicki, G. B., Owocki, S. P. & Castor, J. I. (1990), Astrophys. J. 349, 274.
 Sako et al. (1999) Sako, M., Liedahl, D. A., Kahn, S. M. & Paerels, F. (1999), Astrophys. Journal, Vol. 525, Issue 2, pp. 921934. 525(1997), 921–934.
 Sander et al. (2017) Sander, A. A. C., Fürst, F., Kretschmar, P., Oskinova, L. M., Todt, H., Hainich, R., Shenar, T. & Hamann, W.R. (2017), eprint arXiv:1708.02947 (2016).
 Shakura & Postnov (2017) Shakura, N. & Postnov, K. (2017), eprint arXiv:1702.03393 .
 Stevens (1991) Stevens, I. R. (1991), Astrophys. J. 379, 310.
 Sundqvist & Owocki (2013) Sundqvist, J. O. & Owocki, S. P. (2013), Mon. Not. R. Astron. Soc. 428(2), 1837–1844.
 Sundqvist et al. (2012) Sundqvist, J. O., Owocki, S. P., Cohen, D. H., Leutenegger, M. A. & Townsend, R. H. D. (2012), Mon. Not. R. Astron. Soc. 420(2), 1553–1561.
 Sundqvist et al. (2017) Sundqvist, J. O., Owocki, S. P. & Puls, J. (2017), (2).
 Torrejón et al. (2015) Torrejón, J. M., Schulz, N. S., Nowak, M. A., Oskinova, L., RodesRoca, J. J., Shenar, T. & Wilms, J. (2015), Astrophys. J. 810(2), 102.
 Walter et al. (2015) Walter, R., Lutovinov, A. A., Bozzo, E. & Tsygankov, S. S. (2015), Astron. Astrophys. Rev. 23(1), 2.
 Watanabe et al. (2006) Watanabe, S., Sako, M., Ishida, M., Ishisaki, Y., Kahn, S. M., Kohmura, T., Nagase, F., Paerels, F. & Takahashi, T. (2006), Astrophys. J. 651(1), 421–437.
 Wilms et al. (2000) Wilms, J., Allen, A. & McCray, R. (2000), Astrophys. J. 542(2), 914–924.
 Xia et al. (2017) Xia, C., Teunissen, J., El Mellah, I., Chané, E. & Keppens, R. (2017), submitted .
Appendix A 3D reconstruction
To fully capture the dynamics of the flow around the compact object, three dimensional simulations of the accretion process centered on the accretor and spanning several orders of magnitude in space are required. Time varying outer boundary conditions are naturally provided by the 2D clumpy wind simulations presented in the main text, provided we develop a reconstruction method to expand these results along the second transverse dimension. In order to minimize the alteration of the results from S17, we use the following method. Let us write the information in each pixel of the 2D stripe at a given time (with for the radial direction and the transverse direction). Since the two transverse directions are equivalent, we compute a 3D signal by writing :
(6) 
where the over line indicates that we combine and either with a geometrical mean (for mass density and pressure) or with an arithmetic mean (for the radial velocity component). Concerning the transverse components of the velocity field, they are not retained in the present paper where we discard the shearing of the flow and the orbital effects. It yields the initial signal along the transverse diagonal () and a symmetric matrix for any index (). We then rotate by 45 degrees around the radial axis to retrieve the initial signal in one of the two transverse direction and to make use of the aforementioned matrix symmetry and work only in the upper half of the simulation space. Finally, since the averaging procedure lowers the contrasts in the other transverse direction, we proceed to a histogram correction to retrieve approximately the same histograms for the mass density, the pressure and the speed in both transverse directions.
The benefits of this approach are to conserve the transverse extension of the clumps, their longitudinal correlation length, their overall shape and to retrieve exactly the initial two dimensional signal along one transverse direction. On the contrary, the main drawback which must be acknowledged for is the spurious leakage along the two privileged orthogonal transverse directions associated to the matrices above : crossshaped patterns appear in the transverse slices of the three dimensional reconstructed signals. Overall, this technique tends to produce spherical clumps (see Sundqvist et al. 2012).
a.1 Transverse prolongation
Since the width of the stripes () in the clumpy wind simulations turns out to be smaller than the required simulation space centered on the compact object (of radius , see Section 3.1), we need to enlarge their transverse extension. Only a reduced fraction of the flow, the one with the smallest impact parameter () with respect to the accretor, is expected to actually be accreted at some point (El Mellah & Casse 2015). Beyond a few accretion radii, the flow must still be described within the simulation space but its precise microstructure will not play a major role in the time variability of the inner mass accretion rate. It it thus legitimate to simply repeat the central three dimensional data in both transverse directions. Thanks to the periodic transverse boundary conditions used in the original wind simulations, this does not introduce any discontinuity.
Appendix B Numerical setup
b.1 Multiscale mesh
The mesh introduced here to afford numerical simulations of clumpy wind being accreted is partly displayed in Figure 12 and described in more detail in Xia et al. (2017). The core physical requirements are (i) to be able to follow the flow from down to the shock () and below to resolve the sonic surface (a few hundredths of ) and (ii) to be able to resolve the wind microstructure at any time, in particular in the regions susceptible to be accreted.
b.1.1 Radially stretched
To uniformly probe the flow all along its journey to the accretor, it is not possible to work with a constant radial step. We designed a radially stretched mesh (El Mellah & Casse 2015) to make it affordable and to keep the same cell aspect ratio from the inner boundary of the grid up to the outer boundary. Doing so, we could span almost 3 orders of magnitude in space. It enabled us to properly compute the shock formation (3.1) and to study the accretion process within the shocked region, down to an inner boundary radius of a hundredth of the accretion radius. Larger inner boundaries lead to a significant alteration of the shock which becomes attached to the inner boundary. In physical units, the spatial extension of this inner boundary is approximately 1,000 times the size of the compact object ; for a neutron star, it corresponds to the region where the magnetic field role can no longer be neglected and where the hydrodynamic framework we rely on in this paper would break down.
b.1.2 The pole
Due to the shrinking of the azimuthal spatial step towards the poles and to respect the CourantLewyFriedrich condition, three dimensional spherical simulations suffer strong limitations on the time step at the pole. In Blondin & Raymer (2012), the authors intertwined two spherical meshes rotated by 90 degrees around one of the two equatorial axis, so as to avoid the polar regions of each mesh (the YinYang approach). Here, we simply degraded the resolution around the pole by preventing any AMR to take place within a certain angular region around the pole. Doing so, the grid remains at its first resolution level, typically , and refining the mesh in the equatorial regions does not entail smaller time steps up to an AMR level of 6 : the only additional cost when refining is the higher number of cells, which can easily be handled given the excellent weak scaling of MPIAMRVAC (Porth et al. 2014). The price to pay for such an approach is a lower resolution at the poles and a risk that artificial smoothing of the flow along the pole favors spurious disc formation in the equatorial plane of the mesh. We work with 4 different levels of AMR such that the effective resolution is of .
b.1.3 Selective AMR
The clumps enter the spherically stretched mesh from the upstream hemisphere. Because those clumps are offcentered small scale structures, we need to enable AMR (up to 4 levels) in this hemisphere to resolve them, in particular at the outer edge of the mesh where the radially stretched cells have the largest absolute size on the first AMR level. On the contrary, because we are only interested in the fraction of the flow susceptible to be accreted by the compact object, we inhibit AMR refinement in the downstream hemisphere, except in the immediate vicinity of the accretor (below the stagnation point in the wake). We also prevent excessive refinement in the vicinity of the accretor (no refinement beyond the third level) since the stretching already provides more refined cells. Finally, in the upstream hemisphere, we favor AMR refinement in the accretion cylinder, around the axis of zero impact parameter. All those precautions enable us to follow the flow as it is accreted onto the compact object while still resolving offcentered inhomogeneities.
b.2 Code and solving scheme
We solve the adiabatic equations of hydrodynamics under their conservative form using the finite volume code MPIAMRVAC (Porth et al. 2014). The energy equation contains no cooling nor heating term, and the gravitational field of the compact object is the only source term we consider. We use a shockcapturing HLL solver (Harten et al. 1983), and a Koren slope limiter (Koren 1993). The Courant number is typically set to 0.5. We also wrote a new userdefined type of boundary condition to work with timevarying boundary conditions stored in data files. A testcase of timevarying boundary conditions is now part of the main branch of the code.