© ESO, 2017

1. Introduction

The high-mass X-ray binary (HMXB) Vela X-1 consists of a ~283 s period pulsar (McClintock et al. 1976) in an eclipsing ~9 d orbit with the B0.5 Ib supergiant HD 77581 (Hiltner et al. 1972; Forman et al. 1973). The radius of HD 77581 is 30 R_⊙ and the binary separation is 53.4 R_⊙ or ~1.7 R_{HD 77581} (van Kerkwijk et al. 1995; Quaintrell et al. 2003) meaning that the neutron star is deeply embedded in the strong (mass-loss rate ~ 2 × 10^-6M_⊙ yr^-1, Watanabe et al. 2006) wind of the companion that it accretes from. Newest estimates consistently point towards a low terminal velocity, v_∞, of 600–750 km s^-1 (Krtička et al. 2012; Giménez-García et al. 2016; Sander et al. 2017b). Despite having a modest average X-ray luminosity of ~ 4 × 10³⁶ erg s^-1, with a distance of ~2 kpc (Nagase et al. 1986) Vela X-1 is close enough that it is one of the brightest HMXBs in the sky.

The high inclination of the system (>73°,Joss & Rappaport 1984) naturally allows probing of the accretion and wind geometry through observations at different orbital phases (e.g., Haberl & White 1990; Goldstein et al. 2004; Watanabe et al. 2006; Fürst et al. 2010). Outside of the eclipse, the base level of absorption in the system shows a systematic evolution along the binary orbit, with average column densities ranging from a few times 10²² cm^-2 at φ_orb ≈ 0.25 to 15 × 10²² cm^-2 at φ_orb ≈ 0.75 (Doroshenko et al. 2013). The evolution does not agree with expectations for an unperturbed wind and is attributed to the changing sightline through the complex accretion geometry of the system that includes an accretion wake, a photoionization wake, and possibly a tidal stream (e.g., Eadie et al. 1975; Nagase et al. 1983; Sato et al. 1986; Blondin et al. 1990; Kaper et al. 1994; van Loon et al. 2001; Malacaria et al. 2016).

Superimposed on this change are irregular individual absorption events (Haberl & White 1990; Odaka et al. 2013) and bright flares (Martínez-Núñez et al. 2014) of varying length as well as so-called off-states (Kreykenbohm et al. 2008; Sidoli et al. 2015). Clumpiness of the companion wind has been invoked by several authors to explain the short-term variability through absorption in or accretion of clumps (Kreykenbohm et al. 2008; Fürst et al. 2010; Martínez-Núñez et al. 2014), while Manousakis & Walter (2015) invoke unstable hydrodynamic flows. At least the off-states can also be explained by a quasi-spherical subsonic accretion model (Shakura et al. 2013). The seminal simulations of stellar-wind disruption by a compact source of Blondin et al. (1990) show an oscillating accretion wake structure that also leads to variability in the accretion rate.

High-spectral-resolution observations allow us to address the structure and distribution of material in the system: as plasmas of different densities, temperatures, and velocities reprocess the radiation from the neutron star, they imprint characteristic features onto the X-ray continuum. Previous analyses of the existing Chandra High Energy Transmission Grating Spectrometer (Chandra-HETGS; Canizares et al. 2005) observations concentrated on the global properties of the plasma: emission lines from a variety of ionization stages of Ne, Mg, Si, and Fe, including Kα emission from L-shell ions in some cases, have been detected in eclipse (φ_orb ≈ 0, Schulz et al. 2002) and when looking through the accretion wake (φ_orb ≈ 0.5, Goldstein et al. 2004; Watanabe et al. 2006). These observations allow us to determine the physical conditions in this system, including the properties of the wind. Schulz et al. (2002) showed that a hot optically thin photoionized plasma and a colder plasma that gives rise to the fluorescent Kα lines from L-shell ions, co-exist in the system. Comparing fluorescent emission lines at eclipse and at φ_orb ≈ 0.5, Watanabe et al. (2006) find them an order of magnitude brighter at φ_orb ≈ 0.5. They thus conclude that the lines have to be produced in a region that is occulted during the eclipse and thus lies between neutron star and companion.

Surprisingly, data obtained at φ_orb ≈ 0.25, when our line of sight towards the approaching neutron star is relatively unobscured (Fig. 1), proved more difficult to analyze and interpret. Goldstein et al. (2004) describe emission and absorption lines with a strong continuum contribution. They conclude that these lines must come from different regions, but do not put constraints on the location of the material. Watanabe et al. (2006) do not discuss the spectrum at this orbital phase in detail as they conclude that it is too continuum-dominated for their purposes, but they mention possible evidence for weak absorption lines. In their 3D Monte Carlo simulations of X-ray photons propagating through the wind, the model shows P Cygni profiles at φ_orb ≈ 0.25 that they conclude to not be observable with Chandra given the spectral resolution and sensitivity.

At the same time, recent results with Suzaku (Odaka et al. 2013) and XMM-Newton (Martínez-Núñez et al. 2014) at orbital phases of ~0.25 show strong changes in absorption on timescales as short as ~ks. Similar behavior has also been indicated in historical observations with Tenma and EXOSAT (Nagase et al. 1986; Haberl & White 1990). However, the Chandra observations of Vela X-1 have only previously been analyzed in a time-averaged way, disregarding possible spectral changes within a given observation. This shortcoming prompted us to re-visit an early ~30 ks-long Chandra-HETGS observation of Vela X-1 at φ_orb ≈ 0.25 with the aim to search for variability on shorter timescales using Chandra’s high-resolution capabilities.

We first introduce our data, discuss the overall variability during our observation and motivate the division of the total exposure into periods of low and high hardness, which correspond to low and high absorption and are to be analyzed separately, in Sect. 2. In Sect. 3, we analyze the line content of the resulting spectra. In Sect. 4, we discuss the observations with a focus on the wind structure and the interactions of the compact object with the wind. We close with a summary and an outlook towards future observations and laboratory work necessary to further our understanding of the structured wind and clumpy absorber in Vela X-1 in Sect. 5.

2. Data and variability behavior

We use the ObsID 1928 taken with Chandra-HETGS on 2001-Feb.-05 for ~30 ks at orbital phases φ_orb = 0.21–0.25 (according to the ephemeris of Kreykenbohm et al. 2008, where φ_orb = 0 is defined as mid-eclipse, Fig. 1). This ObsID has previously been included in the analyses of Goldstein et al. (2004), who showed that the spectrum is not affected by pile-up, and Watanabe et al. (2006); both extract one spectrum for the full ObsID only and Watanabe et al. (2006) do not analyze the line content of this observation.

We follow the standard Chandra threads and use CIAO 4.7 using CalDB 4.6.7, but employ a narrower sky mask to improve the flux at very short wavelengths by avoiding some of the overlap between extraction regions. For bright sources, such as X-ray binaries in general and Vela X-1 in particular, Chandra-HETG observations can be considered as essentially background-free. All analyses in this paper were performed with the Interactive Spectral Interpretation System (ISIS) 1.6.2 (Houck & Denicola 2000; Houck 2002; Noble & Nowak 2008). Uncertainties are given at the 90% confidence level for one parameter of interest.

2.1. Light curves and hardness

Figure 2 shows light curves extracted with the approximate time resolution of the neutron star spin period (283 s) in the 0.5–3 keV and 3–10 keV energy bands and the corresponding hardness, defined as the ratio between the harder and the softer band. We can distinguish two distinct behaviors: periods of highly variable high hardness and of stable low hardness. Interestingly, the total count rate changes by a factor of a few during the period of stable low hardness. Such increases in flux can be attributed to increases in the intrinsic X-ray brightness of the neutron star (e.g., Martínez-Núñez et al. 2014). To investigate the origin of the changing hardness, we time-filtered our observation and extracted one spectrum each for the low- and high-hardness periods using user-defined good time intervals corresponding to the two periods highlighted in Fig. 2. Exposures and total counts are listed in Table 1.

2.2. Continuum variability

We first address changes in the overall continuum shape between low- and high-hardness periods. A physical model for X-ray continuum emission in accretion-powered X-ray pulsars is an open question (but see Becker & Wolff 2007; Farinelli et al. 2016; Wolff et al. 2016, for advanced approaches). Thus a multitude of empirical models, usually some combination of absorbed power laws, has been successfully used to describe spectra of accreting pulsars (e.g., Müller et al. 2013).

In Vela X-1, the presence of additional possibly X-ray emitting and reprocessing components such as the stellar wind, the denser material in the accretion, and ionization wakes, as well as the surface of HD 77851, make modeling the soft X-ray emission especially challenging (see Fig. 3 for an overview of the 0.8–9.5 keV spectrum). Below ~2.5 keV, the spectrum is complex and shows strong line features and peculiar continuum components such as a soft excess (Nagase et al. 1986; Hickox et al. 2004; Martínez-Núñez et al. 2014) that depend on orbital phase and absorption (Haberl & White 1990). We therefore define the continuum based on the 2.5–10 keV band, which does not include any strong lines with the exception of the Fe Kα complex (Torrejón et al. 2010; Giménez-García et al. 2015), which we exclude by ignoring the 6.2–6.5 keV band. To improve the signal-to-noise ratio (S/N), we combine the positive and negative first-order High Energy Grating (HEG) and Medium Energy Grating (MEG) data after rebinning to the MEG grid. We then rebin the combined spectrum to a S/N of 5 and a minimum of four MEG channels per bin.

We describe the high- and low-hardness spectra each with a power law, defined by a power law index Γ and a normalization that is modified by absorption modeled with tbnew, the new version of tbabs (Wilms et al. 2000), using wilm abundances (Wilms et al. 2000) and vern cross-sections (Verner et al. 1996). The assumed model is simple and will not result in a perfect description of the continuum, but given the quality of the data it is adequate for our aim to assess continuum-shape changes.

To assess whether changes in absorption, expressed through the equivalent hydrogen column density, N_H, or in power-law slope, Γ, lead to the hardness changes, we perform two fits, with results listed in Table 2. In a first fit, we assume that the absorption column, described by the N_H parameter, is the same for both low and high hardness and tie this parameter between the two hardness modes. No acceptable fit can be obtained in this case. In a second approach, we assume the same power-law shape, that is, we tie Γ, but let the absorption be different. This model results in a better description of the data (Fig. 4).

Since there is still some systematic deviation between data and model above ~7 keV for the high-hardness spectrum, we also fit the two spectra with all model parameters independent of each other (Table 2). In Fig. 5, we show the confidence intervals for Γ and N_H for this case as well as the 90% allowed ranges for Γ and N_H for the case when the power law slope is tied between the two spectra. The confidence contours show the well known degeneracy between the derived values of Γ and N_H (e.g., Suchy et al. 2008). We also see that letting Γ free leads to an even larger difference in hydrogen column density between the two spectra than when Γ is tied. The bulk of the changes in overall spectral shape are thus due to absorption and not due to changes in the shape of the underlying power law continuum, although some changes in spectral shape may also be present.

3. High-resolution spectroscopy

Having established that the changes in overall continuum shape between the low- and high-hardness periods are driven by changes in absorption, we now turn to high-resolution spectroscopy where we focus on regions below 2.5 keV (above 5 Å) as well as on the Fe Kα complex, which were not included in the broadband spectral modeling due to their complexity. We first visually inspect the differences between the low- and high-hardness spectra as presented in Fig. 3 and compare to the results of Goldstein et al. (2004) who analyzed the spectrum of ObsID 1928 without subdividing it into hardness-dependent periods. Because of the larger exposure and flux of the low-hardness data (Table 1), we expect that this period has dominated their analysis. We note a difference in line contributions between the spectra at different hardnesses: features change their relative strengths (e.g., lines of the He-like Si xiii triplet at ~6.7 Å, see Sect. 3.1), change from emission to absorption (e.g, Ne x Ly β at ~10.25 Å, see Sect. 3.3), and new features appear that were not found in the earlier analysis, such as lines from low-ionization stages of Si around ~7 Å (see Sect. 3.1).

In the remainder of this section, we model the low- and high-hardness spectra independently. As in Sect. 2.2, we combine the positive and negative first-order HEG and MEG data after rebinning to the MEG grid, but do not rebin the spectra further. Except where explicitly stated otherwise, we do not assume different line features to originate at the same location or the same medium, especially since multiple gas phases with different temperatures have been suggested as explanations for observations at higher orbital phases (Goldstein et al. 2004; Watanabe et al. 2006). Our analysis of the lines is based on local power-law continuum fits, that is, power-law continua describing a narrow interval around the respective features of interest. Although the slope, Γ, and normalization of such power laws do not have a direct physical meaning, we list the best fit results obtained for all regions discussed in the remainder of this section in order to facilitate reproducibility of results in Table 3. Our choice of wavelength intervals used for such local fits is informed by the previous knowledge of lines observed in Vela X-1 at other orbital phases (Schulz et al. 2002; Goldstein et al. 2004; Watanabe et al. 2006). We consistently use Cash statistics (Cash 1979) in this section as several of the considered intervals require it. We always first fit for the width of any component. However, since none of the detected Gaussian emission components are resolved at MEG resolution (0.023 Å FWHM), we fix their widths, σ, to σ = 0.003 Å, that is, about one third of the MEG resolution.

We use two approaches to search for lines: a blind search based on the Bayesian Blocks algorithm¹ (Scargle et al. 2013) and a conventional, manual approach, where we utilize our knowledge of expected line energies of a given element. The Bayesian Blocks approach for spectral data is introduced, discussed, and benchmarked against other methods in Young et al. (2007). With this latter approach, we test for the existence of lines against a continuum model, starting with a power law that we iteratively fit to the data. At each iteration we extend the model by adding any lines found in the previous step.

The algorithm treats the data in an unbinned manner, but compares the changes in likelihood compared to the applied model² for grouping the data into different numbers of “blocks”, with a single block corresponding to the absence of lines. For a given penalty factor, the algorithm exactly solves the 2ⁿ⁻¹ problem, where n is the number of data bins, of how many blocks to divide the data into. A new block change point can start (or not) between any two bins, hence the 2ⁿ⁻¹ possible divisions of the spectrum. A line is considered to exist between two such change points. Implementations of the Bayesian Blocks algorithm introduce a penalty factor that imposes a requirement for a minimum increase of the log likelihood to consider significance of the addition of a change point to the description of the function (e.g., spectrum or light curve). This penalty factor is controlled by an input parameter denoted by α in the work of Young et al. (2007) and by γ in notation of Scargle et al. (2013). The increase of the log likelihood is approximately equivalent to requiring that the addition of a change point exceeds a significance of p ~ exp(−α). Since a line requires two such change points, choosing an input of α = 1.5 is equivalent to requiring that lines have p< exp(−2α) ≈ 0.05, that is, they are at least 95% significant, or, P> 1−exp(−2α) ≈ 0.95.

Our analysis closely follows that of Young et al. (2007) where we start with a large value of α, apply the algorithm, and decrease α in steps of 0.1 until a line is found. We stop at a lower limit of α = 1.5. We list in Tables 4–8 the value of α at which (and below) the indicated line is found and call it α_sig. The line significance is roughly interpretable as P ≈ 1−exp(−2α_sig), but we note that, for example, He-like triplets are often detected as one excess encompassing all three lines.

3.1. Silicon region

We first address the silicon line region (6.0–7.4 Å). We point out that new laboratory reference data are available for the K-shell transitions of L-shell ions of Si, enabling robust velocity measurements (Hell et al. 2013, 2016).

3.1.1. High-hardness spectrum of the silicon region

The Bayesian Blocks approach finds three emission features, all with α_sig> 3.3. We identify these lines with H-like Si xiv Lyα³, He-like Si xiii f forbidden line $\mathrm{(}\mathrm{1}\mathrm{s}\mathrm{2}\mathrm{s}\hspace{0.17em}{\mathrm{3}}^{}{\mathrm{S}}_{\mathrm{1}}$–$\mathrm{1}{\mathrm{s}}^{\mathrm{2}}\hspace{0.17em}{\mathrm{1}}^{}{\mathrm{S}}_{\mathrm{0}}$, also known as line z in the notation of Gabriel 1972) and a blend of the Kα lines from Si ii–Si iv. We model these features with Gaussians.

A visual inspection of the residuals, however, reveals further wave-shaped structures blueward of both Si xiii f and the Si ii–Si iv blend. We model these structures with one and two Gaussian lines, respectively, and identify them with the He-like Si xiii r resonance line ($\mathrm{1}\mathrm{s}\mathrm{2}\mathrm{p}\hspace{0.17em}{\mathrm{1}}^{}{\mathrm{P}}_{\mathrm{1}}$–$\mathrm{1}{\mathrm{s}}^{\mathrm{2}}\hspace{0.17em}{\mathrm{1}}^{}{\mathrm{S}}_{\mathrm{0}}$, also known as line w in the notation of Gabriel 1972) and blends of Kα transitions in N-like Si viii and O-like Si vii (Fig. 6).

In our best fit model that uses six Gaussian emission lines, large-scale residuals that encompass the area around ~7 Å are still present. We model them with a broad Gaussian (σ ≈ 0.05 Å) in absorption that we explicitly do not interpret as a true line component but as complex continuum shape that a single power law fails to describe. The introduction of the broad Gaussian does not significantly change the obtained line positions. We note especially that the Bayesian Blocks approach relies on a good description of the continuum; the complex continuum shape therefore contributed to the failure of the automated approach to find the Si xiii r line and the Kα lines of Si vii and Si viii. We do not re-run the Bayesian Blocks algorithm on a continuum including the additional curvature that we modeled, as this would be inconsistent with a blind search.

Fitting just one Doppler shift for all detected Si lines does not result in a satisfactory model for the narrow features. In the interest of reducing free parameters, it is reasonable to assume that the lines with a similar ionization stage originate in the same medium. The lack of intermediate-ionization lines naturally leads to two line groups: (1) Si xiv Lyα, Si xiii f, and Si xiii r; and Si viii and Si vii Kα. We therefore fit one line shift to each of these groups (see Table 4 for reference values). Using this approach, we can obtain a good description of the data.

Although we formally fit a Doppler shift for the Si ii–vi Kα line blend, we do not take it into consideration when tying lines or discussing the overall behavior: the contribution of lines to a blend (and thus the measured center of the blend) depends on the properties of the plasma that are certainly different in laboratory and astrophysical environments (Hell et al. 2016). In fact, Hell et al. in their NASA/GSFC electron beam ion trap (EBIT) Calorimeter Spectrometer (ECS) data can distinguish two components in the Si ii–iv Kα lines, whose center they see at 1742.03 ± 0.06 eV: a blend of F-like Si vi Kα and Ne-like Si v Kα at $\mathrm{1742.8}{\mathrm{8}}_{-0.17}^{\mathrm{+}\mathrm{0.15}}$ eV and a blend of Al- to Na-like Si ii–iv Kα at $\mathrm{1740.0}{\mathrm{4}}_{-0.36}^{\mathrm{+}\mathrm{0.27}}$ eV. Shifting the center of the full blend to 1740.04 eV would result in a change in measured velocity of >250 km s^-1. A further complication is the possible presence of the Mg xii Lyβ at 7.106155 Å or 1744.74 eV (Erickson 1977).

The resulting properties of the Si lines are summarized in Table 4. The fits are shown in Fig. 6, where we show them in comparison with low-hardness data and lab measurements, and in Fig. 7, where we also show fit residuals.

3.1.2. Low-hardness spectrum of the silicon region

In the low-hardness spectrum, we find one line with high significance (α_sig> 8.5) at 6.73–6.75 Å and a complex w-shaped structure extending 6.98–7.05 Å at α_sig = 1.5, that is, our lowest acceptable value of α_sig, with the Bayesian Blocks approach.

We identify the Gaussian emission feature as Si xiii f line (Table 4). The w-structure can be empirically described as a combination of emission and absorption lines at the same energy. This results in lines with widths σ_em ≈ σ_abs ≈ 0.013 ± 0.002 Å at a wavelength of 7.022 ± 0.003 Å and precludes a line identification (Figs. 6 and 7): if this were the Si vii Kα line blend, it would have a Doppler velocity of ~−1500 km s^-1, comparable to the one measured in the high-hardness spectrum. However, no other silicon lines except Si xiii f are present in the low-hardness spectrum and it would require extreme fine-tuning to produce the Si vii Kα line alone. The peculiar shape of the feature also remains a problem. Regarding other elements, while the Rydberg series of He-like Mg xi, H-like Mg xii and H-like Al xiii have some Kα transitions surrounding this region, none of them is sufficiently close or strong enough to explain this feature. Both H-like Mg xii Lyβ at 7.10616 Å and H-like Al xiii Lyα at 7.17272 Å (Erickson 1977) are too high in wavelength. While the He-like Mg xi series has several n ≥ 7 → 1 transitions between 7.174 Å (n = 7; Verner et al. 1996) and 7.037 Å (series limit or ionization potential; Drake 1988) there is no indication from lines earlier in the series that the lines should be strong enough to explain the observed feature. Still, a radiative recombination continuum (RRC) feature of Mg xi at 7.037 Å (Drake 1988) could contribute to the observed feature.

3.2. Magnesium region

We next address the 8.0–10 Å region, that is, from longwards in wavelength of Mg xii Lyα to shortwards of the Mg iv Kα blend. Several 4 → 2 transitions from Fe xxi to Fe xxiv fall into this region, but their location noes not agree with any of the features detected in the following, and we do not detect any of the corresponding stronger 3 → 2 transitions from Fe xxi to Fe xxiv (Brown et al. 2002) in Sect. 3.3.

3.2.1. High-hardness spectrum of the magnesium region

We identify the first line found (α_sig = 5.4) by the Bayesian Blocks approach with Mg xii Lyα (8.42101 ± 4.0 × 10^-5 Å, Erickson 1977). The algorithm further finds a very shallow step in the residuals at ~9.5 Å, consistent with the location of the magnesium K edge. However, the step implies an excess where, from an edge, we would expect additional absorption and, consistently, the inclusion of an edge does not improve our model. Indeed, the best-fit optical depth of the edge component is zero, that is, no edge is present, therefore we do not use this component further.

The next detected component consists of a broad excess spanning the whole He-like Mg xi triplet (r: $\mathrm{1}\mathrm{s}\mathrm{2}\mathrm{p}\hspace{0.17em}{\mathrm{1}}^{}{\mathrm{P}}_{\mathrm{1}}$–$\mathrm{1}{\mathrm{s}}^{\mathrm{2}}\hspace{0.17em}{\mathrm{1}}^{}{\mathrm{S}}_{\mathrm{0}}$; i: $\mathrm{1}\mathrm{s}\mathrm{2}\mathrm{p}\hspace{0.17em}{\mathrm{3}}^{}{\mathrm{P}}_{\mathrm{1}}$–$\mathrm{1}{\mathrm{s}}^{\mathrm{2}}\hspace{0.17em}{\mathrm{1}}^{}{\mathrm{S}}_{\mathrm{0}}$ and $\mathrm{1}\mathrm{s}\mathrm{2}\mathrm{p}\hspace{0.17em}{\mathrm{3}}^{}{\mathrm{P}}_{\mathrm{2}}$–$\mathrm{1}{\mathrm{s}}^{\mathrm{2}}\hspace{0.17em}{\mathrm{1}}^{}{\mathrm{S}}_{\mathrm{0}}$, also known as lines x and y in the notation of Gabriel 1972; f: $\mathrm{1}\mathrm{s}\mathrm{2}\mathrm{s}\hspace{0.17em}{\mathrm{3}}^{}{\mathrm{S}}_{\mathrm{1}}$–$\mathrm{1}{\mathrm{s}}^{\mathrm{2}}\hspace{0.17em}{\mathrm{1}}^{}{\mathrm{S}}_{\mathrm{0}}$) and with a stronger narrower excess at the Mg xi forbidden line. We include all three He-line triplet lines in our model, tying the distance between the triplet lines and the Lyα line to reference values from Drake (1988, for He-like) and Erickson (1977, for H-like). The last feature detected with the Bayesian Blocks approach is an emission feature at ~8.95 Å (α_sig = 2) that we cannot identify with an obvious line of magnesium or iron (Figs. 8, lower panel, and 9, left panel).

We parametrize the line strengths of the He-like triplet using the ratios R = f/i and G = (f + i) /r (Gabriel & Jordan 1969; Porquet & Dubau 2000) and model these density and temperature diagnostics directly. We obtain R = 1.5 ± 0.7 and $\mathit{G}\mathrm{=}\mathrm{2.}{\mathrm{0}}_{-0.7}^{\mathrm{+}\mathrm{3.5}}$ and a line shift of −140 ± 80 km s^-1 for the identified H- and He-like Mg lines (see Table 5 for detailed line parameters of the identified lines).

It has been recently shown by Mehdipour et al. (2015) that, under astrophysical conditions, He-like triplet lines of a given element can be affected by absorption from Li-like ions in the same medium. While Mehdipour et al. (2015) discuss this effect for the cases of iron and oxygen, it can easily be extended to other elements. Contamination of the He-like magnesium triplet in emission by absorption lines is clearly seen in the following only for the low-hardness spectrum (Sect. 3.2.2). While we do not detect any absorption lines in the high-hardness spectrum, we cannot exclude some contamination of the triplet.

3.2.2. Low-hardness spectrum of the magnesium region

The Bayesian Blocks approach finds discrete absorption features at ~9.15 Å (α_sig> 7), ~9.29 Å (α_sig> 4.6), ~9.7 Å and ~9.8 Å (α_sig = 2.0). It also finds a wide excess at ~9.1–9.3 Å, consistent with the location of the Mg xi triplet that becomes more pronounced as the narrow absorption features are iteratively taken into account by the model. We model the He-like triplet as in the case of the high absorption spectrum by fixing the relative location of the lines but leaving one line energy and the strengths of the individual lines free. The last component is a discrete absorption feature at ~9.77 Å (α_sig = 2.4).

The resulting model includes eight lines: the He-like triplet lines in emission and five absorption lines (Figs. 8, upper panels, and 9, right upper panel). In spite of suggestive residuals, no significant Mg xii Lyα is found and we leave it out in our best-fit model.

Since the absorption lines are weak and their widths unresolved, we are able to replace the narrow Gaussian absorption components used so far with a simple multiplicative box-function⁴ and thus fit their equivalent widths directly. For the emission lines of the He-like triplet we obtain velocities of −130 ± 60 km s^-1 (He-like triplet line parameters in Table 5).

The locations of the absorption lines are indicative of the different ionization stages of Mg. We use the theoretical reference values for the line energies of He-like Mg xi r to F-like Mg iv Kα lines from Behar & Netzer (2002) to identify the lines (Fig. 8) and calculate the line velocities (Table 6).

The Ne x Lyγ line at 9.70818 Å (Erickson 1977) is very close to N-like Mg vi Kα at 9.718 Å and indeed, Goldstein et al. (2004) have identified this absorption feature with Ne x Lyγ. While we detect Ne x Lyβ (Sect. 3.3.2), this line is blue-shifted by $\mathrm{-}\mathrm{38}{\mathrm{0}}_{-320}^{\mathrm{+}\mathrm{150}}$ km s^-1; a Ne x Lyγ line with a similar shift does not describe the observed absorption feature. We also do not detect the neighboring Ne x Lyδ at 9.48075 Å (Erickson 1977), although our upper limit on the equivalent width of this line is −4 × 10^-3 Å, that is, within the margins of our expectations of the strength of Ne x Lyδ given the equivalent width of Ne x Lyβ. On the other hand, assuming, as done here, that the two prominent absorption features at 9.7–9.8 Å are from Kα transitions in N-like Mg vi and O-like Mg v results in similar strengths for both lines, which both appear blue-shifted, consistent with the detected higher ionization Mg ions, confirming our line identification. We cannot, however, exclude a contribution of Ne x Lyγ to the observed feature we identify as Mg vi Kα in absorption. To test the amount of this contribution, we model the observed absorption feature with two components, one of them at the expected location of blue-shifted Ne x Lyγ. We obtain similar velocity and equivalent width for Mg vi Kα as without the neon contribution (v = −280 ± 90 km s^-1 and $\mathrm{-}\mathrm{5.}{\mathrm{1}}_{-1.6}^{\mathrm{+}\mathrm{1.8}}\mathrm{\times }{\mathrm{10}}^{-3}$ Å) and an equivalent width of $\mathrm{-}\mathrm{2.}{\mathrm{4}}_{-1.6}^{\mathrm{+}\mathrm{2.2}}\mathrm{\times }{\mathrm{10}}^{-3}$ Å for Ne x Lyγ.

The obtained velocity values for the absorption lines of Mg span a wide range of −(300–900) km s^-1 (Table 6) and are inconsistent with each other within their uncertainties. However, the complexity of the spectrum, in particular the overlap of emission and absorption features can lead to modeling problems that the formal errors do not account for. Additionally, inaccuracies in theoretical reference values may be present: while no EBIT laboratory data are yet available for magnesium, Hell et al. (2016) have shown that for silicon and sulfur, laboratory values and theoretical calculations of Behar & Netzer (2002) can disagree at up to 300 km s^-1, especially for lower ionization stages.

To better assess the contributions of the Mg absorption lines, we build a template (Mg xi to Mg iv Kα, Fig. 8 and Table 6) using the Behar & Netzer (2002) values: we assume the same line shift for all lines, but let the line strengths vary. For B-like Mg viii and C-like Mg vii, Behar & Netzer (2002) list two strong transitions. For each ion, both transitions are of comparable oscillator strength but are, respectively, 0.016 Å and 0.024 Å apart, which corresponds to a velocity difference of ~500 km s^-1 and ~750 km s^-1, respectively (cf. Fig. 8 where the line positions are marked). We incorporate these lines in our template using one line with the average line center as listed in Table 6, but note that this can lead to additional uncertainties in the determination of the line Doppler shifts. We also include the He-like Mg xi triplet in emission with parameters as described above and show the resulting best-fit model in Figs. 8, middle panel, and 9, lower-right panel. The line template parameters are listed in Table 6. The template line velocity of −360 ± 40 km s^-1 is in rough agreement with the average velocity of the five independent absorption lines. The equivalent widths of the template lines show that B-like Mg viii and F-like Mg iv Kα are consistent with 0, that is, the lines are not significantly detected, consistent with their non-detection in the Bayesian Blocks approach. The equivalent width of C-like Mg vii Kα may indicate a weak line. The large uncertainties on Li-like Mg x Kα, on the other hand, are due to the blending of the line with the intercombination and forbidden emission lines of the He-like Mg xi triplet.

3.3. Neon region

The region of interest encompasses the wavelengths 10–14 Å, that is, from H-like Ne x Lyβ to the He-like Ne ix triplet. At higher wavelengths, low source counts preclude an analysis and in particular the detection of possible Kα features of lower ionization stages than Ne viii. All identified neon lines in both low- and high-hardness spectra are listed in Table 7 together with their respective reference values.

3.3.1. High-hardness spectrum of the neon region

The Bayesian Blocks approach finds H-like Ne x Lyα with α_sig = 3.4, an excess spanning the location of the He-like Ne ix triplet with α_sig = 2.5, and finally the H-like Ne x Lyβ with α_sig = 2. No further structure is detected with this approach once these five lines are included.

We fit just one line shift for all five lines. Similarly to Sect. 3.2.1, we directly fit the He-like triplet line ratios, R and G (e.g., Porquet & Dubau 2000) and obtain $\mathit{R}\mathrm{=}\mathrm{0.}{\mathrm{8}}_{-0.6}^{\mathrm{+}\mathrm{1.4}}$ and $\mathit{G}\mathrm{=}\mathrm{4.}{\mathrm{3}}_{-2.9}^{\mathrm{+}\mathrm{5.7}}$. The resulting line velocity is −220 ± 70 km s^-1 (Fig. 10).

3.3.2. Low-hardness spectrum of the neon region

The Bayesian Blocks approach finds a complex feature that we identify with the He-like Ne ix triplet lines (α_sig = 5) and model with Gaussian components. The feature shows a dip between the Ne ix f and Ne ix i component and we model it with a narrowline component that we identify with the Li-like Ne viii Kα. The next strongest feature is H-like Ne x at α_sig = 4.6 that we fit with a Gaussian. We also note that the strong Fe xvii $\mathrm{1}{\mathrm{s}}^{\mathrm{2}}\hspace{0.17em}\mathrm{2}{\mathrm{s}}^{\mathrm{2}}\hspace{0.17em}\mathrm{2}{\mathrm{p}}_{\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}\mathrm{2}{\mathrm{p}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{4}}\hspace{0.17em}\mathrm{4}{\mathrm{d}}_{\mathrm{3}\mathit{/}\mathrm{2}}{\mathrm{1}}^{}{\mathrm{P}}_{\mathrm{1}}$–$\mathrm{1}{\mathrm{s}}^{\mathrm{2}}\hspace{0.17em}\mathrm{2}{\mathrm{s}}^{\mathrm{2}}\hspace{0.17em}\mathrm{2}{\mathrm{p}}^{\mathrm{6}}{\mathrm{1}}^{}{\mathrm{S}}_{\mathrm{0}}$ transition (commonly denoted as 4C) is located at 12.124 Å (Brown et al. 1998), that is, close to H-like Ne x, but we do not detect any other strong lines from iron, making an identification of the observed absorption feature with this line doubtful.

The presence of the Ne viii Kα line hints towards a complex interplay of absorption and emission features as also already seen in the low-hardness spectrum of the magnesium region (Sect. 3.2.2).

We tie the centers of the four emission features together and obtain a line velocity of $\mathrm{4}{\mathrm{0}}_{-50}^{\mathrm{+}\mathrm{80}}$ km s^-1. For Ne viii Kα, we obtain v = −220 ± 140 km s^-1, but given the low counts in the spectrum and the interplay of the emission and absorption, the formal errors are likely underestimated.

Because of the presence of H-like Ne x Lyβ in the low-hardness spectrum and suggestive residuals at ~10.2 Å, we add this line in absorption to our model and obtain an equivalent width of $\mathrm{-}\mathrm{0.02}{\mathrm{9}}_{-0.009}^{\mathrm{+}\mathrm{0.011}}$ Å and a line velocity of $\mathrm{-}\mathrm{38}{\mathrm{0}}_{-320}^{\mathrm{+}\mathrm{150}}$ km s^-1. This velocity is inconsistent with that of the emission lines and attempts to tie the velocities do not result in a satisfactory fit. The presence of Ne x Lyβ does not change the modeling results for other lines in the region.

For a single uniform medium, lines of the same series are expected to be all either in absorption or in emission. However, in the case of Vela X-1 we expect a complex geometry, as discussed in Sect. 4, and the sum of emission and absorption line contributions from different regions could result in the observed properties of the Ne x Ly series. The overall behavior however remains puzzling.

3.4. Iron Kα region

The region encompasses wavelengths in the range 1.6–2.5 Å, that is, the locations of the neutral Fe Kα (1.937 Å) and Fe Kβ (1.757 Å) lines as well as the iron edge (1.740 Å) and Ni Kα (1.66 Å) line. The use of a broader region than covered by the line features themselves allows better constraints on the power-law continuum. The results for this section are listed in Table 8 and shown in Fig. 11.

In the high-hardness spectrum, the Bayesian Blocks approach first finds the obviously present Fe Kα line, followed by the Fe Kβ and a step-like feature that we identify as the iron edge. We model the lines with Gaussian components and the edge with the edge model, fixing the edge energy to 1.740 Å (Fig. 11). The resulting line and edge parameters are listed in Table 8. Ni Kα has been previously detected in Vela X-1 (e.g., Martínez-Núñez et al. 2014). However, no feature at the location of the Ni Kα shortwards in wavelength from the edge is found with Bayesian Blocks; there is a strong nearby residual in the low-hardness spectrum, but given the low counts and large uncertainties of this region and the fact that the feature is only one bin wide, we do not consider it to be a physical Ni Kα line.

In the low-hardness spectrum, the Bayesian Blocks approach finds only the Fe Kα and Fe Kβ lines. The inclusion of the edge, however, still improves the fit and we obtain an absorption edge depth of 0.05 ± 0.04, that is, marginally inconsistent with zero.

4. Discussion

4.1. The source of varying absorption column density

Past analyses have established a relatively stable behavior of the central X-ray engine in Vela X-1 (Odaka et al. 2013) even during major flares (Martínez-Núñez et al. 2014). Our modeling in Sect. 2.2 suggests that the rapidly changing absorption is at least one of the sources of the observed changes of spectral shape but we cannot exclude changes in the underlying spectral shape as a further contribution. Physical origin of possible changes in spectral slope has been discussed, for example, by Odaka et al. (2014). Inhomogeneous absorption by clumpy matter has been used in Vela X-1 by Nagase et al. (1986) to explain the transient soft excess. Such absorption changes are expected either if the wind that the neutron star is accreting from is clumpy (review by Martínez-Núñez et al. 2017, and reference therein) or if the inner parts of the accretion flow are highly structured (Manousakis 2011; Manousakis & Walter 2015); both effects could be present simultaneously.

Clumpy winds are expected for O/B stars: line-driven winds of such stars are subject to instabilities (Lucy & Solomon 1970; Lucy & White 1980). Hydrodynamical simulations show that perturbations are present, grow rapidly and cause strong shocks so that dense gas-shells form (Owocki et al. 1988; Feldmeier et al. 1997). Density, velocity, or temperature variations compress the gas further and fragment it into small, overdense clumps embedded in tenuous hot gas (e.g., Oskinova et al. 2012; Sundqvist & Owocki 2013). In the HMXB Cyg X-1, for example, the clumpy wind paradigm has been used to explain the absorption dips, long-term variability of the absorption along the orbit, and the changing line content observed in high-resolution spectra during dips (Hanke et al. 2009; Grinberg et al. 2015; Miškovičová et al. 2016). In Vela X-1, clumpy winds can reconcile the discrepancies between different estimates of the mass-loss rate (Sako et al. 1999). Several lines of evidence also point towards the presence of wind clumping very close to the photosphere of O/B stars (Cohen et al. 2011; Sundqvist & Owocki 2013; Torrejón et al. 2015), well within the ~1.7 R_{HD 77581} distance between the neutron star and HD 77581 in Vela X-1.

To evaluate the impact of the clumpiness along the line-of-sight on the time variability of the column density, we designed a toy model consisting of a point source with a constant X-ray flux orbiting an O/B star at an orbital separation of 1.8 R_⋆ (Fig. 12). The wind is represented relying on the multi-dimensional numerical simulations of Sundqvist et al. (2017). They computed, as a function of time, the mass-density, velocity and pressure profiles in the wind up to 2 stellar radii. Doing so, they witnessed the formation of small overdense regions (clumps), embedded in a lower-density environment, and determined their dimensions and mass. We assumed an edge-on inclination and a self-similar expansion of their results to extend their simulation space to up to 15 R_⋆ and computed the column density along the line of sight at φ_orb = 0.21–0.25. Beyond the expected systematic decrease in column density as the compact object moves towards us, we observe a 30% peak-to-peak spread with a correlation time of at most 1 ks, a few times longer than the self-crossing time of the clumps. The combined absorption of independent clumps along the line-of-sight thus leads to larger and longer increases in column density than a single clump passing by, but not by a factor of a few in amplitude and not for a few kiloseconds. Therefore, given the characteristic size, speed, and distribution of clumps in the wind of isolated massive stars (Sundqvist et al. 2017), the inhomogeneities in a wind unperturbed by the presence of the compact object are not enough to account for the observed enhanced absorption (Sect. 2.2).

However, the compact object can significantly alter the wind structure, either through its gravitational (El Mellah & Casse 2017) or its radiative influence (Blondin et al. 1990; Manousakis & Walter 2015). The emerging structures (tidal stream, trailing accretion and photoionization wake) highly modulate the systematic orbital evolution of the column density (observational evidence in, e.g., Kaper et al. 1994; Fürst et al. 2010; Doroshenko et al. 2013; Malacaria et al. 2016). Yet, their influence should be minimal around φ_orb ≈ 0.25 since most of the aforementioned orbital-scale structures lie essentially along the axis joining the two body (star and neutron star) centers or trail the neutron star (Fig. 1). In this way, they are far enough from the line of sight that we do not expect much of the variation in the column density to be associated to them. Still, closer to the neutron star, Manousakis & Walter (2015) showed that even assuming a smooth upstream wind, the flow is highly structured and variable within the shocked region, possibly giving rise to changing absorption such as that observed here. Further numerical investigations are ongoing and aim at following the clumps as they encounter the shock (El Mellah et al. 2017): the role the clumps play in time-varying absorption would be altered by the deceleration and compression they undergo at the shock. There is also observational support for the presence of such shocks in Vela X-1 since its broadband X-ray count rate has been shown to follow a log-normal distribution (Fürst et al. 2010). This distribution can be described by a multiplicative stochastic model (Uttley et al. 2005; Fürst et al. 2010) and is known to be a possible result of shocks acting on clumpy material (see, e.g., Kevlahan & Pudritz 2009).

Further, newest calculations for stellar-atmosphere models estimate the wind velocity in the inner region, where the neutron star is located, to be as low as ~ 100 km s^-1 (Sander et al. 2017b), that is, on the order of the orbital velocity. In such a case, the wind and accretion structure could be significantly altered (Shapiro & Lightman 1976), possibly leading to the formation of a transient disk whose formation and dissipation could, for example, result in absorption and flaring events. A detailed modeling of such an accretion flow is beyond the scope of this paper, but we point out that, to the knowledge of the authors, no clear disk signatures have yet been detected in any of the X-ray observations of Vela X-1.

4.2. He-like triplet diagnostics: properties of the X-ray emitting gas

We measured the ratios R = f/i and G = (f + i) /r (Gabriel & Jordan 1969; Blumenthal et al. 1972; Porquet & Dubau 2000) for the He-like Mg xi and Ne ix triplets in the high-hardness spectra and obtain R = 1.5 ± 0.7 and $\mathit{G}\mathrm{=}\mathrm{2.}{\mathrm{0}}_{-0.7}^{\mathrm{+}\mathrm{3.5}}$ (Sect. 3.2.1) and $\mathit{R}\mathrm{=}\mathrm{0.}{\mathrm{8}}_{-0.6}^{\mathrm{+}\mathrm{1.4}}$ and $\mathit{G}\mathrm{=}\mathrm{4.}{\mathrm{3}}_{-2.9}^{\mathrm{+}\mathrm{5.7}}$ (Sect. 3.3.1), respectively. The presence of strong resonance lines (Figs. 8 and 10) and the measured best fit value of G< 4 for magnesium imply that the He-like triplets in the high-hardness spectrum may not be consistent with emission from a purely recombination-dominated photoionized plasma. While we do not use the Si xiii triplet to calculate G and R, the two measured lines in the high-hardness spectrum, Si xiii r and Si xiii f, are of comparable strengths, while the intercombination line is too weak to be detected (Figs. 6 and 7 and Table 4), supporting the results obtained for Mg and Ne. We do, however, emphasize the large uncertainties of the R and G values we obtain from our fits: the following discussion of plasma diagnostics should thus be seen as a roadmap for future investigations with missions with a larger collecting area such as JAXA’s X-ray Astronomy Recovery Mission (XARM) or ESA’s Athena (Nandra et al. 2013).

In principle, G< 4 (but note that G = 4 is still within the large measured uncertainties) could be due to the presence of a collisionally ionized component, which is predicted to have a much lower G ratio of ~1. However, it has been shown in a number of other accretion-powered photoionized plasmas that resonance scattering can also contribute to the emission, with the relative contributions of resonance scattering and recombination determining the observed G ratio (Sako et al. 2000; Kinkhabwala et al. 2002, 2003; Wojdowski et al. 2003). Kinkhabwala et al. (2002) calculated spectra for different ionic column densities in the ionization cone of a Type 2 Seyfert galaxy (i.e., viewed from the side, with direct emission from the central source obscured); based on Fig. 7 of their article, we can estimate that an ionic column density of order 10¹⁸ cm² is consistent with the observed G ratios of Ne ix and Mg xi.

To estimate the elemental column densities in the wind, we first calculate the characteristic hydrogen column density of the wind, N_⋆, that is defined as${\mathit{N}}_{\mathit{\star }}\mathrm{\equiv }\frac{\mathit{Ṁ}}{\mathrm{4}\mathit{\pi \mu }{\mathit{m}}_{\mathrm{p}}{\mathit{R}}_{\mathit{\star }}{\mathit{v}}_{\mathrm{\infty }}}\mathit{,}$(1)with Ṁ being the mass-loss rate, μ the mean weight, m_p the proton mass, R_⋆ the radius of the star and v_∞ the terminal velocity of the wind. This is the hydrogen column density that would be observed along a line of sight from the observer to the stellar surface in the case of constant wind velocity, and is a good order-of-magnitude proxy for the typical column density experienced by X-rays escaping along different lines of sight from the neutron star. For the typical values for HDE 77581, that is, Ṁ ~ 2 × 10^-6M_⊙ yr^-1, v_∞ = 750 km s^-1, and R_⋆ = 30 R_⊙, we obtain N_⋆ = 2.4 × 10²² cm^-2. Using the abundances of Wilms et al. (2000), we then estimate typical elemental column densities of 2 × 10¹⁸ for Ne ix, 6 × 10¹⁷ for Mg xi, and 4.5 × 10¹⁷ for Si xiii. These crude estimates for the elemental column densities are of order comparable to the above estimate of the ionic column densities from the G ratios, implying that a substantial fraction of the material along the line of sight is ionized to the He-like charge state for each of these elements. Further, more stringent observational constrains on G ratio at φ_orb ≈ 0.25 will be crucial to confirm or refute this line of argumentation that can currently be cast in doubt by the large measured uncertainties on G.

The measured values of R are not consistent with expectations for a low-density photoionized plasma. This is particularly striking for neon, where we would expect R ≈ 3.1 for a low-density plasma with T = 2 × 10⁵ K and higher R for higher temperatures (Porquet & Dubau 2000), all values above the measured $\mathit{R}\mathrm{=}\mathrm{0.}{\mathrm{8}}_{-0.6}^{\mathrm{+}\mathrm{1.4}}$. One possible interpretation is that the plasma density is high; however, the implied densities of 10¹²–10¹³ cm^-3 (Porquet & Dubau 2000, Fig. 9) are much higher than expected anywhere in the wind, or even in the accretion wake, and should only be found in the hydrostatic part of the stellar atmosphere. Possible structures close to the neutron star such as discussed in Sect. 4.1 are unlikely to lead to such extreme density enhancements; typical high densities reached in the simulations of Manousakis (2011) or Mauche et al. (2008) are of the order of a few 10¹⁰ cm^-3 and below. On the other hand, it is well known from studies of single OB stars (Kahn et al. 2001) that the strong stellar UV field of HDE 77581 will reduce the forbidden line strength near the star by depopulating the metastable $\mathrm{1}\mathrm{s}\mathrm{2}\mathrm{s}\hspace{0.17em}{\mathrm{3}}^{}{\mathrm{S}}_{\mathrm{1}}$ level (Gabriel & Jordan 1969; Blumenthal et al. 1972; Mewe & Schrijver 1978; Porquet et al. 2001).

In Fig. 13, we show predictions for the R ratio as a function of stellar radius for a B-type star with stellar parameters similar to HDE 77581 (T_eff = 25 000 K, log g = 3.0), using a methodology similar to that described in Sect. 2 of Leutenegger et al. (2006) to average over the relevant part of the TLUSTY model stellar continuum, that is, a non-LTE, plane-parallel, hydrostatic model stellar atmosphere (Lanz & Hubeny 2007), to estimate the UV field. In the context of the predicted ratios, and assuming a single radius of formation, the measurements would imply a formation radius of at least a few stellar radii, thus ruling out a substantial contribution of emission from the undisturbed stellar wind near the neutron star at r ≈ 1.7 R_{HD 77581} for Ne ix and Mg xi.

4.3. Presence of multi-component medium

In eclipse (φ_orb ≈ 0), co-existing emission lines from high- and low- ionization species of several elements have previously been detected by Sako et al. (1999) and Schulz et al. (2002) and in both cases interpreted as signatures of a multiphase medium, with cool dense clumps co-existing with hotter plasma. Similarly, the simultaneous presence of emission lines from high- and low-ionized ions of the same element at orbital phase 0.5, when the source is strongly absorbed by the accretion wake (Fig. 1), implies multiple gas phases with different temperatures in the medium (Goldstein et al. 2004; Watanabe et al. 2006).

