X-RAY PHOTOIONIZED BUBBLE IN THE WIND OF VELA X-1 PULSAR SUPERGIANT COMPANION

For supergiant parameters given in Table 1 we first calculated the NLTE wind model with CMF line force neglecting the influence of the companion star. The emergent surface flux is taken from the H–He spherically symmetric NLTE model stellar atmospheres (Kubát 2003). The predicted wind mass-loss rate 1.5 × 10⁻⁶ M_☉ year⁻¹ agrees well with its estimate derived from the observed X-ray spectrum 1.5–2 × 10⁻⁶ M_☉ year⁻¹ (Watanabe et al. 2006). The predicted terminal velocity 750 km s⁻¹ is lower than the observed one 1100 km s⁻¹ (Prinja et al. 1990). This disagreement may stem either from the model simplifications, inaccurate stellar parameters (e.g., metallicity), or from the sensitivity of the wind terminal velocity to a detailed ionization balance in the outer wind regions (Puls et al. 2000), which probably manifests itself in a significant scatter of the ratio of the terminal to the escape velocity for stars with the same effective temperatures (e.g., Crowther et al. 2006).

The velocity structure of the wind model without X-ray irradiation can be approximated by

$$\begin{eqnarray} \tilde{v}(r)&=&\left[{v_1\left({1-\frac{R_*}{r}}\right) +v_3\left({1-\frac{R_*}{r}}\right)^3}\right]\nonumber\\ &&\times\left\lbrace {1-\exp \left[{\gamma \left({\frac{r}{R_*}-1}\right)^2}\right]}\right\rbrace, \end{eqnarray} \tag{ 1 }$$

where

$$\begin{equation} v_1=1042\ {\rm km}\ {\rm s}^{-1},\quad v_3=-297\ {\rm km}\ {\rm s}^{-1},\quad \gamma =-1220. \end{equation} \tag{ 2 }$$

Note that the representation of the velocity law by a polynomial expansion provides better approximation than the ordinary "β-velocity law," because these polynomials may form a functional basis (Krtička & Kubát 2011). The exponential term is included for a better fit of the velocity law close to the sonic point. The calculated X-ray opacity per unit mass averaged for radii 1.5 R_*–5 R_* may be approximated by

$$\begin{equation} \tilde{\kappa }^{\rm X}(\nu)= \left\lbrace \begin{array}{@{}c}a_1(\lambda -b_1)^2,\quad \lambda <\lambda _1,\\ a_2(\lambda +b_2)^3,\quad \lambda >\lambda _1,\\ \end{array}\right. \end{equation} \tag{ 3 }$$

where λ = 10⁸c/ν, a₁ = 0.704 g⁻¹ cm², b₁ = 1.056, a₂ = 4.06 × 10⁻³ g⁻¹ cm², b₂ = 11.41, and λ₁ = 20.18. We stress that λ enters as a non-dimensional parameter here, which has for convenience the same value as the wavelength in units of Å.

The full treatment of the problem would essentially require a complex solution of three-dimensional time-dependent hydrodynamic equations (Friend & Castor 1982; Blondin et al. 1990; Feldmeier et al. 1996). However, because the typical wind flow time $R_*/v_\infty \approx 0.3\,{\rm day{\rm s}}$ is roughly a factor of 30 shorter than the orbital period, we neglect the influence of the orbital motion on the wind structure. Consequently, we assume that the stellar wind of the supergiant component is axisymmetric with respect to the binary axis connecting the centers of components. Moreover, we assume that the stellar wind in the direction given by the inclination ϕ from the binary axis can be modeled by a spherically symmetric wind (see Figure 1). Consequently, in the orbital plane of the binary our two-dimensional wind model consists of sectors of a circle.

Zoom In Zoom Out Reset image size

The wind model in each sector is calculated using our NLTE wind code. Because the radial wind velocity may be non-monotonic in some sectors, we do not use the CMF line force for the calculation of these models directly, however the line radiative force is given by the force calculated in the Sobolev approximation (e.g., Castor 1974) multiplied by the ratio of the CMF to the Sobolev line force derived from the wind model that neglects the radiation from the companion star (Section 2.1).

The influence of the neutron star is taken into account only by inclusion of its X-ray radiation due to the wind accretion on the neutron star. This radiation irradiates the supergiant and interacts with its wind. To describe this effect, we add an additional term J^X(ν, r, d) to the specific intensity J(ν, r, d) in the form

$$\begin{equation} J^{\rm X}(\nu,r,d)=\frac{L^{\rm X}(\nu)}{16\pi ^2d^2}{\it e}^{-\tau (\nu,r,d)}, \end{equation} \tag{ 4 }$$

where the optical depth along the given ray is (z measures the distance along this ray)

$$\begin{equation} \tau (\nu,r,d)=\int _0^d\kappa (\nu,z)\rho (z){\it d} z, \end{equation} \tag{ 5 }$$

L^X(ν) is the X-ray luminosity per unit of frequency, d is the distance from the given point in the supergiant wind region to the surface of the neutron star, κ(ν, z) is the absorption coefficient per unit of mass, and ρ(z) is the wind density. The distribution of emergent X-rays L^X(ν) is approximated by the power law (Watanabe et al. 2006). Energies higher than 3 keV, which are well above the ionization energies of all included ions, were not considered in the model. The absorption coefficient and the density in Equation (5) can be derived from models for individual sectors. However, to simplify our approach, for the calculation of J^X(ν) we use fits following from Equations (1) and (3):

$$\begin{eqnarray} \hspace{-2.2pc}\rho (z)=&\frac{\dot{M}}{4\pi r^2 v(z)}, \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} v(z)=&\min (\tilde{v}(r),\hat{v}(\phi)), \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} \kappa (\nu,z)=&\tilde{\kappa }^{\rm X}(\nu), \end{eqnarray} \tag{ 8 }$$

where the relation between the distance along the ray z and the radius r is derived from the geometry of the problem, and $\hat{v}(\phi)$ is an average velocity in the velocity plateau which occurs due to X-rays (see Figure 3).

The models describing the wind in different sectors with inclination ϕ with respect to the binary axis are calculated for a sequence in ϕ with a step of 10°.

Part of the mass which leaves the primary star via its wind can be accreted onto the compact component. Its potential energy is then released in the form of X-ray radiation. Let us adopt an approximate Hoyle–Lyttleton treatment to estimate the X-ray luminosity. For a moving medium, the radius from which the matter can be accreted onto a point mass is given by the relation r_HL = 2GM_X/v², where M_X is the mass of the compact component (neutron star), v² = v²_wind + v_orb², v_orb is the orbital velocity of the compact component, and v_wind is the wind velocity at the distance of the compact component (denoted by D). Given v_orb = 300 km s⁻¹ and our simulated wind structure, we find that the neutron star is able to collect matter only from a narrow cone defined by the value of ϕ < 15 ° (for ϕ = 15 ° and the distance D we have v_wind = 32.0 km s⁻¹ and r_HL = 0.17D). We can express the relation between the accretion rate $\dot{M}_{\rm acc}$ and the mass-loss rate $\dot{M}$ from the supergiant as (Watanabe et al. 2006) $\dot{M}_{\rm acc}=\dot{M} {r_{\rm HL}^2}/(4D^2)$. Then the X-ray luminosity $L_{\rm X}={GM_{\rm X}\dot{M}_{\rm acc}}/{R_{\rm X}}$ is L_x = 8.8 × 10³⁷ erg s⁻¹. This value is approximately one order of magnitude higher than the observed one (Watanabe et al. 2006). The reason for this difference is the fact that only a small fraction of matter in the accretion cone gets finally accreted and the rest falls back to the surface of the primary as a consequence of the flow stagnation between the supergiant and the neutron star.

The wind equations allow the existence of two types of solutions giving different X-ray luminosities and wind velocities. The solution presented here appears in the case of a strong X-ray source, which significantly affects the wind ionization state and consequently also the radiative force. This results in a slow wind that can be accreted by the neutron star in large amounts, producing a strong X-ray source.

Another type of solution may occur in the case of a weak X-ray source that does not significantly influence the wind ionization state. The radiation force becomes higher, similar to that without any X-ray irradiation. This results in faster outflow roughly corresponding to the "no X-rays" case in Figure 3. Due to the dependence of the accretion rate on the velocity via $\dot{M}_{\rm acc}\sim v^{-4}$ and having roughly two times higher v, the X-ray luminosity is an order of magnitude lower in this case (since $L_{\rm X}\sim \dot{M}_{\rm acc}$). This agrees with the adopted assumption of a weak X-ray source. Consequently, this is also a possible solution of the wind equations.

Different types of solutions may appear in different systems, or even an external perturbation may cause switching between these two wind solutions in a particular binary system, possibly contributing to the variability of X-ray luminosity. This variability might be accompanied by the variability of the distribution of emitted X-rays, if the accretion regime changes with accretion rate (Shakura et al. 2012).

For a slightly lower mass-loss rate or for a slightly higher X-ray luminosity than assumed here the X-rays could penetrate deeply into the stellar wind and significantly influence the ionization state at the wind base, resulting in the decrease of the radiative force. This could lead to the disruption of the stellar wind and a significant decrease of X-ray luminosity, possibly providing another contribution to the overall X-ray variability observed in Vela X-1 (Kreykenbohm et al. 2008).

This means that there exists a maximum X-ray luminosity the compact star can have for a given geometry of the system and wind mass-loss rate. Assuming that the Vela X-1 luminosity is the maximum one, the maximum allowed X-ray luminosities for other HMXBs L^max_X can be derived using the optical depth (compare with the ionization parameter ξ; Tarter et al. 1969) $\tau _{\rm X}=\int _{R_*}^D\kappa \rho {\it d}r \sim \dot{M}(D-R_*)/(v_\infty DR_*)$ and corresponding parameters of Vela X-1 as

$$\begin{equation} L_{\rm X}^{\rm max}e^{-\tau _{\rm X}}= L_{\rm X}(\mbox{Vela X-1})e^{-\tau _{\rm X}(\mbox{\scriptsize Vela X-1})}, \end{equation} \tag{ 9 }$$

or, in scaled quantities

$$\begin{eqnarray} &&\log (L_{\rm X}^{\rm max}/\mbox{1 erg})=32.6+3.9 {\left({\frac{\dot{M}}{1.5\times 10^{-6}\,{M}_{\odot }\,{\rm year}^{-1}}}\right)}\nonumber\\ &&\quad\times {\left({\frac{v_\infty }{750\,{\rm km}\,{\rm s}^{-1}}}\right)^{-1} \left({\frac{(DR_*)/(D-R_*)}{68.5\,{R}_\odot }}\right)}^{-1}. \end{eqnarray} \tag{ 10 }$$

Here $\dot{M}$, v_∞, R_*, and D are the wind mass-loss rate, the terminal velocity, the radius of the luminous component, and the binary separation for individual HMXBs.

The results of observations for other HMXBs collected in Figure 4 clearly support the picture that there exists a maximum allowed X-ray luminosity, which depends on the wind and geometry parameters. Some systems lie close to the boundary of the forbidden area, whereas others have typically higher wind mass-loss rate $\dot{M}\gtrsim 5\times 10^{-6} \,{M}_\odot \,{\rm year}^{-1}$ that keeps them away from the boundary. The position of individual stars in this diagram may vary with time due to the existence of two possible solutions of the wind equations (as discussed in Section 4.2).

Zoom In Zoom Out Reset image size

Our results support the mass-loss rate predictions based on up-to-date wind models in two ways. First, our mass-loss rate prediction of HD 77581, 1.5 × 10⁻⁶ M_☉ year⁻¹, agrees with the value estimated from X-ray spectroscopy, 1.5–2 × 10⁻⁶ M_☉ year⁻¹ (Watanabe et al. 2006). Moreover the wind mass-loss rate cannot be lower than this value, because a decrease of the wind mass-loss rate would lead to lower X-ray opacity, stronger wind X-ray photoionization close to the star, and even more significant reduction of the radiative force. This would finally cause a disruption of the wind and a disappearance of X-ray emission. This imposes a strong lower limit on the observational wind mass-loss rate estimates, which is in agreement with current mass-loss rate predictions (Vink et al. 2001; Krtička & Kubát 2010).

The consistent inclusion of the X-ray irradiation, which is an advantage of our models, also determines their shortcomings. The coupled solution of NLTE and radiative transfer equations is significantly time-consuming even in one dimension. Its inclusion in multidimensional time-dependent simulations is likely beyond the possibilities of current computers. Consequently, while one part of the problem solution is treated in detail, the second one is simplified. The correct picture of the flow in the HMXBs should take into account the results of both approaches: the stagnation of the flow in the direction toward the neutron star, which is studied in this paper, and a complex two-dimensional picture of the wind accretion on the neutron star (Blondin et al. 1990; Blondin & Woo 1995; Feldmeier et al. 1996).

On the other hand, there are effects that are not described by any of the available models. It is well established that the hot star wind is inhomogeneous on small scales (clumped; see Hamann et al. 2008). In the case of HMXBs, clumping (which favors recombination) may affect the region in which the photoionized bubble is formed (Oskinova et al. 2012). However, the hydrodynamical simulations (Feldmeier et al. 1997; Runacres & Owocki 2002) predict that clumping starts above the critical point of the wind solution; consequently, it does not affect wind mass-loss rates and terminal velocities. Thus, we expect that clumping does not significantly influence the results of our models.

The wind inhomogeneities likely cause the high variability of the X-ray source (Fürst et al. 2010; Oskinova et al. 2012). The calculation of the wind ionization should in fact account for such time-dependent X-ray photoionization. However, because the typical flow time of the wind is longer than the typical timescale of X-ray variability, we expect that our models are able to reproduce the mean effect of the X-ray photoionization. On the other hand, the variable X-ray ionization source causes temporal changes of the radiative acceleration, providing external perturbation which may contribute to the natural wind clumping. As a result, the primary wind may be more clumpy than the wind of a similar single supergiant.