Accretion of a clumpy wind in SgXB

# Accretion from a clumpy massive-star wind in Supergiant X-ray binaries

I. El Mellah, J.O. Sundqvist, R. Keppens
Centre for mathematical Plasma Astrophysics, Department of Mathematics, KU Leuven, Celestijnenlaan 200B, B-3001 Leuven, Belgium
KU Leuven, Instituut voor Sterrenkunde, Celestijnenlaan 200D, B-3001 Leuven, Belgium
E-mail: ileyk.elmellah@kuleuven.be
Accepted XXX. Received YYY; in original form ZZZ
###### Abstract

Supergiant X-ray Binaries (sgxb) host a compact object, often a neutron star (NS), orbiting an evolved O/B star. Mass transfer proceeds through the intense line-driven 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 X-ray emission from sgxb. They also display peak-to-peak variability of the X-ray flux by a factor of a few 10 to 100, along with changes in the hardness ratios possibly due to varying absorption along the line-of-sight. We use recent radiation-hydrodynamic simulations of inhomogeneities (aka clumps) in the non-stationary wind of massive hot stars to evaluate their impact on the time-variable 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 time-variability of the intrinsic mass accretion rate is severely tempered by the crossing of the shock, compared to the purely ballistic Bondi-Hoyle-Lyttleton estimation. We also account for the variable absorption due to clumps passing by the line-of-sight 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 X-ray flux and the hardness ratio in Vela X-1.

###### keywords:
accretion, accretion discs – X-rays: binaries – stars: neutron, supergiants, winds, outflows – methods: numerical
pubyear: 2017pagerange: Accretion from a clumpy massive-star wind in Supergiant X-ray binariesB.2

## 1 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 line-driven 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 X-ray Binaries (hmxb), systems where a stellar-mass 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 X-ray Binaries, where it essentially proceeds through Roche Lobe Overflow (Paczynski, 1971). Until now, the Bondi-Hoyle-Lyttleton (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 sub-class of classic Supergiant X-ray Binaries (sgxb), the evolved Sg companion undergoes dramatic mass losses through a fast, dense and isotropic line-driven wind : because of resonant line absorption of ultra-violet photons by partly ionized metal ions, the outer stellar layers are provided net outwards linear momentum (Lucy & Solomon, 1970; Castor et al., 1975). Such line-driven winds are known to be the siege of an intrinsic and remarkably strong instability (the line-deshadowing instability, LDI, Lucy & White, 1980; Owocki & Rybicki, 1984) whose development leads to internal shocks within the wind. The non-linear growth of this instability has been computed first in one-dimensional simulations (see e.g. Owocki et al., 1988; Feldmeier et al., 1997) but it is only in the 2000’s that the first multi-dimensional simulations were performed (Dessart & Owocki, 2005). These first multi-dimensional 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 two-dimensional pseudo-planar 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 X-ray 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 line-force based on the so-called 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 uni-dimensional 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, Cruz-Osorio 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 two-dimensional 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 time-variability 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 physically-motivated three-dimensional wind derived from full radiation-hydrodynamic 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 X-ray luminosity.