Here, at φ_orb ≈ 0.25, we detect Kα transitions of L-shell ions of both silicon in emission in the high-hardness spectrum, that is, at higher overall absorption (Sect. 3.1.1), and of magnesium in absorption in the low-hardness spectrum, that is, at lower overall absorption (Sect. 3.2.2). In both cases, the L-shell ions co-exist with H- and He-like ions from the same elements. In the case of magnesium, we detect absorption from Mg xi r and Kα lines from Li-like Mg x, Be-like Mg ix, N-like Mg vi, O-like Mg v and, tentatively, C-like Mg vii. These lines imply the presence of cooler material in the system outside of the large-scale accretion structure, which should not be crossing our line of sight at φ_orb ≈ 0.25 (Fig. 1), as the comparatively modest measured absorption also implies. Clumps in the structured wind would provide a natural explanation for such material but so would a structured flow close to the neutron star. While detailed theoretical considerations of the interaction of a clumpy wind with a compact object and its radiation are beyond the scope of this work, we undertake two simplified estimates of whether or not the observed signatures can be produced in a single-zone absorber assuming a photoionized wind (Sect. 4.3.1) and whether or not the observed L-shell transitions may come from ions present in the stellar wind itself (Sect. 4.3.2).

Further hints of the presence of multiple components come from the detection of Ne x Lyβ in absorption and Ne x Lyα in emission in the low-hardness and thus low-absorption spectrum (Sect. 3.3.2). Lines of the same element in both emission and absorption cannot be produced in a single uniform medium and the observation can only be explained through the presence of multiple components whose sum results in the observed spectrum. Both, Ne x Lyα and Lyβ are detected in emission in the high-hardness and thus high-absorption spectrum (Sect. 3.3.1).

4.3.1. Photoionization simulations

First, to simulate a photoionized wind medium, we approximate the broadband spectrum of the irradiating continuum as the sum of the emission from the neutron star and from the companion. Since our data do not extend above 10 keV, we rely on the broadband (~3–70 keV) fits to NuSTAR observations obtained by Fürst et al. (2014) for the description of the neutron star continuum. They analyze two observations and find that a power-law model with a Fermi-Dirac cut-off (FDCUT; Tanaka 1986) gives a good description of the data: $\mathit{F}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{=}\mathit{A}{\mathit{E}}^{\mathrm{-}\mathrm{\Gamma }}{\left(\mathrm{1}\mathrm{+}\exp \left(\frac{\mathit{E}\mathrm{-}{\mathit{E}}_{\mathrm{cut}}}{{\mathit{E}}_{\mathrm{fold}}}\right)\right)}^{-1}\mathit{,}$(2)where A is the overall normalization, Γ the photon index, E_cut the cut-off energy , and E_fold the folding energy. We use Γ = 0.99, consistent with our measurement in Sect. 2.2 and with the average Γ obtained by Fürst et al. (2014) for their two observations, and set the normalization to 0.277 keV^-1 cm^-2 s^-1 (at 1 keV), the value we measure. Using their average values, we also set E_cut = 20.8 keV and E_fold = 11.9 keV. We model the contribution of HD 77581, which dominates the UV-continuum, as a blackbody using the parameters from Nagase et al. (1986): distance of 2 kpc, stellar radius of 31 R_⊙ and temperature of 25 000 K.

For our purposes here, it is sufficient to determine whether or not we can produce the observed range of Kα absorption lines of magnesium ions, from Mg v to Mg xi, in a single-zone absorber. With the above continuum model, we use XSTAR v2.31⁵ to calculate the ionization balance in the illuminated material for three different densities: 10⁴,10¹⁰, and 10¹⁶ cm^-3. Using line optical depths from XSTAR, we then generate approximate synthetic spectra from 8 to 10 Å, at the ionization parameters⁶ξ₁, where Mg v and Mg xi have comparable ion fractions, and ξ₂, where Mg v peaks. Regardless of density, at ξ₁ we overpredict the strength of the Mg xii line at 8.42 Å, and significantly underpredict the Mg v line. In general, at both ξ₁ and ξ₂, we expect lines from Mg vi to Mg viii that are comparable to or stronger than our Mg v/Mg xi lines of interest.

In short, it is difficult to produce both Mg v and Mg xi in a single absorption zone without producing other strong Mg lines between 8 and 10 Å. While such lines are not negligible in our data, they are generally weaker than we would expect based on our XSTAR calculations; further XSTAR runs show that the discrepancy between predictions and observations is exacerbated if we suppose that the absorber is shielded from the stellar continuum and only illuminated by the neutron star. Thus, although we have not performed any detailed ionization modeling of the absorber, it is unlikely that the observed lines originate in a single ionization zone. For the continuum spectra and ionization parameters considered here, the gas appears to be thermally stable (Krolik et al. 1981), indicating that multiple distinct absorbers (rather than a multi-phase gas) may be required.

4.3.2. Stellar atmosphere modeling

In a second approach, we use a detailed wind model for the donor star to see which ions are to be expected and how this changes due to the presence of a strong X-ray source. We employ the hydrodynamically consistent stellar atmosphere model calculated for Vela X-1 (Sander et al. 2017b) using the most recent version of the Potsdam Wolf-Rayet (PoWR) model atmosphere code that is able to obtain the mass-loss rate and wind stratification self-consistently for a given set of stellar parameters (Sander et al. 2017a). This model takes into account density inhomogeneities in the wind in the form of optically thin clumps but does not include changes in the wind geometry due to the presence of the neutron star.

For the purpose of our work, the donor star models from Sander et al. (2017b) were recalculated to account for higher Mg ions up to Mg viii and we refer the reader to this work for the details of the calculations. In these models, Si is covered up to Si vi. While we do see higher ions in our X-ray observations, they would not be expected from a single massive star wind (Fig. 14, blue curves). For the model without X-rays, only lower ions of both elements are found. For silicon, Si iv is leading in the wind until Si iii takes over in the outer part. Close to the donor star, Si v also plays a role, but already Si vi is essentially not populated in the wind. For magnesium, the situation is even simpler with Mg iii being the leading ion in the wind and all higher ions quickly becoming depopulated beyond the surface of the donor.

The situation changes when X-rays are included. Although the treatment of the X-ray irraditaion in the atmosphere models is only a rough approximation, as these are 1D models and the X-rays are modeled as Bremsstrahlung, the changes in the population numbers illustrate that the X-ray illumination due to the accretion onto the neutron star leads to a more complex structure and to the presence of several higher ions. Si vi and Si v are now only about an order of magnitude less populated than the leading Si iv; they remain important until the wind is essentially transparent and Si vi becomes the leading ion in the outer wind (Fig. 14, left panel). For magnesium, the situation is roughly similar, but more complex as more ions are present (Fig. 14, right panel). At the distance of the neutron star, the previously unimportant ions Mg iv to Mg vii are now almost similarly populated. At larger distances, Mg vii becomes the leading ion and only Mg vi remains significantly populated. Interestingly, the relative population of Mg viii increases outwards, but then drops down again. However, the later could be an artifact from our restriction to 1D models with monotonic wind velocity laws.

In spite of the limitations of our models, we can conclude that the ionization situation around Vela X-1 is complex. The PoWR models have demonstrated that the donor wind as such produces very low ionized metals (e.g., C iii, Mg iii, Si iv, Fe iv), but the illumination due to X-rays changes this situation significantly and even in the 1D approach, distinct regions with different leading ions are present. However, the model also demonstrated that some ionization stages can exist in the same region. Given the fact that due to the orbit of the neutron star only a part of the wind is strongly illuminated, while other parts are less or even not at all affected by X-rays, it is safe to conclude the observed X-ray spectrum stems from a multitude of layers with different ionization situations.

4.4. Iron fluorescence

