Stellar wind in state transitions of highmass Xray binaries
^{1}
Astronomical Institute, Academy of Sciences,
Boční II 1401/1, 141 00
Praha 4, Czech
Republic
email: cechura@astro.cas.cz; had@sunstel.asu.cas.cz
^{2}
Faculty of Mathematics and Physics, Charles
University, Prague,
Czech Republic
Received: 19 July 2014
Accepted: 30 November 2014
Aims. We have developed a new code for the threedimensional timedependent raditation hydrodynamic simulation of the stellar wind in interacting binaries to improve models of accretion in highmass Xray binaries and to quantitatively clarify the observed variability of these objects. We used the code to test the influence of various parameters on the structure and properties of circumstellar matter.
Methods. Our code takes into account acceleration of the wind due to the Roche effective potential, Coriolis force, gas pressure, and (CAK) radiative pressure in the lines and continuum of the supergiant radiation field that is modulated by its gravity darkening and by the photoionization caused by Xray radiation from the compact companion. The parameters of Cygnus X1 were used to test the properties of our model.
Results. Both two and threedimensional numerical simulations show that the Coriolis force substantially influences the mass loss and consequently the accretion rate onto the compact companion. The gravitational field of the compact companion focuses the stellar wind, which leads to the formation of a curved conelike gaseous tail behind the companion. The changes of Xray photoionization of the wind material during Xray spectralstate transitions significantly influence the wind structure and offer an explanation of the variability of Cygnus X1 in optical observations (the Hα emission).
Key words: accretion, accretion disks / hydrodynamics / methods: numerical / stars: winds, outflows / Xrays: binaries
© ESO, 2015
1. Introduction
Highmass Xray binary systems (HMXBs) are interacting binaries in which a compact companion, either a neutron star or a black hole, orbits a massive earlytype star, typically an OB supergiant. This type of stars is characterized by an enhanced massloss rate of the order of ~10^{6} M_{⊙} yr^{1}. The compact companion is immersed in the stellar wind and accretes material from it, giving rise to a strong Xray flux. The interaction between the companion and the stellar wind was a subject of many numerical studies, for example, Blondin et al. (1990), Hadrava & Čechura (2012), or Manousakis et al. (2012), which revealed that the wind of the massive star is severely disrupted by the gravity and photoionization of the companion. In earlytype stars, the stellar wind is predominantly driven by the line absorption of radiation from the primary by the wind material. This mechanism can impart sufficient momentum to accelerate the material to velocities of ~1500 km s^{1}.
In this study, we present an enhanced version of our hydrodynamic code introduced in Hadrava & Čechura (2012), which we have used to investigate the properties and dynamics of the stellar wind in Cygnus X1. First, we investigate in twodimensional simulations the role of various physical parameters that influence the interplay between the stellar wind and the compact companion: the mass ratio of the binary components, and parameters of the linedriven wind model. Then, we show the results of threedimensional simulations, revealing the importance of Xray photoionization. In Sect. 2, we describe in detail the physical model we used in our simulations and provide a summary of all physical effects and phenomena involved in the model. We specify the numerical hydrodynamic code in Sect. 3. The revised results of the anisotropic stellar wind in HMXBs can be found in Sect. 4 together with the description of the multitude of simulations computed for specific values of selected physical parameters. All calculations show the formation of a gaseous tail with complicated dynamics behind the compact companion. This tail is a product of a process that resembles BondiHoyleLyttleton accretion (Bondi & Hoyle 1944). We discuss the results and draw our conclusions in Sect. 5.
2. Physical model of the stellar wind
The most commonly used model of stellar winds from massive earlytype stars has been introduced by Castor et al. (1975, hereafter CAK). This model describes the mass loss driven by line absorption and scattering of the supergiant radiation field. It is based on a simple parametrization of the line force (by parameters α and k), which represents the contribution of spectral lines to the radiative acceleration by a powerlaw distribution function.
In a rapidly expanding stellar wind, where the radial velocity gradient is assumed to be large, the line optical depth τ_{L} in the radial direction can be reduced to a purely local quantity (Castor 1974). Following CAK, (1)where ρ and v are the density and the velocity of the wind material. v_{th} is the thermal velocity of the ion, which is generally different for each element, and κ_{L} is the monochromatic line opacity per unit mass of the particular transition of the ion. It is inversely proportional to the broadening of the line by v_{th}. The contribution of each line to the overall radiative drag is proportional to the oscillator strength, that is, to the frequencyintegrated κ_{L} and hence to the product κ_{L}v_{th}. Two new variables are introduced in CAK: the local optical depth parameter t, which is independent of the line strength and is given in the case of an expanding atmosphere by (2)where the quantity σ_{e} is the electronscattering coefficient and v_{th} is the reference thermal speed of an ion, usually of the hydrogen atom at the temperature 10^{4} K. The second variable is the ratio of the line opacity coefficient κ_{L} to the electronscattering coefficient σ_{e}(3)Then, as a consequence, (4)Following Castor (1974), under the assumption of the Sobolev approximation, the total force due to lines can be approximated as (5)where L_{∗} is the luminosity of the primary star, and (6)is the force multiplier function of the local optical depth parameter t, which is a convenient means of parametrizing the line force that is often used in wind calculations. α represents the slope and k the amplitude of M(t), at t = 1 (i.e., k = M(1)).
Abbott (1982) improved the CAK theory by calculating the line force considering the contribution of the strengths of the hundreds of thousands of lines. He also included a third parameter (δ) that takes into account the changes in ionization throughout the wind. Despite this immense effort to give a more realistic representation of the line force, evident discrepancies with observational data still remained.
To elucidate the grounds of the CAK theory from the first principles of physics and to avoid a dependence of the parameter k on an arbitrarily chosen value v_{th}, Gayley (1995) introduced an alternative notation in which k is replaced by a dimensionless linestrength parameter , the value of which (only slightly varying around 10^{3}) can be found from the metalicity of the wind. We refer to Gayley for a detailed explanation, but because his alternative approach is equivalent to the standard CAK theory, we adhere here to the standard CAK notation.
2.1. mCAK model and the finitedisk correction factor
In the CAK formalism, the authors adopted the pointlike source approximation of a stellar radiation field. This assumption is, however, a poor one close to the stellar photosphere, where the wind is rapidly accelerated. Later improvements to the theory reported by Friend & Abbott (1986) and Pauldrach et al. (1986) extended the CAK formalism by adding the effect of an outward centrifugal acceleration to onedimensional models of the wind outflow in the equatorial plane. Both papers independently derived a modified CAK model (mCAK) that relaxes the CAK simplifying approximation of a pointlike star and properly accounts for the finite cone angle subtended by the star. Assuming a uniformly bright spherical source of radiation, they introduced a multiplicative finitedisk correction factor K_{FDCF} (commonly referred to as the FDCF) to the force multiplier M(t). The FDCF is attained by adopting the exact optical depth rather than the radial one, (7)where for their stellar radius R_{∗}, and σ is given by Castor (1974)(8)
2.2. Photoionization
By its nature, the dynamics of the wind can be strongly influenced by the ionization structure of the medium. The heating and photoionization by the Xray flux depopulate the electron levels available to absorb the momentum of radiation from the primary and decrease the radiative drag on the wind.
Assuming an optically thin gas in local ionization and thermal balance, irradiated by a pointlike source of Xrays of a given spectral shape, the ionization structure of the medium is determined solely by the ionization parameter ξ (Tarter et al. 1969) given as (9)where L_{x} is the Xray luminosity of the source, n is the nucleon number density of the gas, and r_{x} is the distance from the Xray source. While in reality, optical depth effects probably play a role in the wind dynamics in HMXBs, the difficulty of realistically calculating the radiative transfer is too challenging for the scope of this paper. For the purposes of our calculations, the gas was therefore assumed to be optically thin.
The preliminary results of our model support the hypothesis that the Xray ionization tends to slow down the wind material in the immediate vicinity of the compact companion and thus to increase the overall accretion rate (Hatchett & McCray 1977). However, if the zone of full ionization extends to the proximity of the surface of the donor where the wind does not yet reach the escape velocity, the outflow is obstructed immediately at the base of the wind. Therefore, the Xray feedback effectively cuts out the accretion process from additional material. The dependence of the line force on ξ is very complicated owing to the many different ions responsible for the variation. The value of ξ at which the cutoff occurs depends on the nature of process responsible for the acceleration of the wind. For a conservative estimate for a typical abundance of the absorbing element relative to hydrogen of 10^{4} we adopt a cutoff value of log ξ ≃ 2 for C 4 and O 6 (Hatchett et al. 1976).
2.3. Parameterizing the force multiplier
To capture the flattening of M(t) for small t, we followed Owocki et al. (1988), and modified Eq. (6) in such a manner that the force multiplier becomes constant for small t (as it must be in the optically thin limit), (10)where τ_{max} = tη_{max}(ξ) and η_{max} is a cutoff to the maximum line strength. Equation (10) now explicitly indicates that M depends on both t and ξ. Consequently, we describe the influence of ξ on the stellar wind dynamics via the CAK parameters k(ξ) and η(ξ). Allowing for ξdependence in these quantities captures two systematic changes in M(t,ξ) with ξ: decreasing η_{max}(ξ) with increasing ξ allows the turnover in M(t) to shift to higher t with increasing ξ, while a reduction in k(ξ) at large ξ describes the overall decrease in M for larger ionization parameters.
We followed Stevens & Kallman (1990) in choosing that α does not explicitly depend on ξ. Although it is possible to allow α to vary, by doing that, we only add unwarranted complexity to the parametrization. We chose a value of α that allowed us to reproduce the correct observed terminal velocity for a single star of the appropriate spectral type in the standard CAK theory. Stevens & Kallman (1990) fitted M(t) with two relevant parameters k and η_{max} as a function of ξ. They found that in the regime of , the fitting parameters k and η_{max} can in turn be fitted with the following exponential functions of ξ alone: (11)and (12)In this way, the complicated behaviour of M(t) with ξ can be approximated by an analytic function of ξ. Using the analytic formulae for k and η_{max} slightly reduces the accuracy of the representation of M(t), but for the typical wind model we employed in our simulations, the accuracy of the fits is still within a factor of 2 for the appropriate values of ξ and t.
Taking into account Eqs. (2) and (10), we finally derive to an expression for the line force in the framework of the mCAK model with the inherent dependency on the ξ parameter (13)
2.4. Limb and gravitydarkening
Fig. 1 Ratio of the tangential to radial component of the radiative force in the plane of stellar centres and rotational axis. Mass ratio q = 2 / 3, volume of the Roche equipotential (drawn as a light blue line) 0.95 of the critical Roche lobe. Limb and gravitydarkening (0.0, 0.0) and (0.6, 0.25) are plotted in the left and righthand panel. The black isocontours correspond to the values of − 3, − 2, and − 1. 

Open with DEXTER 
In the mCAK model, the radiation force is often calculated assuming a uniformly bright finitesized spherical star. As Cranmer & Owocki (1995) and Pelupessy et al. (2000) showed for the rapidly rotating stars or Hadrava & Čechura (2012) for the tidally distorted surface of the primary, the changes of the shape of the stellar surface induce the gravitydarkening effect (von Zeipel 1924) as a function of latitude or of both latitude and longitude. For either a rotating or nonrotating star the decrease of the temperature outwards in the photosphere produces a limbdarkening effect that also modifies the finitedisk correction factor. The theoretical formalism for computing a selfconsistent radiation force for nonspherical rotating stars, including the effects of stellar oblateness and limb and gravitydarkening, was developed by Cranmer & Owocki (1995). However, to distinguish the effects of each one of these competing processes upon the wind structure, these authors presented a semiquantitative analysis and estimated that the limbdarkening effect increased the massloss rate (Ṁ) by about 11% to 13% over the uniformly bright models. However, this higher mass loss would imply a reduction in the wind terminal speed.
According to von Zeipel’s theorem (von Zeipel 1924), the distribution of the radiative flux F and the effective temperature T_{eff} across the tidally or rotationally distorted surface of a star should be given by the local gravity acceleration g as (14)where β ≡ 0.25 in the standard von Zeipel formula. Because the radiative force f_{rad} given by Eq. (13) is proportional to the flux, it enhances the wind more in the polar region of the star where g reaches its highest value than it does at the equator and especially in the line joining the component stars, where it has its lowest value. For the parameter values chosen in our calculations, g is higher by 23% at the poles and lower by 15% and 40% in the directions towards the Lagrangian points L_{2} and L_{1}, respectively, than its value in the perpendicular direction. This effect is the opposite of the direct modulation of the wind by the effective potential. In the region where the distance from the star is similar to its radius, this effect is additionally enhanced by the ellipticity of the star. In evaluating of the radial component of the radiation line force f_{rad} in Eq. (13), we therefore used a directionally dependent radiative flux F that is modulated by the gravity and limbdarkening effects.
The assumption of purely radial radiative drag given by Eq. (13) is accepted here in agreement with the mCAK models of spherically symmetric stellar winds. This approximation neglects the tidal distortion of the binary component and the consequent anisotropy of its radiation. To estimate an error caused by this simplification, we computed an auxiliary model of a radiation field in the vicinity of a precise Roche equipotential with a prescribed limb and gravitydarkening. In each node of a Cartesian coordinate mesh around the donor star, we used this model to integrate the momenta of photons radiated from all visible parts of the stellar surface. The ratio of the tangential to radial component of the radiative force projected with respect to the radiusvector from the centre of the star is depicted in Fig. 1 in the x − z plane (spanned by directions towards the companion and the rotational axis from the centre of the primary). This ratio is highest at about 23% closely above the Roche equipotential between the pole and the L_{1} point where the stellar surface is appreciably skewed by the tidal distortion. The radiative drag pushes the stellar wind here toward a higher equipotential rather than radially from the centre, but the nonradial component rapidly decreases below 10% (drawn by the innermost isocontour in Fig. 1). The nonradial component is partly diminished by the gravitydarkening in the space between the components because it dims the tidal bulge toward the L_{1} point. The values of the limb (u = 0.6) and gravitydarkening (β = 0.25) used in this calculation are, however, only a rough guess based on simplifying assumptions (as discussed in our previous paper, Hadrava & Čechura 2012) and should be replaced by frequencyweighted means corresponding to the lines causing the radiative drag. Moreover, this model does not take into account the radiation reprocessed in the wind. We therefore did not include it into the following hydrodynamic calculations.
3. Radiation hydrodynamic simulations
3.1. Numerical hydrodynamics
The results presented here were obtained using the secondgeneration radiationhydrodynamic model of the stellar wind in HMXBs (Hadrava & Čechura 2012). Simulations were conducted on a threedimensional Cartesian coordinate grid corotating with the modelled components of a binary. The model uses the timedependent equations of Eulerian hydrodynamics. The relevant equations for mass, momentum, and energy conservation are
respectively, where (18)is the total energy density and e is the specific internal energy. We adopted an ideal gas equation of state, P = (γ − 1)ρe, where γ = 5 / 3 is the ratio of specific heats. The model incorporates the physical effects discussed in Sect. 2. The gravity of the primary and compact component reduced by the radiative pressure gradient in continuum (Thomson scattering) as well as centrifugal force due to the orbital motion are represented by the first term on the righthand side of Eq. (16) as a gradient of the effective potential Φ_{eff}. The second term, corresponding to the Coriolis effect, accounts for the noninertial nature of the assumed Cartesian reference frame, which corotates with the components of the binary with angular velocity ω. The last term represents the radiation line force exerted on a unit volume of the medium. f_{rad} is a radial vector with magnitude f_{rad}. The term in Eq. (17) is the net heating/cooling rate, which is properly defined in Sect. 3.3.
For all threedimensional simulations, we adopted an equidistantly spaced grid of 207 × 157 × 157 cells in x, y and z direction, respectively. Similarly, for the twodimensional case, we used a resolution of 407 × 307 cells in x and y direction.
In each integration step we computed a timestep Δt that satisfies the CourantFriedrichsLewy condition (19)where Δx, Δy, and Δz are the distances between the neighbouring grid nodes in the x, y, and zdirection, respectively, and c_{s} represents the local isothermal speed of sound. In each timestep, all physical variables evolve in accordance with the conservation equations. The employed numerical technique is based on an explicit Eulerian version of the piecewise parabolic method developed by Colella & Woodward (1984). To optimize the calculations, we attempted to find the highest value of C_{0} that allows for the longest possible timestep but still prevents the simulation from breaking down. We successfully ran our model for C_{0} = 0.7 without any apparent effect on the resultant condition. Another simulation for C_{0} = 0.9, however, broke down shortly after its initialization. To find a compromise between performance and stability, we settled down for a more conservative value of 0.5 which we adopted in all simulations presented in this paper.
3.2. Boundary conditions
The grid was initialized with a smooth, steady wind outflow coming from a specific Roche equipotential that represents the surface of the primary. Density was kept uniform across the whole inner boundary condition and constant at the value of ρ_{0}. We estimated ρ_{0} by assuming an initial massloss rate of Ṁ = 2 × 10^{6} M_{⊙} yr^{1}, which is consistent with value obtained from our early radial stellar wind model (Hadrava & Čechura 2012) and with setting the radial outflow velocity from the surface to the isothermal speed of sound. The inner boundary condition was specified as inflow. The velocity of the material that is incoming from the inner boundary thus always points inward of the computationally active area. Its tangential component was set to 0. Its magnitude, however, was allowed to vary so that the amount of inflowing material is dependent on the conditions in the computational active area. When necessary, the magnitude of the velocity vector can drop to zero to accommodate the case of zero outflow from the donor. In this case, the wind material may even reverse to fall down onto the surface of the donor. For the outer boundaries, we adopted an outflow boundary condition by setting gradients of all quantities to 0, permitting the matter to freely flow out of the grid. Additionally, the velocity at the outer boundary was prevented from pointing inward of the computationally active area.
As a boundary condition for the compact companion, we used the sinkparticle treatment. A sinkparticle is a region that accreted incoming material but has no internal structure. It is defined by the accretion radius r_{a}. Any matter within the accretion radius of a sinkparticle was removed from the computational grid, and its mass was added to the mass of the sinkparticle. If a computational cell was located only partially within the accretion radius or if the accretion radius was smaller than a cell that a sinkparticle resides in, then only a proportional amount of matter was accreted onto the sinkparticle. The position of the sinkparticle was not arbitrarily fixed within the computational grid, but since we neglected the amount of angular and linear momentum that is being transferred from the accreted gas, the sinkparticle, representing the compact companion, remained stationary for the whole duration of the simulation.
3.3. Radiative cooling and Xray heating
The radiative cooling and Xray heating were computed under the assumptions of optically thin gas illuminated by an isotropic photon flux with a bremsstrahlung spectrum of 10 keV (cf. Kallman & McCray 1982, for details of the photoionization calculation). We approximated these rates with an analytic expression as a function of the local gas density, temperature, and ionization parameter – cf. Blondin (1994). The net heating/cooling rate for our approximate formula is given by , where Here T_{x} is the characteristic temperature of the bremsstrahlung radiation. The energy sources were calculated in a separate step from the hydrodynamics and were computed under the assumption that the gas is in ionization equilibrium with a constant isotropic Xray photon flux from the compact companion. The recombination timescale was assumed to be sufficiently fast for the ionization state of the gas to be given solely by the local density, temperature, and photon flux (Fransson & Fabian 1980). At each point in the wind, the heating rates Γ_{x} due to Xrays and Γ_{Comp} due to Compton heating and the cooling rate Λ due to locally emitted radiation were calculated as a function of the local temperature and ionization parameter ξ. The effects of the primary were accounted for by preventing the wind temperature from dropping below the photospheric temperature of the primary T_{eff}, which simulates the UV heating. Although the temperature of the wind is probably slightly lower because of adiabatic cooling, the results of our simulations depend only very weakly on this temperature.
4. Results for Cygnus X1
Simulation parameters.
Fig. 2 Series of three twodimensional simulations of the stellar wind in the equatorial plane of Cygnus X1 with an increasing α parameter. These simulations disregard the effects of the Xray ionization on the dynamics of the wind. Additionally, the k parameter is kept constant at 0.25. A clear trend of a shrinking gaseous tail and decreasing matter captured in the vicinity of the compact companion with increasing α is evident. The three upper panels shows the density distribution after the solution reached the quasistationary state, ~3P_{orb}. The lower panels represent corresponding velocity magnitude plots. (This figure is available in colour in the electronic version of the paper.) 

Open with DEXTER 
The parameters of the following simulations were chosen in accordance with the observed parameters of Cygnus X1/HD 226868 (see Table 1) – a representative of the HMXBs and wellknown black hole candidate. The compact component of the system is a stellarmass black hole orbiting a massive O9.7Iab supergiant. In the past four decades, many estimates have been made of the masses of the binary components – CaballeroNieves et al. (2009) for instance inferred wide observational limits stating the mass of and for the primary and the compact component of the system, respectively. However, most of these estimates are unreliable because they are based on anunsatisfactory determination of the distance to Cygnus X1. Recently, Orosz et al. (2011) were able to constrain some of the principal parameters of Cygnus X1. By using an unprecedentedly precise distance from a trigonometric parallax measurement for Cygnus X1 (Reid et al. 2011), which is , they found masses of M_{∗} = 19.16 ± 1.90 M_{⊙} and M_{x} = 14.81 ± 0.98 M_{⊙} for the Ostar and black hole, respectively. These estimates are considerably more direct and robust than the previous ones, owing largely to the new parallax distance. However, Ziółkowski (2014) recently showed that the mass of the supergiant is inconsistent with the evolutionary models for the massivecorehydrogen burning stars. Based on the evolutionary models, the mass of the supergiant is most likely in the range of 25 to 35 M_{⊙}. The corresponding mass of the black hole is in the range of 13 to 23 M_{⊙}. If, as a result of rotationinduced mixing, the hydrogen content of HD 226868 is equal to about 0.6 (as suggested by some observational data), then its present mass may be somewhat lower, ~24 M_{⊙}.
For the purpose of our model, we adopted values of M_{∗} = 24 M_{⊙} and M_{x} = 16 M_{⊙}. The orbital period is 5.6 days and the corresponding orbital separation of the components is 45.4 R_{⊙}. The critical Roche lobe of this system has a mean radius of 18.9 R_{⊙}: the Lagrangian point L_{1} is at a distance of 24.5 R_{⊙}, while the intersection of the critical equipotential with the rotational axis of the supergiant is 17.7 R_{⊙} from the supergiant centre. We set the surface of the donor star to a mean radius of 18.5 R_{⊙} (i.e. at the equipotential of size 22.3 R_{⊙} towards the companion and 17.5 R_{⊙} towards the pole). The gravity acceleration g varies across this equipotential from 1.5 to 3.1 × 10^{3} cm s^{2} towards the L_{1} and the pole, respectively, that is, the mean value log g = 3.3 is close to the observational limit 3.00 ± 0.25 found by CaballeroNieves et al. (2009). We assumed that the primary is tidally locked with the companion and exhibits synchronous rotation.
The interplay between a stellar wind and the gravitational field of a compact object in HMXBs can be complicated by the interaction of several competing processes. To isolate and understand the importance of each of the physical effects that influence the gas flow, we present a series of simulations in which only a specific parameter is being changed. In the following subsections we show results of restricted simulations we ran to discuss the effects of α parameter, mass ratio of the compact companion to the primary, and the Xray ionization on the overall solution. Because of the number of simulations and their high computational time requirements, we were forced to restrict our model to a twodimensional grid. But here we are only interested in qualitative changes that take place in the gas flow with varying parameters. The general conclusions we make can be easily extended to the threedimensional case. Later, we examine the full threedimensional case that incorporates all of the physics available in our numerical model.
In the twodimensional simulations, we used an equidistantly spaced grid of 407 × 307 cells in x and y direction. The computational volume has a range of x = [ −1,1.66 ] × D and y = [ −1,1 ] × D, where , yielding a spatial resolution of the computational grid dl = 2.1 × 10^{8} m. The materialdonating primary star is centred at [ x,y ] = [ 0,0 ], while the position of the compact companion is [ x,y ] = [ D,0 ]. The timestep Δt is adjusted every computational step, satisfying the condition given by Eq. (19) and typically reaches ~10^{4}P_{orb}. Most of our simulations were run for the total duration of up to three orbital periods. In all cases, the simulations quickly converged on a quasistationary solution, taking less than one orbital period to establish a dynamical balance.
4.1. Dependency on α parameter
Our first twodimensional simulation was of the stellar wind from the corotating primary in which the only influence of the compact companion on the gas flow is via its gravity. The Xray luminosity of the source was set to 0 and the k parameter was kept constant at 0.25. The parameter α took progressively higher values of 0.5, 0.6, and 0.7. Other parameters of the simulations were taken from Table 1. Because of the lack of the effects of the Xray ionization, this simulation is expected to be closest to the approximation of uniform, axisymmetric accretion – BondiHoyleLyttleton (BHL) accretion.
The density distribution and the velocity magnitude of the wind in the orbital plane, after the quasistationary solution was reached (typically ~3P_{orb}), are displayed in Fig. 2. The red and blue disks in the upper and bottom panels represent the surface of the supergiant with a preset value of density ρ_{0}. Small tangential variations in density are caused by the imperfection of the inner boundary condition, more specifically, by the fact that we tried to capture a round object within the Cartesian coordinate grid. In principle, these artefacts may be avoided by adopting a computational grid that follows the symmetry of the problem more closely – e.g. a spherical coordinate grid. But our intention is to investigate the complex interaction between the stellar wind and compact companion. As we saw so far, these interactions may involve a dense matter stream between the components of the binary, strong shocks in the wind medium, and extensive accretion structures around the compact companion. In general, these structures do not have to be spherically symmetric at all, which justifies the need for a more universal computational grid. Nonetheless, the variations are small and do not influence the overall picture. We can use them, on the other hand, as streamline markers. By doing so, we reveal spiralling trajectories of the wind out of the system – clearly the effect of the Coriolis force.
For the ideal gas with the adiabatic index γ = 5 / 3, capturing of the stellar wind in the vicinity of the compact companion leads to the formation of an accretion disk. The size of the disk is predominantly determined by the amount of material captured within the gravity well of the compact companion, which is inversely related to the velocity magnitude of the wind passing in close proximity of the compact object. When the supersonic wind meets an obstacle in form of the accretion disk, it creates a standing bow shock extending to a taillike structure of hot turbulent slowly moving gas. As the bow shock extends downstream, it is affected by the Coriolis force, which causes the flow to curve clockwise. However, the shape of the bow shock remains roughly axisymmetric.
We note that increasing values of α lead to higher velocities of the wind passing near the compact companion. This effect has two consequences. Firstly, the size of the accretion disk shrinks as much less material is captured. The obstacle profile is smaller, therefore the bow shock is less pronounced. Secondly, the angle of the bow shock “wings” decreases, leaving the shock cone narrower. The higher values of α also cause the gas within the tail to be less turbulent. This is a consequence of several effects that the higher velocity profile has on the gas dynamics. The occasional oscillations seen in the first two columns of Fig. 2 are a result of the nonzero angular momentum in the accreting stellar wind with respect to the compact object. The exact axisymmetry of the accretion flow is violated by the orbital velocity of the system, even for a corotating companion. Because the wind approaches the compact object at a finite angle with respect to the line of centres, wind on the trailing side of the compact object (lagging behind in orbit) is closer to the primary and therefore denser than the wind on the leading side due to the 1 /r^{2} divergence of the wind. Initially, the asymmetry in the flow leads to a higher ram pressure on the highdensity side that drives the bow shock off towards the lowdensity side. As the bow swings over to this side, it begins to intercept more of the wind from this side than from the highdensity side and is eventually pushed back to the highdensity side, where the whole process repeats itself in a quasiperiodic fashion. The higher value of α tends to homogenize the outflow of material in the direction of the L_{1} point. For α = 0.5, we see origins of the formation of the focused stellar wind, whilst for the higher values of α the wind becomes more isotropic.
It is also instructive to view these results in the perspective of Gayley’s notation of the CAK process, as mentioned in the beginning of Sect. 2. Substituting T_{eff} = 30 000 K and into formula (59) by Gayley (1995), which yields , we find k ≃ 0.54, 0.13 and 0.034 for α = 0.5, 0.6 and 0.7, respectively. Vice versa, if the value of k is kept fixed for different α like in the above presented computations, we find from the inverse transformation the value of varying in a wide range, namely , 5. × 10^{3} and 780 × 10^{3}. For fixed k, a change of the lineopacity distribution given by α thus correlates with a variation of the overall opacity in lines characterized by .
Fig. 3 Integrated column density for the simulation shown in Fig. 2 as a function of orbital phase Φ timeaveraged over one orbital period. Each line represents the column density profile for different values of α – 0.5 (dashed line), 0.6 (dashdotted line), and 0.7 (solid line). 

Open with DEXTER 
To illustrate the density distribution in the system, we calculated the column density along the line of sight to the Xray source N_{H} in the orbital plane. Figure 3 shows the absorbing column density as a function of orbital phase Φ derived from the simulation results in Fig. 2 timeaveraged over one orbital period. All profiles, each corresponding to a different α parameter, exhibit some common characteristics. At early phases, the smooth wind component dominates, followed by a rise in N_{H} starting from orbital phase ~0.25. This increase represents the bow shock and the tail forming behind the accretion disk. However, the shape of the central part of the profile for various values of α differs significantly. For α equal to 0.6 and 0.7, we obtain a doublepeak structure, which is a clear sign of a welldefined bow shock that accumulates most of the gas. In the case of α = 0.5, much more gas is captured in the accretion disk and the distribution of matter within the tail is more isotropic. The column density profile is consequently more homogeneous, with a growing trend towards the peak value around orbital phase ~0.67. It declines from here to another local minimum at a phase of about ~0.78.
4.2. Effects of the companion mass
Fig. 4 Similarly as in Fig. 2, this series of six twodimensional simulations in the equatorial plane shows the dependency of the solution on the mass ratio M_{x}/M_{∗} and the mean radius of the primary. The masses of the primary and the compact companion for all three cases can be found in Table 2. The shape of the primary, represented by an equipotential surface, is a function of M_{∗} and M_{x}. The white index in the upperright corner of each panel represents the volume of the primary as a fraction of volume of the critical Roche equipotential. We disregard the effects of the Xray ionization in these simulations. α is set to 0.6 and the k parameter is kept constant at 0.25. Changes among the simulations become more pronounced as we take into account the scale differences. For a given value of P_{orb}, D grows with increasing overall mass of the system. The panels show the density distribution after the solution reached the quasistationary state, ~3P_{orb}. (This figure is available in colour in the electronic version of the paper.) 

Open with DEXTER 
As stated at the beginning of this section, there is no clear consensus on the masses of the components of Cygnus X1. We wish to investigate the effects that different masses and, consequently, ratios of M_{x}/M_{∗} might have on the wind structure and dynamics. The gravitational potential of the compact companion affects the wind massloss rate and its distribution across the surface of the primary. The massloss rate is enhanced along the line of centres of the binary system by two related effects. The first one is the tidal distortion of the primary, which lifts the base of the wind farther away from the centre of the primary in the direction of the compact companion. The second related effect is the weakening of the gravitational force along the line of centres so that the wind has to overcome less gravity in this direction. In this way, the massloss rate is strongly enhanced in the direction of the companion, giving rise to the highly anisotropically modulated stellar wind. This effect, however, is reduced by taking into account the gravitydarkening of the primary photosphere.
The value of the mean radius of the primary, derived from the observations, is even less well established than in the case of its mass. Thus, the need of investigating how the mean radius affects the wind structure and dynamics becomes apparent. The photosphere of the primary in our simulations is set to coincide with a selected surface of the effective equipotential with a given volume. The inner boundary condition is thus defined by specifying its volume as a fraction of the volume of the critical Roche lobe. We used two values of this parameter throughout all six simulations – 90% and 95% of the critical Roche equipotential.
Figure 4 presents a series of the twodimensional simulations showing the density distribution of the stellar wind in the orbital plane. We used three sets of masses of the primary and the compact companion (see Table 2) with the unique mass ratio indicated on the top of each column. The first column refers to one of the limiting cases of the range given by CaballeroNieves et al. (2009) with the corresponding mass ratio M_{x}/M_{∗} = 0.36. The second column represents the outcomes of the evolutionary models of the massivecorehydrogen burning stars by Ziółkowski (2014) – low hydrogen content case – with M_{x}/M_{∗} = 0.67. And the last column uses values found by Orosz et al. (2011) with M_{x}/M_{∗} = 0.77. In the upper panels, the volume of the primary corresponds to 90% of the volume of the critical Roche equipotential. In the bottom panels the value is 95%.
Fig. 5 Twodimensional simulations of the stellar wind in the equatorial plane of Cygnus X1 including the effects of the Xray ionization. Panels show the density distribution ρ, the velocity magnitude , and the ionization parameter ξ. α is set to 0.6, and k = k(ξ). The wind launching from the hemisphere facing the Xray source is noticeably slower than the wind coming from the opposite hemisphere, which lies in the Xray shadow (the dark conelike region in the right panels). A strong shock forms in the region where the fast and slow wind collide. The accretion of the slow wind is much more efficient, giving rise to an extensive disk that is enveloped by a strong shock front. In the region where log ξ ≥ 2, the linedriven acceleration mechanism of the material of the wind is switched off, which causes the wind coming from the hemisphere facing the Xray source to slow down. (This figure is available in colour in the electronic version of the paper.) 

Open with DEXTER 
Similarly to the previous simulations, we used the ideal gas with an adiabatic index γ = 5 / 3 and did not take into account the effects of the Xray ionization, ξ = 0. The value of k was kept constant at 0.25, and α was set to 0.6.
The results of the simulations in Fig. 4 were acquired with a varying spatial scale, with D as a unit of distance. We determined D by specifying the orbital period P_{orb} and the masses of the binary components. By increasing the overall mass of the system, we also increased the distance unit D in the simulations. D defines the physical dimension of the computational region, which was always set to (2.66 × 2)D in x and y direction. Therefore, while D changes, the relative distance between the binary components within the computational grid remains constant in all simulations. The size of the primary, on the other hand, is scaled up and down depending on the actual value of D. By changing the mass ratio M_{x}/M_{∗}, we also change the shape of the effective potential and, consequently, the shape of the surface of the primary. To comply with our initial condition settings, we need to adjust ρ_{0} with changing size of the primary to obtain the initial preset massloss rate Ṁ = 2 × 10^{6} M_{⊙} yr^{1}.
In Fig. 4, we conclude from the outcomes of our simulations that the effects of the arbitrarily chosen mean radius on the wind structure and dynamics are rather small. There are no qualitative changes in the solutions corresponding to the cases when the volume of the primary equals 90% or 95% of the critical Roche equipotential. Even quantitatively, the solutions are virtually identical. The smallscale fluctuation are caused by numerical instabilities and generally random turbulent motion of the gas. Therefore, we adopted 95% of the volume of the critical Roche equipotential as the volume of the primary in all following simulations.
The dependency of the solution on the mass ratio M_{x}/M_{∗} seems to be of more importance. By increasing the mass of the compact companion M_{x}, we also increase the amount of matter captured in the accretion disk, which results in the bow shock growing larger. The matter is accumulated within the gravity well of the compact companion and then released in a quasiperiodic fashion. The interaction of these released clumps with the rest of the shock gives rise to the turbulent environment of the tail.
4.3. Effects of the Xray ionization
Until now, we neglected all effects of the strong Xray source within the system on the wind from the primary. In the following twodimensional simulations, we expanded our physical model by identifying the source of the Xray luminosity with the position of the compact companion. We have used parameters from Table 1 amended by the value of α = 0.6, which makes the results comparable with the previous simulations in Sect. 4.1 (the case of α = 0.6), and in Sect. 4.2 (the case of M_{x}/M_{∗} = 0.67). The ionization parameter ξ is now, in contrast to the previous physical model, computed as a function of the local quantities – gas density, Xray luminosity, and the distance from the Xray source – as indicated in Eq. (9). The CAK parameters k and η no longer have fixed values. They are calculated as k = k(ξ) and η = η(ξ) according to Eqs. (11) and (12). We assume that the Xray radiation field is blocked by the surface of the supergiant. The primary therefore casts an Xray shadow, leaving the region behind it relatively uninfluenced by the Xray source.
In Fig. 5, we present the results of our simulations for two values of the the Xray luminosity. The upper panels which represent the case with the Xray luminosity equal to 1.9 × 10^{37} erg s^{1}, correspond to the low/hard Xray state of Cygnus X1. The lower panels, then, show the case of the Xray luminosity equal to 3.3 × 10^{37} erg s^{1}, which corresponds to the high/soft Xray state. Each column in Fig. 5 represents a different physical quantity. From the left, this is the density distribution ρ, the velocity magnitude  v , and the ionization parameter ξ of the stellar wind in the orbital plane of Cygnus X1.
It is evident, in comparison with the corresponding simulations in the previous subsections, that the Xray radiation field has a profound effect on the structure and dynamics of the wind. Ionizing the gas material severely limits the efficiency of the CAK linedriven mechanism and decreases the radiative drag on the wind. Even in the low/hard state, the outflow from the hemisphere facing the Xray source is noticeably slower than in the previous simulations. This is in direct contrast with the wind launching from the opposite hemisphere lying in the Xray shadow (depicted in the right panel of Fig. 5 as the dark conelike region) where the wind is relatively unaffected by the strong Xray source. The shadow cast by the primary provides an environment where the wind can be accelerated without interruption, leading to velocities that are an order of magnitude higher than in the slowwind region where the acceleration mechanism is impaired right at the base of the wind. Similarly to the simulation in Fig. 2, slowing down the wind caused the broadening of the bow shock that forms in front of the accretion disk around the compact object. The amount of the material captured in the accretion disk also increases.
The intensity of the stellar wind is highly anisotropic. In directions where the stellar wind coming from the facing hemisphere is additionally rarefied by the divergence of the streamlines, the ξ parameter reaches its highest values and effectively cuts off the outflow of the wind in these directions. The gas from these regions condensates into a narrow dense stream and passes in the vicinity of the L_{1} point in the direction of the compact companion. The massloss rate in the direction of the L_{1} point is thus significantly enhanced. This focused stellar wind, which resembles the Roche lobe overflow, interacts with the bow shock and the outer layers of the accretion disk. The region undergoes quasiperiodic changes in density. First, the material of the focused wind is accumulated by the bow shock. A part of it is accreted by the compact companion, but most of it sticks together in a form of a dense clump. When a critical density is reached, the whole clump is released and flows downstream out of the computational region.
Next to the interplay between the stellar wind and the compact companion, there is an additional source of hydrodynamic shocks in this simulation. This shock is a direct consequence of the Xray source in the system. The streamlines of the slow rarefied parts of the wind from the facing hemisphere are bent more by the Coriolis force and are brought to the region dominated by the fast less dense wind originating from the shadowed hemisphere. In the region that precedes the primary in orbit, both types of wind join, giving rise to a strong shock. There is no hydrodynamic shock in the region that lags behind the primary in orbit because the streamlines of the fast and slow wind diverge.
For the high/soft state (bottom panels in Fig. 5), the increase in Xray luminosity of the compact object effectively cuts off the outflow of the material from the entire facing hemisphere. The most prominent consequence of the increased Xray luminosity is the disruption of the narrow dense stream between the two binary components. The mass transfer from the primary to the accretion disk is therefore interrupted and the accretion disk shrinks as the matter within it is accreted onto the compact companion. Lacking the incoming flow of gas to compress the accretion disk, the gas from the disk is scattered in the wide region around the compact companion. Without any prevailing velocity field, the movement of the gas in this region is mostly slow and turbulent. The bow shock around the accretion disk completely disappears. The hydrodynamic shock preceding the primary in orbit transforms into a contact discontinuity. The same phenomenon appears in the region that lags behind the primary in orbit.
Fig. 6 Integrated column density for the simulation shown in Fig. 5 as a function of orbital phase Φ timeaveraged over one orbital period. The solid and dashed lines correspond to the high/soft and low/hard Xray state of Cygnus X1, respectively. 

Open with DEXTER 
Corresponding simulated timeaveraged column densities for both Xray states of Cygnus X1 are shown in Fig. 6. The profile of column density of the low/hard state is similar to that in Fig. 3 for α = 0.5. At early phases, the smooth wind component dominates, followed by a rise in N_{H} in orbital phase ~0.2. Similarly, the column density in the central part of the profile gradually grows, although the peak is not as sharp as in the case of α = 0.5. When we increase the Xray luminosity to simulate the transition to the high/soft state, the column density decreases as the gas is accreted onto the compact object or leaves its vicinity. The relative flatness of the central part of the profile suggests that, neglecting smallscale fluctuations in density, the gas around the compact object is distributed rather isotropically.
4.4. Threedimensional model
Fig. 7 Threedimensional simulation of the stellar wind of HD 226868 with an Xray luminosity equal to 1.9 × 10^{37} erg s^{1} corresponding to the low/hard Xray state of Cygnus X1. α is set to 0.6. η = η(ξ), and k = k(ξ) are functions of ξ. Columns show various cross sections throughout the computational volume. All mutually perpendicular planes are centred on the compact companion and are defined as x = D, y = 0, z = 0. The displayed quantities are the density distribution (the top panels), the velocity magnitude (the middle panels), and the ionization parameter ξ (bottom panels). (This figure is available in colour in the electronic version of the paper.) 

Open with DEXTER 
Fig. 8 Threedimensional simulation of the stellar wind of HD 226868 with an Xray luminosity equal to 3.3 × 10^{37} erg s^{1} corresponding to the high/soft Xray state of Cygnus X1. α is set to 0.6. η = η(ξ), and k = k(ξ) are functions of ξ. Columns show various cross sections throughout the computational volume. All mutually perpendicular planes are centred on the compact companion and are defined as x = D, y = 0, z = 0. The displayed quantities are the density distribution (the top panels), the velocity magnitude (the middle panels), and the ionization parameter ξ (bottom panels). (This figure is available in colour in the electronic version of the paper.) 

Open with DEXTER 
Fig. 9 Evolution of the density distribution in the orbital plane of Cygnus X1 showing gradual decay of the highdensity flow in the proximity of L_{1} point. The lefthand side panel corresponds to the quasistationary solution representing the low/hard Xray state of Cygnus X1. The righthand panel shows the outcome of the simulation when the new equilibrium is reached after adjusting the Xray luminosity of the compact object to the level corresponding to the high/soft state. The central panel represents a transient state roughly in the middle of the transition. The time index in the upper right corner of each panel is expressed in units of P_{orb}. (This figure is available in colour in the electronic version of the paper.) 

Open with DEXTER 
After developing and extensively testing our code in various twodimensional simulations, we now proceed with the more complex and extensive threedimensional simulation of Cygnus X1. So far, we have studied the roles of particular parameters and their influence on the shape of the solution separately. In this section, we present the results of the threedimensional simulations that encompass the entire physical model described in Sect. 2. We employ an equidistantly spaced grid of 207 × 157 × 157 computational cells in x, y, and z direction, respectively. Similarly to the simulations discussed in Sect. 4.3, we used parameters appropriate for Cygnus X1 from Table 1. We set α equal to 0.6. The ionization parameter ξ was calculated from Eq. (9) as a function of the local gas density, Xray luminosity and distance from the Xray source. The quantities k = k(ξ) and η = η(ξ) are functions of ξ according to Eqs. (11) and (12). The computational volume has a range of x = [ −1,1.66 ] × D, y = [ −1,1 ] × D, and z = [ −1,1 ] × D, where D = 3.16 × 10^{10} m, yielding a spatial resolution of the computational grid dl = 4.2 × 10^{8} m. The primary star is centred on [ x,y,z ] = [ 0,0,0 ], while the position of the compact companion is [ x,y,z ] = [ D,0,0 ]. The timestep Δt is adjusted every computational step, satisfying the condition given by Eq. (19), and typically reaches ~10^{4}P_{orb}.
In the first simulation, we set the Xray luminosity of the compact companion equal to 1.9 × 10^{37} erg s^{1}. This value corresponds to the low/hard state of Cygnus X1, making the results directly comparable with the first simulation in Fig. 5. Figure 7 shows various cross sections throughout the computational volume. All planes, defined as z = 0, y = 0, x = D, are mutually perpendicular and centred on the compact companion. The results of the threedimensional simulation resemble those acquired with the twodimensional model. The most prominent feature in the orbital plane is a dense gas stream in the direction of L_{1} point. The stream is brought to the proximity of the compact companion where it interferes with an extensive accretion disk enveloped by a bow shock. As the bow shock extends downstream, it creates a tail that has a similarly turbulent nature as the one we saw earlier in the twodimensional case in Fig. 5. We also have additional information about the vertical structure of the accretion disk and the tail. The accretion disk, especially its central region, is rather thick. This is likely caused by the insufficient grid resolution of the accreting region. One computational cell corresponds to over 4000 Schwarzschild radii, and several cells are needed to create a large enough pressure gradient to counter the gravitational force of the compact companion. The tail is mostly symmetric in the x − z and y − z plane with a fine turbulent structure.
We also note the second hydrodynamic shock in the region that precedes the primary in orbit where the slow wind hits the fast parts of the wind launched from the hemisphere in the Xray shadow. The shock is less pronounced than in the twodimensional case and also slightly shifted counterclockwise with respect to the projection of the surface of the primary in the orbital plane. Similarly, we note a spread of the region of the wind characterized by the fast solution. These effects are caused by the coarser computational grid we adopted in the threedimensional case where the distance between two neighbouring grid nodes is roughly twice the distance in the simulations in Sect. 4.3. As a consequence, the subsonic region that lies in the immediate proximity of the surface of the primary and where the material of the wind is strongly accelerated is only poorly resolved. This allows the wind to gain higher velocities before it leaves the Xray shadow cast by the primary. There is no hydrodynamic shock in the region lagging behind the primary in orbit because, similarly to the twodimensional case, the streamlines there diverge.
The density of the gas stream and, and consequently, the density in the inner parts of the accretion disk, is higher by about a factor of 2 than the twodimensional case. This increase is caused by the gravitational force that focuses the streamlines of the wind launched from the higher latitudes of the surface of the primary into the orbital plane and thus brings more material into the stream. There is also an overall increase of velocity of the wind in the region of the Xray shadow. This increase originates from the different geometrical formulation of the problem. In the twodimensional simulation, gas pouring into the computational region from the inner boundary condition is diluted as 1 /r^{2} as it advances towards the outer boundary. In the threedimensional case this dilution factor is 1 /r^{3}. The density of the wind therefore generally decreases more quickly with the distance from the primary in the threedimensional case. As a result, the linedriven force is enhanced since Eq. (13) is inversely related to the local density. This effect becomes stronger farther away from the surface of the supergiant in the region of the Xray shadow where it is not influenced by the ionization. The velocity of the outflowing gas in these regions is increased by approximately 10%. On the other hand, the lower density increases the ξ parameter. In the regions close enough to the Xray source, the drop of the k parameter with growing ξ in Eq. (13) completely nullifies and in some regions even surpasses the effect of lower density on the velocity structure of the wind. As in the two dimensional case, there are also numerical artefacts due to the staircasing of the inner boundary condition, which is more prominent in the threedimensional case because of the lower resolution used. This artefact is caused by the roughness of the inner boundary condition and grows stronger as the computational grid becomes less refined. The staircasing affects the smoothness of the wind solution and is responsible for the small tangential variations in density and velocity of the wind. Especially the streamlines launched from the transient region between areas in the Xray shadow and those that are completely exposed to the Xray source are strongly affected. These artefacts are visible in Figs. 7 to 9 as radial rays and are particularly pronounced at angles where the stellar surface is skewed with respect to the coordinate grid.
The simulation evolves for 2.5 P_{orb} and reaches a quasistationary state in ~1P_{orb}, after which we observe only quasiperiodic releases of material accumulated in the vicinity of the compact companion. At 2.5 P_{orb} after the beginning of the simulation, we abruptly increased the Xray luminosity of the compact object to 3.3 × 10^{37} erg s^{1} corresponding to the high/soft state of Cygnus X1. A new equilibrium is reached at around 0.25 P_{orb}. Since the stellar wind is highly supersonic for the most parts, this time roughly corresponds to the time needed for the gas launched from the facing hemisphere of the donor to leave the computational area. The results of the simulation after the new quasistationary solution is reached are shown in Fig. 8. The planes defining the cross sections of the computational volume are the same as in the previous case.
We note that the parametric change of Xray luminosity significantly influences gas dynamics in the vicinity of Cygnus X1. While the wind, launched from the parts of the surface of the primary that lie in the Xray shadow, is relatively unaffected, the wind originating from the facing hemisphere experiences dramatic changes because the CAK linedriven mechanism is seriously impaired. As the ξ parameter grows, the efficiency of the momentum transfer between the material of the wind and the radiation field declines. At this value of the Xray luminosity, the bubble of full ionization (approximated by the condition log ξ ≥ 2) almost reaches the surface of the supergiant and prevents the wind from achieving the escape velocity. The wind falls back and the outflow of the material in the direction of L_{1} point is interrupted. The disruption process of the dense stream of gas between the binary components is depicted in Fig. 9, where we show the density distribution in the orbital plane. The decline of the dense gas stream is a quick process, it takes only around 0.2 P_{orb} since the increase of the Xray intensity before it vanishes.
As the gas in the accretion disk is accreted and the transfer of the new material from the primary is obstructed, the accretion disk becomes less dense and shrinks. Eventually, the whole material entrapped in the gravity well of the compact companion would be accreted and the accretion disk would cease to exist. This would, however, require a considerably longer duration of the simulation, which is beyond our computational capabilities.
4.5. 3D visualization
For illustrative purposes, we produced an interactive figure (Fig. 10) containing an isodensity surface for ρ = 2.5 × 10^{9} kg m^{3}. The isodensity surface corresponds to the low/hard Xray state of Cygnus X1 from the simulation presented in Fig. 7. This approach allows us to visualize the global distribution of the accretion disk and the hydrodynamic shocks.
Fig. 10 Interactive threedimensional isodensity surface for ρ = 2.5 × 10^{9} kg m^{3}. 

Open with DEXTER 
5. Discussion and conclusions
We have presented our enhanced radiation hydrodynamic model of the stellar wind in HMXBs, which we used to simulate the circumstellar environment of Cygnus X1. First, in a series of twodimensional simulations, we investigated the role of various parameters on the distribution and dynamics of the stellar wind. Then we made use of the capabilities of our code to perform a more complex threedimensional analysis of the problem. We find that the wind parameters and its ionization structure established by the presence of a strong Xray source, have a significant impact on the structure and dynamics of the stellar wind.
Both the two and threedimensional simulations show the formation of an extensive bow shock enveloping the accretion disk and the compact companion. This conelike structure is curved by the Coriolis force as it advances beyond of the computational area, but it remains roughly axisymmetric and resembles BHL accretion. The position and orientation of the shock depend on the relative velocity of the compact object through the stellar wind, mass of the compact companion, and the Xray luminosity. On the other hand, the shock is relatively insensitive to the changes in the mean radius of the primary.
The outcomes of the twodimensional simulations support the importance of the Xray feedback in HMXBs. The Xray ionization tends to slow down the wind material in the immediate vicinity of the compact companion and thus to increase the overall accretion rate (Hatchett & McCray 1977). However, if the zone of full ionization extends to the proximity of the donor surface where the wind does not yet reach the escape velocity, the outflow can be obstructed immediately at the base of the wind, which effectively leaves the accretion process cut out from the additional material.
We investigated the properties of the stellar wind in Cygnus X1 in low/hard and high/soft Xray state and simulated the transition between these two states. The threedimensional simulation revealed dramatic changes of the material outflow from the primary as a response to the increased Xray luminosity. The dense gas stream in the proximity of the L_{1} point, which resembles the Roche lobe overflow in semidetached binaries, completely subsided. The increase of the Xray radiation from the compact companion ionized the base of the wind in the stellar photosphere, rendering the linedriven mechanism incapable of accelerating the wind. The material cannot reach the escape velocity and falls back onto the surface of the primary. The duration of the decline of the gas stream (~0.25 P_{orb}) is relatively short and corresponds to the time needed for the gas launching from the facing hemisphere of the primary to leave the accreting region.
The suppression of the dense gas stream as the system switches from low/hard state to high/soft state fully agrees with earlier optical observations of Cygnus X1 (Hadrava & Čechura 2013), where a decomposed spectral line component of Hα corresponding to the circumstellar matter anticorrelates with the soft Xray emission.
To minimize the computation cost, we made several simplifying approximations in our threedimensional numerical model. Probably the most drastic one was neglecting the optical depth effects in the optical and Xray region which would have required solving the radiative transfer in the stellar wind simultaneously with the hydrodynamic calculations. These effects probably play a role in the wind dynamics in HMXBs, but realistically calculating the radiative transfer is beyond the scope of this paper. The situation is less severe in the case of Xray radiation, for which the medium of the wind remains mostly transparent. The typical column densities of the high/soft and low/hard state are 10^{24} and 10^{25} particles per cm^{2}, as indicated in Fig. 6. For the solar abundances, 10 keV bremsstrahlung spectrum and known Xray cross sections of heavier elements, we can estimate the optical depth τ of O, C and N atoms to be in the range of 10^{1}−10^{2} and 10^{0}−10^{1} for the high/soft and the low/hard state, respectively. Inversely, following the same assumptions, O, C and N become optically thin for the Xray radiation at the energies of 3.5, 2, 8 and 2, 3 keV in the high/soft state, and 12, 8 and 6 keV in the low/hard state.
The assumption of a smooth mCAK wind is an idealization that is obviously violated in the turbulent flow between the component stars indicated by our present results as well as by other similar hydrodynamic simulations. Clumping, porosity, and vorosity arising due to instabilities even in the symmetric winds may influence the dynamics of the outflow (cf. e.g. Owocki 2014). It should be thus taken into account in a selfconsistent treatment of radiation hydrodynamics of the circumstellar matter in binaries. The radiative transfer of both the optical radiation of the donor star and the Xray radiation from the compact companion and the interplay of these two radiations through their interaction with the mass flow should be treated in a better approximation.
Acknowledgments
The authors are grateful to S. Owocki and R. Wünsch for useful comments and suggestions. This work was supported in part by the Czech Science Foundation (GAČR 1437086G) – Albert Einstein Center for Gravitation and Astrophysics, and grant SVV260089.
References
 Abbott, D. C. 1982, ApJ, 259, 282 [NASA ADS] [CrossRef] [Google Scholar]
 Blondin, J. M. 1994, ApJ, 435, 756 [NASA ADS] [CrossRef] [Google Scholar]
 Blondin, J. M., Kallman, T. R., Fryxell, B. A., & Taam, R. E. 1990, ApJ, 356, 591 [NASA ADS] [CrossRef] [Google Scholar]
 Bondi, H., & Hoyle, F. 1944, MNRAS, 104, 273 [NASA ADS] [CrossRef] [Google Scholar]
 CaballeroNieves, S. M., Gies, D. R., Bolton, C. T., et al. 2009, ApJ, 701, 1895 [NASA ADS] [CrossRef] [Google Scholar]
 Castor, J. L. 1974, MNRAS, 169, 279 [NASA ADS] [CrossRef] [Google Scholar]
 Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Cranmer, S. R., & Owocki, S. P. 1995, ApJ, 440, 308 [NASA ADS] [CrossRef] [Google Scholar]
 Fransson, C., & Fabian, A. C. 1980, A&A, 87, 102 [NASA ADS] [Google Scholar]
 Friend, D. B., & Abbott, D. C. 1986, ApJ, 311, 701 [NASA ADS] [CrossRef] [Google Scholar]
 Gayley, K. G. 1995, ApJ, 454, 410 [NASA ADS] [CrossRef] [Google Scholar]
 Hadrava, P., & Čechura, J. 2012, A&A, 542, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hadrava, P., & Čechura, J. 2013, in IAU Symp. 290, eds. C. M. Zhang, T. Belloni, M. Méndez, & S. N. Zhang, 219 [Google Scholar]
 Hatchett, S., & McCray, R. 1977, ApJ, 211, 552 [NASA ADS] [CrossRef] [Google Scholar]
 Hatchett, S., Buff, J., & McCray, R. 1976, ApJ, 206, 847 [NASA ADS] [CrossRef] [Google Scholar]
 Kallman, T. R., & McCray, R. 1982, ApJS, 50, 263 [NASA ADS] [CrossRef] [Google Scholar]
 Manousakis, A., Walter, R., & Blondin, J. M. 2012, A&A, 547, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Orosz, J. A., McClintock, J. E., Aufdenberg, J. P., et al. 2011, ApJ, 742, 84 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. 2014 [arXiv:1409.2084] [Google Scholar]
 Owocki, S. P., Castor, J. I., & Rybicki, G. B. 1988, ApJ, 335, 914 [NASA ADS] [CrossRef] [Google Scholar]
 Pauldrach, A., Puls, J., & Kudritzki, R. P. 1986, A&A, 164, 86 [NASA ADS] [Google Scholar]
 Pelupessy, I., Lamers, H. J. G. L. M., & Vink, J. S. 2000, A&A, 359, 695 [NASA ADS] [Google Scholar]
 Reid, M. J., McClintock, J. E., Narayan, R., et al. 2011, ApJ, 742, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Stevens, I. R., & Kallman, T. R. 1990, ApJ, 365, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Tarter, C. B., Tucker, W. H., & Salpeter, E. E. 1969, ApJ, 156, 943 [NASA ADS] [CrossRef] [Google Scholar]
 von Zeipel, H. 1924, MNRAS, 84, 665 [NASA ADS] [CrossRef] [Google Scholar]
 Ziółkowski, J. 2014, MNRAS, 440, L61 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Ratio of the tangential to radial component of the radiative force in the plane of stellar centres and rotational axis. Mass ratio q = 2 / 3, volume of the Roche equipotential (drawn as a light blue line) 0.95 of the critical Roche lobe. Limb and gravitydarkening (0.0, 0.0) and (0.6, 0.25) are plotted in the left and righthand panel. The black isocontours correspond to the values of − 3, − 2, and − 1. 