In Section 2, we describe how we used the most recent simulations on the non-linear 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 radiation-hydrodynamic simulations of a line-driven wind generally show that the non-linear growth of the LDI leads to high-speed 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 radiation-hydrodynamic simulations within 2 stellar radii () which demonstrate how these characteristic LDI-shells 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 pseudo-planar approach. After a few wind dynamical flow time-scales, the multi-dimensional structures develop into a complex but statistically quite steady flow, characterized by localized high-density 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 line-force. The obtained overdensities range from being sometimes roughly bow-shaped to sometimes more isotropic, have typical sizes of a few at and are denser than the inter-clump 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 X-1) 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 MPI-AMRVAC 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 time-variable 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 :

 Racc=2GMv2∞ (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 accretor111Not 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 line-force 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.

### 3.2 Ionization surface

What about the radiative influence of the X-rays 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 ultra-violet 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 :

 ζ=LXnr2 (2)

where is the X-ray 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 line-force 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ínez-Núñez et al., 2017), describes entirely the local ionization state of the flow. Figure 2 is plotted for ergs, consistent with the low-luminosity 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 time-varying 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 time-varying X-ray 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 three-dimensional 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 two-dimensional 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 axially-symmetric relaxed state is resampled on a three-dimensional 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 time-varying 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 two-dimensional simulations, we do not see any significant flip-flop 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 two-dimensional numerical simulations and later confirmed and explored in more detail for different parameters (see e.g. Blondin & Pope, 2009). However, three-dimensional 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 flip-flop instability in the orbital plane possible. Within one accretion radius from the accretor, we do witness transient disc-like 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 (wide-configuration, 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 :

 ˙MBHL=0.8πR2accρ∞v∞ (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 time-averaged measured. However, such smooth-wind 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 time-variability.

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 Mach-1 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 key-feature 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 time-variability 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 X-ray 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 one-to-one 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 line-of-sight (LOS) :

 NH=∫LOSn⋅dl (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 time-varying for an edge-on system using the following process. We put a constant X-ray 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 on-the-go 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 pseudo-planar 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 self-similar 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 semi-transparent 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) :

 v(r)=v0(1−R∗r)βwhere{v0∼2,240km⋅s−1β∼0.89 (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 peak-to-peak spread of approximately cm and a correlation time of at most 1 kilosecond. It is longer than the self-crossing 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 micro-structure 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 X-ray 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 X-ray 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 X-ray 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 X-1. 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 edge-on 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 Gimenez-Garcia 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 X-1 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 X-1, 4U 1700+37, OAO 1657-415 and 4U 1907+09. These results are extracted from Figure 8 of Walter et al. (2015) where the eclipses and the off-states 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 log-normal distributions (Fürst et al., 2008, 2010; Walter et al., 2015), although cut at high activity. A serious discrepancy with observations concerns the restricted time-variability we see without absorption : for the close-configuration, varies by a factor of approximately 10 while for the wide-configuration, this factor is . The time-variability of the hard X-rays 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 self-coherent 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 time-variability 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 OB-supergiant in the Galaxy and not been specifically tuned to the parameters of Vela X-1 (see discussion below). Other potential mechanisms for enhancing the time-variability include varying X-ray ionizing feedback at the orbital scale (Manousakis & Walter, 2015), magnetic-centrifugal gating at the outer rim of the NS magnetosphere (Illarionov & Sunyaev, 1975; Bozzo et al., 2008) or quasi-spherical 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 time-variability 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 close-configuration, the maximum accretion luminosity (with a maximum accretion efficiency of 1) obtained in Figure 5 amounts to ergs. For the wide-configuration, it hardly reaches ergs. Since the conversion of into X-ray 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 X-1 is closer to a couple of ergs, we likely underestimate the maximum by an order of magnitude. In the close-configuration, a 2 solar masses accretor and a 30% slower wind (Sander et al., 2017) could explain this discrepancy. Since the time-variability 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 time-variability. 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 pilot-study of clumpy wind accretion, it also somewhat underestimates the level of clumping. Indeed, in the preliminary rotating 2D LDI-models by S17, the small-scale clumps are embedded in larger-scale structures emerging directly from the stellar surface. Future work will examine if considering these structures will bring up the the time-variability to the observed levels. Finally, more realistic lower wind speeds for systems like Vela X-1 (Gimenez-Garcia 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 time-variability 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 X-1, OAO 1657-415 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 time-variability observed at low luminosity in sgxb is essentially due to clumps generated in the line-driven wind of the donor star, whereas further studies are required to determine if this could be true also for the high-luminosity events.

Even if we do not expect significant absorption in the photo-energy 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 X-rays (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 log-normal distributions, with an enrichment of the low luminosity end. To compute a time-variability from the X-ray emission which accounts for absorption, we rely on the time and phase-dependent 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; Gimenez-Garcia 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 edge-on inclination. To visualize the impact of those uncertainties, we consider low and high absorption configurations with a photoionization cross-section of respectively cm and cm (Wilms et al., 2000). We also compare those results to a zero-absorption case, representative of the hard X-ray emission (above 10keV). Once corrected for absorption, we obtain effective mass accretion rates or "pseudo-luminosities" 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 time-variability 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 disc-like structure before accretion of matter with opposite net angular momentum triggers effective accretion.

Comparing to the X-ray luminosity diagrams at high energy, we showed that the wind micro-structure computed by the clumpy wind models considered in this paper is not sufficient, per se, to retrieve the time-variability levels observed in classical sgxb such as Vela X-1, 4U 1700-37, OAO 1657-415 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 X-ray ionizing feedback in Manousakis & Walter, 2015). However, as discussed above, while the LDI simulations used here provide a perfect test-bed 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 line-driven clumpy wind outflow, future work will thus use LDI simulations more specifically tuned to Vela X-1 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łodowska-Curie 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 X-ray 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/2015-014 and by the Interuniversity Attraction Poles Programme by the Belgian Science Policy Office (IAP P7/08 CHARM). The simulations were conducted on the Tier-1 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., Simon-Diaz, 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., Christensen-Dalsgaard, 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.
• Cruz-Osorio et al. (2017) Cruz-Osorio, A., Sanchez-Salcedo, F. J. & Lora-Clavijo, 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.
• Gimenez-Garcia et al. (2016) Gimenez-Garcia, A., Shenar, T., Torrejon, J. M., Oskinova, L., Martinez-Nunez, S., Hamann, W.-R., Rodes-Roca, J. J., Gonzalez-Galan, A., Alonso-Santiago, J., Gonzalez-Fernandez, 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., Martinez-Nunez, 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ínez-Núñez et al. (2017) Martínez-Núñez, S., Kretschmar, P., Bozzo, E., Oskinova, L. M., Puls, J., Sidoli, L., Sundqvist, J. O., Blay, P., Falanga, M., Fürst, F., Gímenez-García, Á., Kreykenbohm, I., Kühnel, M., Sander, A., Torrejón, J. M. & Wilms, J. (2017), Sp. Sci Rev .
• Martins et al. (2016) Martins, F., Simon-Diaz, S., Barba, R. H., Gamen, R. C., Ekstroem, S., Simón-Dí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 0035-8711) 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. 2820-2831. 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(3-4), 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. 921-934. 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., Rodes-Roca, 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 :

 Si,j,k=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯si,j⋅si,k (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 : cross-shaped 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 micro-structure 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 Multi-scale 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 micro-structure at any time, in particular in the regions susceptible to be accreted.

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 Courant-Lewy-Friedrich 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 Yin-Yang 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 MPI-AMRVAC (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 off-centered 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 off-centered inhomogeneities.

### b.2 Code and solving scheme

We solve the adiabatic equations of hydrodynamics under their conservative form using the finite volume code MPI-AMRVAC (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 shock-capturing 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 user-defined type of boundary condition to work with time-varying boundary conditions stored in data files. A test-case of time-varying boundary conditions is now part of the main branch of the code.

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters