© The Authors 2025

1. Introduction

High-mass X-ray binaries (HMXBs) are a sub-class of X-ray binaries in which a compact object accretes matter from a massive (> 10 M_⊙) main-sequence stellar companion predominantly via the capture of stellar wind. 4U 0900−40 (henceforth Vela X−1) was discovered with the Uhuru satellite in 1967 (Chodil et al. 1967). It is the prototypical super-giant high-mass X-ray binary (sg-HMXB) pulsar. Vela X−1 comprises a neutron star with a mass of M_NS ∼ 1.77 M_⊙ (Quaintrell et al. 2003) and a B0.5Ib optical companion that has a mass of M_B ∼ 23 M_⊙ and a radius of R_* ∼ 30 R_⊙ (van Kerkwijk et al. 1995a). The neutron star has a spin period of P_s ∼ 283.5 s (Rappaport 1975). It is an eclipsing HMXB with an orbital period of ∼9 days (van Kerkwijk et al. 1995b). The distance to the source has been calculated to be $1.{99}_{-0.11}^{+0.13}$ kpc (Kretschmar et al. 2021). The neutron star of Vela X−1 orbits the companion in a slightly eccentric ∼0.1 (Bildsten et al. 1997) orbit with an apastron distance of ∼50 R_⊙ (Bildsten et al. 1997) from the companion.

The companion star in the Vela X−1 binary system is known to have a mass loss rate of ∼10⁻⁶ M_⊙ yr⁻¹ (Nagase et al. 1986) with a wind terminal velocity of v_∞ ∼ 1100 km s⁻¹ (Watanabe et al. 2006). The neutron star accretes some of the stellar wind material, leading to X-ray emission. The binary environment populated by the stellar wind of the companion reprocesses part of the emitted X-rays via scattering and fluorescence emission. The signatures of the reprocessing region are encoded in the X-rays that reach us in the form of prominent Iron Kα emission lines (Becker et al. 1978; Goldstein et al. 2004; Watanabe et al. 2006; Aftab et al. 2019; Rikame et al. 2024). Some of the emitted X-rays are photo-electrically absorbed by the material crossing our line of sight (Verner et al. 1996). The variation in the absorption column density in the binary orbit tells us about the distribution of the absorbing medium. Vela X−1 exhibits variability in both absorption (Nagase et al. 1986; Haberl & White 1990; Watanabe et al. 2006; Malacaria et al. 2016; Kretschmar et al. 2021) and overall flux (Kreykenbohm et al. 2008; Fürst et al. 2010; Kretschmar et al. 2021) on different timescales, from months (Haberl & White 1990; Kreykenbohm et al. 2008) down to seconds (Fig. 1a of Inoue et al. 1984, Fig. 4 of Boerner et al. 1987, Fig. 5 in Kreykenbohm et al. 2008). Vela X−1 is also known to have off-states occurring outside of eclipses with minimal or undetectable X-ray emissions. These off-states have been attributed to diverse factors such as low-density regions or accretion processes hindered by magnetospheric effects (Inoue et al. 1984; Kreykenbohm et al. 1999; Sidoli et al. 2015; Kretschmar et al. 2021). Several studies have revealed that the stellar wind of the companion in Vela X−1 is not symmetric (Blondin et al. 1990; Doroshenko et al. 2013; Diez et al. 2022, 2023). Hence, the smooth-stellar wind model (Castor et al. 1975) cannot explain the variation in the absorption column density with the orbital phase. Such an effect was also seen in some other HMXBs, such as IGR J17252−3616 (Manousakis et al. 2012) and 4U 1538−52 (Rodes-Roca et al. 2015). It has been suggested that this effect is due to the formation of a photoionisation wake, an accretion wake, or both because of the influence of a neutron star on the stellar wind. A photoionisation wake is formed when the X-rays emitted from the neutron star slow down the ejected wind by photoionising the medium, which reduces the scattering cross-section in the region; hence, a dense region within the strömgren sphere trailing the neutron star is formed (Blondin et al. 1990; Kaper et al. 1994). An accretion wake is formed as the gravity of the neutron star focuses the stellar wind onto its centre of gravity, thereby creating an extra-dense structure that trails the neutron star. Recently, when our work was under review, Abalo et al. (2024) reported results from their study of orbit-to-orbit variability in the flux of Vela X−1 using more than 14 years of MAXI/GSC data. They used the hardness ratio to quantify the X-ray-absorbing material and concluded that a combination of an accretion wake, a photoionisation wake, and the clumpy nature of the stellar wind, along with the inherent variability of the source, contributes to these variations stochastically.

To explain the asymmetric variation in the absorption column density in Vela X−1, Doroshenko et al. (2013) proposed a toy model to include the contribution of the photoionisation wake. In this work, we have studied how the asymmetric variation changes across three different intensity levels – low, mid, and high – by performing intensity-and-orbital-phase resolved spectroscopy. For this analysis, we have used the long-term data acquired with MAXI/GSC in the last ∼13 years. We have also estimated the wind parameters governing the stellar wind with the absorption column density obtained from the orbital-phase resolved spectroscopy with a similar wake model.

The paper is organised as follows. In Sect. 2, we present the instrument and data analysis techniques used. The timing analysis of the overall light curve (Sect. 3.1) and the intensity-resolved light curves (Sect. 3.2) are presented in Sect. 3. In Sect. 4, we perform time-averaged spectroscopy (Sect. 4.1), orbital-phase resolved spectroscopy (Sect. 4.2), intensity-resolved spectroscopy (Sect. 4.3), and intensity-and-orbital-phase resolved spectroscopy (Sect. 4.4). In Sect. 5, we describe a toy model to explain the variation in the observed absorption column density. Finally, the implications of the results are discussed in Sect. 6.

2. Instrument and observation

Monitor of All-Sky X-ray Image (MAXI) is an X-ray monitoring mission operated from the International Space Station (Matsuoka et al. 2009). The MAXI has two X-ray instruments; namely, the Gas Slit Camera (GSC) (Mihara et al. 2011) and the Solid-state Slit Camera (SSC) (Tomida et al. 2011). The MAXI/GSC instrument consists of 12 position-sensitive gas proportional counters that operate in the energy range 2.0−30.0 keV, and MAXI/SSC has charge-coupled devices (CCDs) that operate in the energy range 0.5−12 keV. These two instruments cover the entire sky in every ISS orbital period (∼90 minutes) with fields of view of 1$\stackrel{°}{.}$5 × 160° and 1$\stackrel{°}{.}$5 × 90° for GSC and SSC, respectively.

For the current work, we have utilised MAXI/GSC data for Vela X−1 from MJD 55200 to MJD 59791. The spectra and light curves were obtained from the MAXI on-demand website¹. The timing resolution of all the light curves was taken as 0.0625 d (∼90 min ∼one ISS orbit). For orbital-phase resolved and intensity-resolved spectral analysis, we downloaded the data using good time interval files to select the relevant phase and intensity level, respectively.

3. Timing analysis

3.1. Overall light curve

The long-term light curve of Vela X−1 obtained in the 2.0−20.0 keV range is plotted in Fig. 1 with a bin size of one orbital period. The light curve shows long-term variations (Fig. 1) with no apparent periodicity. The orbital period of Vela X−1 was determined from the MAXI/GSC light curve using the epoch folding tool, efsearch² to be 8.96432 ± 0.00029 days, consistent with Kreykenbohm et al. (2008). The other orbital parameters used in this paper are given in Table 1. The orbital intensity profile in 2.0−20.0 keV energy band is plotted in brown in Fig. 2. The orbital profile of Vela X−1 is asymmetric in nature, with a significantly higher count rate before an orbital phase of ∼0.4 (Doroshenko et al. 2013).

Fig. 1. ∼13-year long-term MAXI/GSC light curve of Vela X−1 with a bin size of one orbital period. The red, green, and blue arrows mark the low (< 0.36 counts s−1), mid (0.36−0.46 counts s−1), and high (> 0.46 counts s−1) intensity levels. The horizontal lines (cyan) separate the three intensity levels.

To study the variation in the orbital profile with energy, we downloaded the source light curves in soft and hard energy bands of 2.0−4.0 keV (soft) and 12.0−20.0 keV (hard), respectively. The soft and hard energy bands orbital profiles are plotted in Fig. 2 with 64 bins per period, and the hardness ratio between the two bands is also shown. The count rate in the soft energy band decreases after an orbital phase of 0.4, whereas the hard X-ray profile is nearly flat. The hardness ratio also shows higher values during eclipse-ingress and eclipse-egress. There is an excess of hard photons after an orbital phase of 0.4. A similar variation in the hardness ratio was reported by Diez et al. (2022, 2023) with NuSTAR (2022) and XMM-Newton (2023) observations. A similar variation in the averaged hardness ratio between the energy bands 2.0−4.0 keV and 4.0−10.0 keV was shown by Abalo et al. (2024), even though the hardness ratio varies from orbit-to-orbit.

Fig. 2. Energy-dependent orbital profiles 2.0−4.0 keV (purple), 12.0−20.0 keV (olive), and 2.0−20.0 keV (brown) are plotted with the bin size of 64 bins per interval. The hardness ratio (black) is plotted between 2.0−4.0 keV (soft) and 12.0−20.0 keV (hard) energy bands. The vertical lines (cyan) divide the orbital profiles into eclipse, egress, OOE, and ingress.

3.2. Intensity resolved light curves

The long-term light curve of Vela X−1 in Fig. 1 suggests that the source intensity varies from orbit to orbit. The long-term MAXI/GSC data has covered approximately 512 Vela X−1 binary orbits during the observation period. We are interested in intensity variation in the out-of-eclipse (OOE) region; hence, we selected approximately the phase range 0.15 − 0.80 (where ϕ = 0 is the mid-eclipse) to study the orbit-to-orbit intensity variation. The average count rate per orbit for all orbits in the aforementioned phase range was calculated, and its histogram is plotted in Fig. 3. The distribution suggests that the average intensity indeed varies by a few factors from orbit to orbit; the variable nature of Vela X−1 over different timescales has been reported in the literature (Kreykenbohm et al. 2008; Fürst et al. 2010; Malacaria et al. 2016). The orbits of Vela X−1 were divided into three intensity levels. While dividing, we ensured that the sum of the averaged count rate of all the orbits in each intensity region was comparable. In the present case, the sum is around 65 counts. Under this condition, the following threshold values for count rates were obtained for the different intensity levels: low, < 0.36 counts s⁻¹ with 222 orbits; mid, 0.36 − 0.46 counts s⁻¹ with 161 orbits; and high, > 0.46 counts s⁻¹ with 111 orbits, respectively. The orbital profiles in the energy band of 2.0−20.0 keV in the low, mid, and high-intensity levels are plotted in Fig. 4 with 32 bins per period. The orbital intensity profile in the high-intensity level has some sub-structures. It is not known what causes the multiple peaks in the high-intensity level, but such peaks are also observed in the orbital intensity profile of Vela X−1 made in the hard X-ray band using the INTEGRAL/ISGRI (17−150 keV) light curve (Fig. 3 of Falanga et al. 2015). In the mid-state, there is a dip between phase 0.4−0.5 on top of a decaying profile. Meanwhile, in the low state, there is a gradual decay of the count rate. The intensity resolved hardness ratio plot is given in Fig. 5; the hardness ratio between 2.0−4.0 keV and 12.0−20.0 keV in the low-intensity level rises slowly compared to the high-intensity level.

Fig. 3. Histogram of total intensity per orbit grouped into three intensity levels. The X axis represents the MAXI/GSC counts rate per orbit. The Y axis represents the number of orbits. The counts rate range for the low (red), mid (green), and high (blue) intensity levels are < 0.36 counts s−1, 0.36 − 0.46 counts s−1, and > 0.46 counts s−1, respectively. The total number of orbits in the three intensity levels are 222, 161, and 111, respectively.

Fig. 4. 2.0−20.0 keV orbital profiles of Vela X−1 are plotted for low (red), mid (green), and high (blue) intensity levels with the bin size of 32 bins per interval. The vertical lines (cyan) divide the orbital profiles into eclipse, egress, OOE, and ingress.

Fig. 5. Hardness ratio plots between the soft (S; 2.0−4.0 keV) and hard (H; 12.0−20.0 keV) energy bands for the low (red), mid (green), and high (blue) intensity levels. The inset in the upper right corner shows that the hardness ratio in the low-intensity level rises slowly compared to the high-intensity level in the phase range 0.4−0.6.

4. Spectral analysis

The variation in the hardness ratio throughout the orbit and between different intensity levels indicates changes in the spectral parameters. We performed time-averaged spectroscopy and orbital-phase resolved spectroscopy to probe the orbital variation in spectral parameters in the overall spectrum and in different intensity levels. The errors on best-fit spectral parameters in Table 2 are quoted at a 90% confidence level. The spectral analysis was carried out using XSPEC v: 12.13.0c (Arnaud 1996). The elemental abundance was taken from Wilms et al. (2000) and the photoelectric cross-sections from Verner et al. (1996).

4.1. Time-averaged spectroscopy

Different phenomenological models have been used to study Vela X−1. In particular, some have described its continuum emission by a cut-off power law (e.g. Malacaria et al. 2016) or a power-law with a Fermi Dirac cut-off (e.g. Diez et al. 2022, 2023). In this work, we use a physically motivated Comptonised continuum model, which allows us to compare our results with the ones reported earlier by Doroshenko et al. (2013). Hence, our final best-fit model is tbabs*pcfabs(CompST + Gaussian), where CompST is the inverse-comptonisation of cool photons by the hot electrons (Sunyaev & Titarchuk 1980), tbabs (Wilms et al. 2000) models the photoelectric absorption by the Galactic ISM along the line of sight, pcfabs was used for incorporating the contribution from the variable local absorption by the stellar wind, and there is a prominent Fe Kα line at ∼6.4 keV that was modelled with a single Gaussian emission line (Gaussian). The Galactic absorption component modelled by tbabs was frozen to the Galactic value (∼0.37 × 10²² cm⁻²) along the line of sight (HI4PI Collaboration 2016). The best-fit parameter values along with error estimates at a 90% confidence level are given in Table 2 and are consistent with results reported by Doroshenko et al. (2013) within the error bars. The corresponding spectrum, best-fit model, and residuals for the best-fit model are shown in Fig. 6.

Fig. 6. Time-averaged spectrum of Vela X−1 in 2.0−20.0 keV energy band obtained from MAXI/GSC. Top panel: Data (black crosses) and best-fit model (dashed line) components. Bottom panel: Deviation of the data from the best-fitted model spectrum.

4.2. Orbital-phase resolved spectroscopy

In order to study the orbital variation in the hardness ratio (Figs. 2, 5) further, we performed orbital-phase resolved spectroscopy. We divided the Vela X−1 orbital profile (Fig. 2) into three regions, namely the eclipse-egress (egress hereafter), OOE, and eclipse-ingress (ingress hereafter), with corresponding phase ranges of 0.09−0.16, 0.16−0.83, and 0.83−0.91, respectively, after visual inspection of the hardness ratios and orbital profiles. The three regions were further divided into multiple orbital phase bins, with the OOE into 14 orbital phase bins, and the egress and ingress into two orbital phase bins each, ensuring that each phase bin contains approximately the same total number of source photons. Using orbital ephemeris (given in Table 1), we performed joint fitting of all 18 spectra with the model tbabs*(CompST + Gaussian). Here, we removed pcfabs, considering the reduced total counts for orbital-phase resolved spectroscopy, and thus avoided over-fitting the spectra, similar to what was done in Doroshenko et al. (2013). We let tbabs alone account for both local and Galactic absorption. We fixed the temperature (T) and optical depth (τ) of CompST and the line centre and width of the iron emission line to their respective time-averaged values (given in Table 2), as they could not be constrained in individual orbit phases due to limited photon statistics. The variation in different spectral parameters with the orbital phase is shown in Fig. 7. The top panel in Fig. 7 shows the variation in the absorption column density with the orbital phase. The absorption column density starts increasing after phase 0.4 as in the case of the hardness ratio. Such a variation in Vela X−1 was previously reported in the literature (Doroshenko et al. 2013; Diez et al. 2022, 2023).

Fig. 7. Variation in spectral parameters (black) obtained from orbital-phase resolved spectroscopy. The error bars are at the 90% confidence level. The orbital profile (brown) is plotted in the 2.0−20.0 keV energy band. The vertical lines (cyan) are drawn to denote the four different regions (eclipse, egress, OOE, and ingress).

4.3. Intensity resolved spectroscopy

The asymmetric distribution observed in the absorption column density depicted in Fig. 7 and varying rise of the hardness ratio across different intensity levels in Fig. 5 prompts us to investigate similar variations across different intensity levels. Additionally, the fluctuating count rates from orbit to orbit illustrated in Fig. 4 suggest potential differences in spectral parameters across different intensity levels.

We did intensity-resolved spectroscopy by fitting each spectrum to the same spectral model, tbabs*pcfabs*(CompST + Gaussian), as the time-averaged spectrum (Sect. 4.1). Similar to time-averaged spectroscopy, we fixed tbabsN_H to the Galactic absorption column density along the line of sight, while allowing the local absorption parameters (pcfabs) to vary. Due to constraints imposed by reduced photon counts, the emission line centre and width were fixed to their respective values determined from the time-averaged spectrum. The rest of the parameters were kept free. The best-fit parameters’ values obtained from the fits for three different intensity levels are given in Table 2 and the corresponding spectra are shown in Fig. 8. From Fig. 8, the three intensity states do not show significant changes in spectral shape, but the total emission strength has increased from the low- to the high-intensity level. The iron line is present in all three spectra. In Table 2, along with an increase in the intensity levels, normalisations of the continuum (CompST) and the iron emission line (Gaussian) are found to increase.

Fig. 8. Three intensity-resolved spectra in the 2.0−20.0 keV energy band.

4.4. Intensity-and-orbital-phase resolved spectroscopy

We also performed orbital-phase resolved spectroscopy on each intensity level to probe the variation in patterns in the spectral parameters as a function of the orbital phase. The net count rate was reduced in each intensity level; hence, the orbital OOE region was subdivided into six spectral bins, along with one bin each for egress and ingress. We ensured an approximately equal number of counts for all phase bins. We conducted joint fitting of eight spectra of the three intensity levels using the same spectral model utilised for orbital-phase resolved spectroscopy described in Sect. 4.2. The temperature (T) and optical depth (τ) of CompST, as well as the centre and width of iron emission lines, were fixed to their respective intensity resolved values. The results, depicted in Fig. 9, illustrate how spectral parameters vary with the orbital phase across the three intensity levels. Notably, the asymmetric variation in the absorption column density persists across all three intensity levels. There is a distinct shift in absorption column density around the orbital phase ∼0.4. As can be seen in the top panel of Fig. 9, the increase in the absorption column density happens earlier as we move from low to high intensity levels.

Fig. 9. Variation in spectral parameters obtained from intensity-and-orbital-phase resolved spectroscopy for the low (red), mid (green), and high (blue) intensity levels. The error bars are at a 90% confidence level. The orbital profiles for low (red), mid (green), and high (blue) intensity levels are plotted in the 2.0−20.0 keV band.

5. Wind modelling

The asymmetric distribution of the absorption column density after phase 0.4 in both orbital-phase resolved (Fig. 7) and intensity-and-orbital-phase resolved (Fig. 9) spectroscopy is indicative of the neutron star getting obscured behind dense absorbing matter after that phase. The variation in the overall hardness ratio in Fig. 2 also suggests the same. In the following sections, we try to understand how the changing geometry of the absorbing medium composed of stellar wind and a trailing photoionisation wake can lead to the observed asymmetric variation in the absorption column density at different intensity levels.

5.1. Smooth stellar wind with photoionisation wake: A toy model

The wind velocity profile, v^SSW(r), of the radiatively driven smooth-stellar wind (Castor et al. 1975) at a radial distance, r, from the centre of the companion star is given as

$$\begin{array}{c}\hfill v^{\mathrm{SSW}}(r)=v_{\infty }(1-\frac{R_{\ast }}{r}{)}^{{\beta }^{\mathrm{SSW}}},\end{array}$$(1)

where, v_∞ is the terminal velocity of the stellar wind, β^SSW is the velocity exponent, and R_* is the radius of the companion star. The hydrogen number density for the smooth-stellar wind, $n_{\mathrm{H}}^{\mathrm{SSW}}(r)$, is given as

$$\begin{array}{c}\hfill n_{\mathrm{H}}^{\mathrm{SSW}}(r)=\frac{{\chi }_{\mathrm{H}}\dot{M}}{4\pi r^2{v}_{\mathrm{SSW}}(r)},\end{array}$$(2)

where χ_H is the conversion factor from total mass density to effective hydrogen particle density and Ṁ is the mass loss rate from the companion star. The absorption column density, ($N_{\mathrm{H}}^{\mathrm{SSW}}$), was obtained by integrating Eq. (2) along the line of sight from the position of the neutron star (r′_NS) in the orbit to the observer as

$$\begin{array}{c}\hfill N_{\mathrm{H}}^{\mathrm{SSW}}={\displaystyle {\int }_{r_{\mathrm{NS}}^{\prime }}^{\infty }}{n}_{\mathrm{H}}^{\mathrm{SSW}}(r^{\prime })\mathrm{d}{r}^{\prime }.\end{array}$$(3)

To construct the orbit of Vela X−1, we used the orbital ephemeris given in Table 1. The inclination of the binary plane is i > 74° (Table 1), but for simplicity we assumed the binary plane to be edge-on; that is, i = 90° (similar to Doroshenko et al. 2013).

Symmetric stellar wind gives rise to a symmetric variation in the absorption column density. Therefore, to account for the asymmetric variation in the absorption column density (top panel of Figs. 7 and 9), we adopted a simple toy model similar to the one of Doroshenko et al. (2013), consisting of a photoionisation wake formed as a result of intense X-ray emission of the neutron star (Fig. 10). In our model, the photoionisation wake is a cone with an opening angle (Θ). The axis of the cone is bent away from the direction of motion of the neutron star by a constant factor called the phase-lag parameter (ϵ) so that it trails the neutron star. Boundary 1 (bd1) and 2 (bd2) of the wake are described by the following co-ordinates:

Fig. 10. Sketch of Vela X−1 with the photoionisation wake (in grey). NS (in black) is the neutron star at phase 0.4. The companion (in blue) is at the focus of the orbit. The angles θ and Θ are the immersion and the opening angles of the neutron star and the wake, respectively. The line AB is the line of apsides, and the line of sight is along the positive X axis.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill x_{\mathrm{wake}}^{\mathrm{bd}1},y_{\mathrm{wake}}^{\mathrm{bd}1}& =rsin(ϵr),rcos(ϵr)\hfill \\ \hfill x_{\mathrm{wake}}^{\mathrm{bd}2},y_{\mathrm{wake}}^{\mathrm{bd}2}& =rsin(\mathrm{\Theta }+ϵr),rcos(\mathrm{\Theta }+ϵr).\hfill \end{array}\end{array}$$(4)

Within the wake, we reduced the velocity of the stellar wind by a constant factor called the velocity reduction factor (f) to make the stellar wind denser than the region outside the wake. We also assumed a different velocity exponent (β^wake) for the wake. The neutron star sits fully immersed in the wake, making an angle, θ, with respect to the leading edge of the wake, which we call the immersion angle (θ) (Fig. 10). Therefore, the velocity, v^wake, and the hydrogen number density, $n_{\mathrm{H}}^{\mathrm{wake}}(r)$, of the wake are given by

$$\begin{array}{c}\hfill v^{\mathrm{wake}}(r)=\frac{v_{\infty }}{f}(1-\frac{R_{\ast }}{r}{)}^{{\beta }^{\mathrm{wake}}}\end{array}$$(5)

and

$$\begin{array}{c}\hfill n_{\mathrm{H}}^{\mathrm{wake}}(r)=\frac{{\chi }_{\mathrm{H}}\dot{M}}{4\pi r^2{v}^{\mathrm{wake}}(r)}·\end{array}$$(6)

Both Eqs. (4) and (6) hold for a radial value of r < r₀, where r₀ is the length scale chosen at a distance where the hydrogen number density inside the wake becomes comparable to the hydrogen number density of smooth stellar wind.

The equivalent hydrogen number density, n_H(r), is given by

$$\begin{array}{c}\hfill n_{\mathrm{H}}(r)=n_{\mathrm{H}}^{\mathrm{SSW}}(r)+n_{\mathrm{H}}^{\mathrm{wake}}(r).\end{array}$$(7)

The total absorption column density, (N_H), becomes

$$\begin{array}{c}\hfill N_{\mathrm{H}}={\displaystyle {\int }_{r_{\mathrm{NS}}^{\prime }}^{r_{\mathrm{w}0}}}{n}_{\mathrm{H}}^{\mathrm{wake}}(r^{\prime })\mathrm{d}{r}^{\prime }+{\displaystyle {\int }_{r_{\mathrm{w}0}}^{\infty }{n}_{\mathrm{H}}^{\mathrm{SSW}}}(r^{\prime })\mathrm{d}{r}^{\prime },\end{array}$$(8)

where r_w0 is the extent of the wake along the line of sight.

5.2. Fitting of observed absorption column density

5.2.1. Stellar wind and wake structure for overall Vela X−1 observation

The four parameters that govern the structure of the wake are the opening angle of the wake (Θ), the phase lag parameter (ϵ), the neutron star immersion angle with respect to the leading edge of the wake (θ), and the velocity reduction factor (f). The variable parameters for the wind model are the ratio of the mass loss rate to the terminal velocity (Ṁ/v_∞), the smooth wind velocity exponent (β^SSW), and the wake velocity exponent (β^wake). The radius of the companion star is taken as 30 R_⊙ (van Kerkwijk et al. 1995a), while the conversion factor (χ_H) for the absorption column density was calculated as 0.7 (Wilms et al. 2000).

For a given set of four wake parameters, we setup the wake structure and calculate the absorption column density at a given orbital phase value using Eq. (8). Then we vary the three wind parameters to fit the observed absorption column density.

The wake parameters were varied about their best-fit values to estimate the errors from the chi-square distribution. The errors on the wind parameters were calculated using a Markov chain Monte Carlo sampler (Python: emcee; Foreman-Mackey et al. 2013; Goodman & Weare 2010) and their posterior distribution is shown in Fig. 11. The best-fit values for both wind and wake parameters, along with errors estimated at a 1σ confidence level, are given in Table 3. The corresponding best-fit absorption column density predicted by our model is shown in red in Fig. 12. The residuals indicate that the model parameters reproduce the data well. From the posterior distribution (Fig. 11), it is evident that there is a strong negative correlation between β^wake and Ṁ/v_∞. A similar but less strong correlation is observed between β^SSW and Ṁ/v_∞.

Fig. 11. Posterior distribution of the three wind parameters Ṁ/v∞, βsmooth, and βwake obtained from fitting the observed absorption column density to the model predicted absorption column density. The two extreme vertical lines are drawn at a 1σ confidence level. The contours correspond to two-dimensional 1σ, 2σ, and 3σ confidence levels.

Fig. 12. Top panel: Fitting of the orbital-phase resolved absorption column density (NH) to the modified smooth stellar wind model. Black: Observed absorption column density from MAXI/GSC with error bars at a 1σ confidence level. Red: Absorption column density with wake. Green: Absorption column density without wake. Bottom panel: Deviation of the data from the best-fitted model.

5.2.2. Intensity dependence of wake structure in Vela X−1

The wind parameters from the companion star are not expected to change over a timescale of decades; therefore, we assumed the same wind parameters for the different intensity levels. We only estimated the wake parameters by freezing the wind parameters to the one obtained from fitting of the observed orbital-phase resolved N_H (Table 3). The best-fit wake parameters were estimated based on the minimisation of chi-squares for the three intensity levels. The best wake parameters and their errors estimated at a 1σ confidence level are given in Table 3. The best-fit model and the observed absorption column density for the three intensity levels are shown in Fig. 13. To check the validity of our assumption of freezing the wind parameters to the overall best-fit value, we also varied the wind parameters for the three intensity levels. We found that the resulting best-fit parameters remained within the 1σ confidence interval of the value obtained from the overall Vela X−1 data.

Fig. 13. Top panel: Fitting of the absorption column density derived from the low-, mid-, and high-intensity orbital-phase-resolved spectroscopy to the modified smooth stellar wind model. Black: Observed absorption column density from MAXI/GSC with error bars at 1σ confidence level. Red: Absorption column density with wake. Green: Absorption column density without wake. Bottom panel: Deviation of the data from the best-fitted model.

6. Discussion

In this work, we have studied the sg-HMXB, Vela X−1, with almost thirteen years of observations by MAXI/GSC, which has helped us to better investigate the long-term orbital phase dependence of the spectral parameters. The stellar wind from OB-type stars can be characterised by the mass loss rate (Ṁ), the terminal velocity (v_∞), and the velocity exponent (β) (Eqs. 1, 2). In Vela X−1, the variation in the absorption column density with orbital phase shows a deviation from the symmetric behaviour expected from an isotropic stellar wind (top panels of Figs. 7, 9). The asymmetric variation in the absorption column density in Vela X−1 has been observed with different X-ray telescopes (Doroshenko et al. 2013; Kretschmar et al. 2021; Diez et al. 2022, 2023). Also, there are hydrodynamic simulations that have been performed to understand the origin of asymmetric variation in the N_H in Vela X−1 (Blondin et al. 1990; Manousakis et al. 2012). A wake-like structure was proposed by Doroshenko et al. (2013) to explain the asymmetric behaviour of the absorption column density. The onset of a wake-like structure has been reported by Diez et al. (2022, 2023) from pointed XMM-Newton and NuSTAR observations. With 13 years of MAXI/GSC data, we confirm the presence and the persistent nature of the wake-like structure (Fig. 7) in Vela X−1, which also exhibits an asymmetric variation in the absorption column density in our line of sight as was previously reported (see Sect. 1). We proposed a wake model in which the geometry of the wake is characterised by an opening angle, immersion angle, phase lag parameter, and a velocity reduction factor (Sect. 5.1). With a longer dataset, we were able to constrain the wind parameters and also quantitatively calculate the wake parameters, which was not attempted in Doroshenko et al. (2013). The best-fit wake and wind parameters obtained from our analysis are given in Table 3, where the errors are estimated at a 1σ confidence level. The wake parameters are consistent with Doroshenko et al. (2013). The wind parameters agree with previous reports (Nagase et al. 1986; Prinja et al. 1990; Watanabe et al. 2006; Falanga et al. 2015).

Another characteristic is the intensity fluctuation over the long term exhibited by Vela X−1 in the timescale of the orbital period (Fig. 1). In other pulsars such as LMC X-4 and Her X-1, the change in intensity on the timescale of a few orbital periods is periodic in nature, varying with a super-orbital period (Paul & Kitamoto 2002 for LMC X-4 and Katz 1973 for Her X-1) that is associated with the precession of binary plane and precession of the accretion disk, respectively. However, no such super-orbital periodicity in Vela X−1 has been reported so far.

The longer MAXI/GSC data allows us to study the variation in the wake structure with the intensity level. Orbital-phase resolved spectroscopy in the three different intensity states was performed. The main finding of the current work is that the accretion and/or photoionisation wake is present at all intensity levels (Fig. 13), which is understood from the variation in both the absorption column density with orbital phase (the top panel of Fig. 9) and the hardness ratios (Fig. 5). The increase in N_H (Fig. 9) and the hardness ratio (Fig. 5) happens earlier at a high intensity, indicating that the geometry of the wake is intensity-dependent. The varying wake parameters for the three intensity levels are given in Table 3 with the errors quoted at a 1σ confidence level. The consistent increase in the wake opening angle with the intensity level suggests that the wake structure prefers to be wider for higher intensity levels. Similar values of the velocity reduction factor suggest that the wake remains almost equally dense at all three intensity levels. We do not see such a trend for the phase lag parameter and the immersion angle.

Abalo et al. (2024), whose findings were published during the review process of this paper, have used the hardness ratio between 2.0−4.0 keV and 4.0−10.0 keV energy bands to study the X-ray absorbing medium along the line of sight at the individual orbit level. Similar to our analysis, they also show that large-scale structures in Vela X−1, such as the accretion wake trailing the neutron star, significantly influence the observed orbital variation in the absorption patterns, which is consistent with existing literature (Doroshenko et al. 2013; Malacaria et al. 2016). Their work shows that a considerable number of orbital hardness ratio profiles deviate from the average hardness ratio towards both hard and soft bands, and they attribute this to stochastic effects such as the clumpy nature of the stellar wind and the inherent variability of the source. Complementarily to Abalo et al. (2024), the current work focuses on the study of the variability of the accretion wake structure with changes in source luminosity.

Currently, because of the edge-on assumption (i = 90°), we have disregarded the third spatial axis. The immediate modification to the model that one could consider is it being in three dimensions. One could also consider a smooth transition of density at the wake boundaries instead of assuming an abrupt density change, which is beyond the scope of this work.

Acknowledgments

The authors would like to thank the comments and suggestions from the referee, which helped improve the quality of this manuscript. This research has made use of MAXI data provided by RIKEN, JAXA and the MAXI team.

References

All Tables

All Figures

	Fig. 1. ∼13-year long-term MAXI/GSC light curve of Vela X−1 with a bin size of one orbital period. The red, green, and blue arrows mark the low (< 0.36 counts s−1), mid (0.36−0.46 counts s−1), and high (> 0.46 counts s−1) intensity levels. The horizontal lines (cyan) separate the three intensity levels.
In the text

	Fig. 2. Energy-dependent orbital profiles 2.0−4.0 keV (purple), 12.0−20.0 keV (olive), and 2.0−20.0 keV (brown) are plotted with the bin size of 64 bins per interval. The hardness ratio (black) is plotted between 2.0−4.0 keV (soft) and 12.0−20.0 keV (hard) energy bands. The vertical lines (cyan) divide the orbital profiles into eclipse, egress, OOE, and ingress.
In the text

Fig. 3. Histogram of total intensity per orbit grouped into three intensity levels. The X axis represents the MAXI/GSC counts rate per orbit. The Y axis represents the number of orbits. The counts rate range for the low (red), mid (green), and high (blue) intensity levels are < 0.36 counts s−1, 0.36 − 0.46 counts s−1, and > 0.46 counts s−1, respectively. The total number of orbits in the three intensity levels are 222, 161, and 111, respectively.

In the text

	Fig. 4. 2.0−20.0 keV orbital profiles of Vela X−1 are plotted for low (red), mid (green), and high (blue) intensity levels with the bin size of 32 bins per interval. The vertical lines (cyan) divide the orbital profiles into eclipse, egress, OOE, and ingress.
In the text

	Fig. 5. Hardness ratio plots between the soft (S; 2.0−4.0 keV) and hard (H; 12.0−20.0 keV) energy bands for the low (red), mid (green), and high (blue) intensity levels. The inset in the upper right corner shows that the hardness ratio in the low-intensity level rises slowly compared to the high-intensity level in the phase range 0.4−0.6.
In the text

	Fig. 6. Time-averaged spectrum of Vela X−1 in 2.0−20.0 keV energy band obtained from MAXI/GSC. Top panel: Data (black crosses) and best-fit model (dashed line) components. Bottom panel: Deviation of the data from the best-fitted model spectrum.
In the text

	Fig. 7. Variation in spectral parameters (black) obtained from orbital-phase resolved spectroscopy. The error bars are at the 90% confidence level. The orbital profile (brown) is plotted in the 2.0−20.0 keV energy band. The vertical lines (cyan) are drawn to denote the four different regions (eclipse, egress, OOE, and ingress).
In the text

	Fig. 8. Three intensity-resolved spectra in the 2.0−20.0 keV energy band.
In the text

	Fig. 9. Variation in spectral parameters obtained from intensity-and-orbital-phase resolved spectroscopy for the low (red), mid (green), and high (blue) intensity levels. The error bars are at a 90% confidence level. The orbital profiles for low (red), mid (green), and high (blue) intensity levels are plotted in the 2.0−20.0 keV band.
In the text

	Fig. 10. Sketch of Vela X−1 with the photoionisation wake (in grey). NS (in black) is the neutron star at phase 0.4. The companion (in blue) is at the focus of the orbit. The angles θ and Θ are the immersion and the opening angles of the neutron star and the wake, respectively. The line AB is the line of apsides, and the line of sight is along the positive X axis.
In the text

	Fig. 11. Posterior distribution of the three wind parameters Ṁ/v∞, βsmooth, and βwake obtained from fitting the observed absorption column density to the model predicted absorption column density. The two extreme vertical lines are drawn at a 1σ confidence level. The contours correspond to two-dimensional 1σ, 2σ, and 3σ confidence levels.
In the text

	Fig. 12. Top panel: Fitting of the orbital-phase resolved absorption column density (NH) to the modified smooth stellar wind model. Black: Observed absorption column density from MAXI/GSC with error bars at a 1σ confidence level. Red: Absorption column density with wake. Green: Absorption column density without wake. Bottom panel: Deviation of the data from the best-fitted model.
In the text

Fig. 13. Top panel: Fitting of the absorption column density derived from the low-, mid-, and high-intensity orbital-phase-resolved spectroscopy to the modified smooth stellar wind model. Black: Observed absorption column density from MAXI/GSC with error bars at 1σ confidence level. Red: Absorption column density with wake. Green: Absorption column density without wake. Bottom panel: Deviation of the data from the best-fitted model.

In the text