Open with DEXTER  
In the text 
Fig. 2 Series of three twodimensional simulations of the stellar wind in the equatorial plane of Cygnus X1 with an increasing α parameter. These simulations disregard the effects of the Xray ionization on the dynamics of the wind. Additionally, the k parameter is kept constant at 0.25. A clear trend of a shrinking gaseous tail and decreasing matter captured in the vicinity of the compact companion with increasing α is evident. The three upper panels shows the density distribution after the solution reached the quasistationary state, ~3P_{orb}. The lower panels represent corresponding velocity magnitude plots. (This figure is available in colour in the electronic version of the paper.) 

Open with DEXTER  
In the text 
Fig. 3 Integrated column density for the simulation shown in Fig. 2 as a function of orbital phase Φ timeaveraged over one orbital period. Each line represents the column density profile for different values of α – 0.5 (dashed line), 0.6 (dashdotted line), and 0.7 (solid line). 

Open with DEXTER  
In the text 
Fig. 4 Similarly as in Fig. 2, this series of six twodimensional simulations in the equatorial plane shows the dependency of the solution on the mass ratio M_{x}/M_{∗} and the mean radius of the primary. The masses of the primary and the compact companion for all three cases can be found in Table 2. The shape of the primary, represented by an equipotential surface, is a function of M_{∗} and M_{x}. The white index in the upperright corner of each panel represents the volume of the primary as a fraction of volume of the critical Roche equipotential. We disregard the effects of the Xray ionization in these simulations. α is set to 0.6 and the k parameter is kept constant at 0.25. Changes among the simulations become more pronounced as we take into account the scale differences. For a given value of P_{orb}, D grows with increasing overall mass of the system. The panels show the density distribution after the solution reached the quasistationary state, ~3P_{orb}. (This figure is available in colour in the electronic version of the paper.) 

Open with DEXTER  
In the text 
Fig. 5 Twodimensional simulations of the stellar wind in the equatorial plane of Cygnus X1 including the effects of the Xray ionization. Panels show the density distribution ρ, the velocity magnitude , and the ionization parameter ξ. α is set to 0.6, and k = k(ξ). The wind launching from the hemisphere facing the Xray source is noticeably slower than the wind coming from the opposite hemisphere, which lies in the Xray shadow (the dark conelike region in the right panels). A strong shock forms in the region where the fast and slow wind collide. The accretion of the slow wind is much more efficient, giving rise to an extensive disk that is enveloped by a strong shock front. In the region where log ξ ≥ 2, the linedriven acceleration mechanism of the material of the wind is switched off, which causes the wind coming from the hemisphere facing the Xray source to slow down. (This figure is available in colour in the electronic version of the paper.) 

Open with DEXTER  
In the text 
Fig. 6 Integrated column density for the simulation shown in Fig. 5 as a function of orbital phase Φ timeaveraged over one orbital period. The solid and dashed lines correspond to the high/soft and low/hard Xray state of Cygnus X1, respectively. 

Open with DEXTER  
In the text 
Fig. 7 Threedimensional simulation of the stellar wind of HD 226868 with an Xray luminosity equal to 1.9 × 10^{37} erg s^{1} corresponding to the low/hard Xray state of Cygnus X1. α is set to 0.6. η = η(ξ), and k = k(ξ) are functions of ξ. Columns show various cross sections throughout the computational volume. All mutually perpendicular planes are centred on the compact companion and are defined as x = D, y = 0, z = 0. The displayed quantities are the density distribution (the top panels), the velocity magnitude (the middle panels), and the ionization parameter ξ (bottom panels). (This figure is available in colour in the electronic version of the paper.) 

Open with DEXTER  
In the text 
Fig. 8 Threedimensional simulation of the stellar wind of HD 226868 with an Xray luminosity equal to 3.3 × 10^{37} erg s^{1} corresponding to the high/soft Xray state of Cygnus X1. α is set to 0.6. η = η(ξ), and k = k(ξ) are functions of ξ. Columns show various cross sections throughout the computational volume. All mutually perpendicular planes are centred on the compact companion and are defined as x = D, y = 0, z = 0. The displayed quantities are the density distribution (the top panels), the velocity magnitude (the middle panels), and the ionization parameter ξ (bottom panels). (This figure is available in colour in the electronic version of the paper.) 

Open with DEXTER  
In the text 
Fig. 9 Evolution of the density distribution in the orbital plane of Cygnus X1 showing gradual decay of the highdensity flow in the proximity of L_{1} point. The lefthand side panel corresponds to the quasistationary solution representing the low/hard Xray state of Cygnus X1. The righthand panel shows the outcome of the simulation when the new equilibrium is reached after adjusting the Xray luminosity of the compact object to the level corresponding to the high/soft state. The central panel represents a transient state roughly in the middle of the transition. The time index in the upper right corner of each panel is expressed in units of P_{orb}. (This figure is available in colour in the electronic version of the paper.) 

Open with DEXTER  
In the text 
Fig. 10 Interactive threedimensional isodensity surface for ρ = 2.5 × 10^{9} kg m^{3}. 

Open with DEXTER  
In the text 