The location of the detected Fe Kα emission features (Table 8) is consistent with the Kα₁ and Kα₂ lines of neutral or lowly ionized iron ions below Fe xii (House 1969). Interestingly, PoWR model atmosphere calculations for Vela X-1 (Sander et al. 2017b, see also Sect. 4.3.2) imply that the main iron ion species present in the wind without X-ray irradiation should be Fe iv, followed by Fe iii. In the case of strong X-ray irradiation, Fe v becomes more strongly populated and surpasses Fe iii. Furthermore, in this case, higher ionization stages (Fe vii and above) can also be reached in the outer wind.

The flux ratio between Fe Kβ and Fe Kα is consistent with the theoretical expectations of ~0.13–0.14 (Palmeri et al. 2003). We note, however, the large uncertainties on Fe Kβ which are exacerbated by the presence of the Fe K edge.

Kallman et al. (2004) predict the curve of growth for the equivalent width of Fe Kα, EQW(Fe Kα), in units of eV, and the hydrogen equivalent absorption column density of the reprocessing material, ${\mathit{N}}_{\mathrm{H}}^{\mathrm{22}}$, in units of 10²² cm^-2, to be $\mathrm{EQW}\mathrm{\left(}\mathrm{Fe}\mathrm{K}\mathit{\alpha }\mathrm{\right)}\backsimeq \mathrm{3}{\mathit{N}}_{\mathrm{H}}^{\mathrm{22}}$ (see also Inoue 1985). This has been observationally verified for HMXBs by Torrejón et al. (2010) using a Chandra-HETG sample that included Vela X-1. To compare our results, we calculate the equivalent width of the Fe Kα and Kβ lines in eV (Table 8) and use the values for N_H obtained from continuum fits (Sect. 2.2). The observed lines are much stronger than would be predicted by the above correlation; however, Torrejón et al. (2010) obtain similar results for the same data: in their analysis, Vela X-1 clearly lies above the correlation.

Several effects can contribute to the disagreement between the measured values and the theoretically predicted curve of growth: on one hand, N_H measurements are generally dependent on the continuum model and our choice of a power law to describe the 2.5–10 keV may lead to systematic deviations, although we note that our results roughly agree with the time averaged analyses of the same Chandraobservation by other authors (Goldstein et al. 2004; Torrejón et al. 2010). On the other hand, Kallman et al. (2004) compute the above relationship for a spherical geometry that is clearly not given in Vela X-1, given the influence of the compact object on the large-scale wind structure. The surface of the companion star can also contribute significantly to the observed FeKα emission as shown by Watanabe et al. (2006) who, however, also show that this cannot account for the full observed equivalent width. We further implicitly assume that the iron fluorescence happens in the same material that is absorbing the X-rays along the line of sight. For a clumpy absorber this does not have to be the case and the neutron star may be hidden behind some material without the ambient source of fluorescent emission being affected (Inoue 1985). A cool cloud, for example the accretion wake, irradiated by the neutron star could account for the excess emission (Watanabe et al. 2006). Three reprocessing sites – stellar wind, companion atmosphere, and cold material in the proximity of the neutron star – have also been suggested by Sato et al. (1986). We note that such a geometry would imply a more complex continuum model, with multiple differently absorbed continua as has been employed, for example, by Martínez-Núñez et al. (2014) to describe XMM-Newton-pn observations of Vela X-1, but is not useable in our case due to data quality.

5. Summary and outlook

We have shown that at φ_orb ≈ 0.25, when our line of sight to the compact object is expected to be least affected by the large-scale accretion structure, Vela X-1 shows spectral changes on timescales of ~hours (~a few ks). These changes are at least partially driven by fluctuating absorption. The variable absorption cannot be explained by unperturbed clumpy stellar winds, but the large-scale accretion structure (accretion wake, photoionization wake, accretion stream) is also not expected to vary on such short timescales or to cross our line of sight at this orbital phase at all.

We have detected variable emission and absorption line features from neon, magnesium and silicon at φ_orb ≈ 0.25, that is, at an orbital phase when few line features have been reported previously. Overall, the observational picture is consistent with a complex geometry, with multiple phases co-existing in the system and contributing different components to the overall observed spectrum. He-like triplet diagnostics imply that a substantial fraction of the material along the line of sight further away from the compact object is ionized to He-like charge state of magnesium and neon. The L-shell ions of Si and Mg are a signature of the presence of cooler material. The strong Fe Kα implies the presence of cool material as source for the fluorescent emission.

We speculate that dense regions could arise in shocks close to the compact object and have longer crossing times and higher equivalent hydrogen column density than unperturbed wind clumps, thus explaining the observed variability. The possible origin of such dense regions through the interaction of wind clumps at shocks is a notion to be explored in upcoming simulations of wind accretion in HMXBs that can then be compared to the observations presented here.

We show that time-resolved spectroscopy is crucial to disentangle the contributions of different components. However, the presented Chandra-HETG observation does not allow us, for example, to explore possible fluctuations in absorption during the high-hardness periods, to perform pulse-resolved spectroscopy of the line features, or to measure the possible response of the ionized gas to changes in the X-ray flux of the neutron star. More observations at φ_orb ≈ 0.25 will allow to accumulate more data and provide insight into the finer structure of the structured absorber. The major leap forward, however, will only be possible with future instruments that have a higher effective area. Future work would also profit from an improvement in line reference values of L-shell ions for magnesium, both to enable a better analysis of the currently available high-resolution data and for forthcoming missions with their higher throughput and/or higher resolution.

In particular, Vela X-1 was among the planned targets for Astro-H/Hitomi (Kitamoto et al. 2014) and is thus also likely to be among the main sources observed in the future by approved microcalorimeter missions such as Athena (Nandra et al. 2013) and XARM⁷ as well as planned high-resolution spectroscopy grating missions such as Arcus (Smith et al. 2014), a Medium-Class Explorer mission that has recently been selected for a concept study, and X-ray Surveyor/LYNX⁸. Their higher sensitivity will allow a detailed study of changing absorption – but our results also imply that we will need to carefully consider ks-timescale variability when drawing any conclusions from these future observations.

Note added in proof. As pointed out by K. Postnov (priv. comm.), an alternative explanation for the observed timescales is subsonic accretion onto the NS magnetosphere by forming a convective plasma shell, as suggested in Shakura et al. (2012). In this case the characteristic variability time is ${\mathit{R}}_{\mathrm{B}}^{\mathrm{3}\mathit{/}\mathrm{2}}$/$\sqrt{\mathit{GM}}$, where R_B is the Bondi radius, which is about an hour for an assumed stellar wind velocity of 500 km s^-1 in Vela X-1. A quantitative estimation of absorption variations inside such a shell is, however, not easily possible.

---

Acknowledgments

V.G. is thankful for support of her work through the ESA internal research fellowship. Part of this work was performed under the auspices of the US Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. Support for this work was provided by NASA through the Smithsonian Astrophysical Observatory (SAO) contract SV3-73016 to MIT for Support of the Chandra X-Ray Center (CXC) and Science Instruments; CXC is operated by SAO for and on behalf of NASA under contract NAS8-03060. I.E.M. acknowledges funding from the European Union’s Horizon 2020 research and innovation programme under the Marie SkłlCodowska-Curie grant agreement No. 665501 with the Research Foundation Flanders (FWO). A.A.C.S. is supported by the Deutsche Forschungsgemeinschaft (DFG) under grant HA 1455/26. J.N. acknowledges support from NASA through the Hubble Postdoctoral Fellowship Program, grant HST-HF2-51343.001-A. S.M.N. acknowledges support by research project ESP2016-76683-C3-1-R. For the initial data exploration, this research has made use of the ChandraTransmission Grating Data Catalog and Archive⁹ (tgcat; Huenemoerder et al. 2011). This research has made use of ISIS functions (isisscripts)¹⁰ provided by ECAP/Remeis observatory and MIT and of NASA’s Astrophysics Data System Bibliographic Service (ADS). We thank John E. Davis for the development of the slxfig¹¹ module used to prepare most of the figures in this work. Some of the color schemes used were based on Paul Tol’s palettes and templates¹². I.E.M., A.A.C.S., F.F., P.K., M.K., S.M.N. and J.W. are grateful for the hospitality of the International Space Science Institute (ISSI), Bern, Switzerland which sponsored a team meeting initiating a tighter collaboration between stellar wind and X-ray specialists. And finally, we thank the anonymous referee whose comments and suggestions helped us to improve this manuscript.

References

All Tables

All Figures