© The Authors 2025

1. Introduction

Vela X-1 is one of the best-studied high-mass X-ray binary pulsars, since it was detected early (Chodil et al. 1967), shows persistent X-rays with clear pulsations (McClintock et al. 1976), and has a rich phenomenology in a wide range of wavelengths. The following parameters are taken from the recent review by Kretschmar et al. (2021). The system lies at a distance of ∼2 kpc with an intrinsic luminosity of ∼4 × 10³⁶ erg s⁻¹. It is formed by the B0.5 Ia supergiant HD77581 (M_⋆ ≈ 20 − 30 M_⊙, R_⋆ ≈ 30 R_⊙) and a neutron star (M_NS ≈ 1.7 − 2.1 M_⊙, R_NS ≈ 11 − 12 km). The neutron star orbits the supergiant in an almost circular orbit (e ≈ 0.0898) with a period of ∼8.964 days and a mean distance of ∼1.7R_⋆. X-ray eclipses of ∼20% of the orbit indicate that the orbital plane is observed near edge-on.

Being deeply embedded in the acceleration zone of the supergiant’s strong stellar wind (Martínez-Núñez et al. 2014; Kretschmar et al. 2021), the X-ray emission from the neutron star is significantly affected by photoelectric absorption at all orbital phases. In general, the observed absorption clearly decreases at eclipse egress and then rises again to higher values around orbital phase 0.5 (where mid-eclipse is phase zero), but with strong variations around this mean pattern (Sato et al. 1986; Doroshenko et al. 2013; Kretschmar et al. 2021; Abalo et al. 2024).

The X-ray spectrum of Vela X-1 is similar to that of other accreting X-ray pulsars. It has been mostly described by a powerlaw modified by variable strong photoelectric absorption at lower energies and an exponential cutoff beyond ∼20 keV, further modified by two cyclotron lines at ∼25 keV and ∼54 keV (Fürst et al. 2014). In addition, a strong iron fluorescence emission line is present.

Fürst et al. (2010) carried out a study of hard X-ray data from INTEGRAL. They found a log-normal distribution of observed fluxes, with a median absolute luminosity ⟨L_X⟩ = 5.1 × 10³⁶ erg s⁻¹ and multiplicative standard deviation σ ≈ 2.

Vela X-1 pulsates with a period of ∼283 s. The period varies on all timescales, from days to years, described by a random walk in pulse frequency (Boynton et al. 1984; Deeter et al. 1989). Within the history of pulse period determinations, the range of variability is about ±0.5 s (Kretschmar et al. 2021, Fig. 8).

The broadband pulse profile shows two broad pulses (see, e.g., Fig. 1 in Staubert et al. (1980) and Sect. 3.3), which have been related to the contributions from emission regions near the two magnetic poles (McClintock et al. 1976; Kraus et al. 1996), similar to many other accreting pulsars. The overall pulse profile shape has been found to be very stable over decades – compare, e.g., the profiles in McClintock et al. (1976), Staubert et al. (1980), Raubenheimer (1990), La Barbera et al. (2003), or Alonso-Hernández et al. (2022). The pulse profile is strongly energy dependent. Below 6 keV the shape is complex, with up to five visible peaks. At higher energies (≥10 keV) it evolves into a double-peaked profile. The relative strength of these two peaks also changes with energy (Raubenheimer 1990; La Barbera et al. 2003).

Variations of the X-ray emission in accreting X-ray pulsars from one pulse cycle to the next are commonly observed but relatively rarely studied more systematically. For Vela X-1, Staubert et al. (1980) compared individual pulse cycles and the mean pulse profile during a hard X-ray balloon observation, finding significant pulse-to-pulse variations, but without distinguishing between overall flux changes and shape variability. Using Suzaku/PIN data, Kretschmar et al. (2014) applied a rolling scale factor to the mean profile before the comparison to allow for flux changes. They found deviations between the scaled profile and the actual light curve on the order of a few 10%, but also regular structures in the residuals, and did not study changes in the shape of single pulse cycles as function of flux. Variations at short timescales have also been reported for other sources, e.g., 1A 0535+262 (Frontera et al. 1985; Klochkov et al. 2011), Cen X-3 (Müller et al. 2011), GX 301−2 (Fürst et al. 2011), or 4U 0114+65 (Sanjurjo-Ferrín et al. 2025).

Over the years, there have been multiple attempts to model the pulsed emission from Vela X-1, as well as that from other accreting X-ray pulsars (e.g., Sturner & Dermer 1994; Bulik et al. 1995; Leahy & Li 1995; Laycock et al. 2025). But these different studies have used different assumptions and approximations, and have arrived at very different solutions for the emission geometry. We, therefore, still lack a reliable base to explain the pulse profiles and their variations based on an emission model. Thus, we focus on an empirical study, which can be later used, for example, as benchmark to test emission models against.

In the following, we present a detailed analysis of pulse profile and pulse-to-pulse variations of Vela X-1 in the 1–10 keV range, based on three long observations of the source with XMM-Newton – see Sect. 2 below. Section 3 presents the energy-resolved light curves and pulse profiles obtained from these observations. In Sect. 4, we discuss the time- and flux-dependence of the pulse profiles. We analyze variations of individual pulse cycles in Sect. 5. The overall results are discussed in Sect. 6, and a summary and future outlook are given in Sect. 7.

2. Observations and data reduction

2.1. Observations

In this study, we analyzed the only three XMM-Newton observations of Vela X-1 with public data available in the archive as of late 2024, referred to hereafter as Obs 1, Obs 2, and Obs 3, respectively. In Table 1 and Fig. 1, we provide an overview of these observations and the orbital ranges they covered.

Fig. 1. Top view of the Vela X-1 system during the three XMM-Newton observations analyzed in this article. The red circles represent the periastron (filled) and the apastron (empty). The observer is located facing the system at x = 0 and at minus infinity along the y-axis. The gray zone represents the eclipse range.

For Obs 1, at a later phase of the orbit, no detailed study has been published so far. Obs 2, studying the egress from eclipse, has been studied by Martínez-Núñez et al. (2014) and Lomaeva et al. (2020). Obs 3 successfully aimed at observing the accretion wake moving into the line of sight (Diez et al. 2023).

Slightly different combinations of the X-ray instruments were used in each observation. In this publication we focus exclusively on data from the European Photon Imaging Camera pn-CCDs (EPIC-pn; Strüder et al. 2001) in timing mode, due to its high temporal resolution and ability to mostly avoid saturation and mitigate pile-up for this bright X-ray source.

2.2. Data reduction

To generate the event lists, we used the Science Analysis System (SAS) software v21.0 with the Current Calibration Files as of March 2024, starting from the Observation Data Files level, running epproc with default calibration. Although Diez et al. (2023) proposed a method in their Appendix to correct for an offset of instrumental and physical spectral lines toward high energies observed for Vela X-1 with EPIC-pn in timing mode, we extracted the events with the default Rate-Dependent PHA calibration, as this work focuses on timing analysis. We checked for flaring particle background in all observations, but no filtering was needed. In addition, we did not extract background for any of the observations, as the source was so bright that it illuminated the entire CCD plane. The event times were barycentered and deleted if found on bad pixels. We removed the outermost part of the point spread function wings to reduce the influence of background and possible dust scattering effects, leading us to extract events in the source region RAWX in [34:40] for Obs 1, RAWX in [31:43] for Obs 2, and RAWX in [32:44] for Obs 3.

The arrival times were corrected for the binary orbit. We specifically used the BinaryCor function from the Interactive Spectral Interpretation System (ISIS) v1.6.2-51 (Houck & Denicola 2000) scripts (ISISscripts), provided by ECAP/Remeis Observatory and MIT¹, and the ephemeris from Kreykenbohm et al. (2008).

All subsequent analysis was restricted to the 1–10 keV range, due to very low count rates below 1 keV and the relatively limited effective area of the EPIC-pn above 10 keV (see Figs. 11 and 12 in Sect. 3.2.2.1 of the XMM-Newton Users Handbook²). Since extraction regions differ between observations, and the instrument response evolves over time, the effective area varies across observations. Thus, we calculated the ratio between the average effective areas of the three observations (Table 2) and we applied them to the derived count rates (see Sect. 3.2 for more details).

2.3. Further software used

Besides the use of SAS to generate the event files, the data analysis and visualization of this paper were performed with Python v3.11.10 64-bit and, among others, the Astropy (Astropy Collaboration 2013) and Scipy (Virtanen et al. 2020) packages. We also used different functions from the spectral-timing software package Stingray³ (Huppenkothen et al. 2019; Bachetti et al. 2024) v2.1⁴, specifically pulse.search.epoch_folding_search to determine the pulse periods, Lightcurve.make_lightcurve for light curve generation, and pulse.pulsar.fold_events for pulse profiles.

3. Light curves and pulse profiles per observation

3.1. Pulse periods and pulse cycles

For consistency, we redetermined pulse periods for all three observations, using 256 as the number of bins and as oversampling factor in the search. For Obs 2 and Obs 3, we used the published periods from Martínez-Núñez et al. (2014) and Diez et al. (2023) as references, rounded to 283.39 s and 283.45 s, respectively. For Obs 1, we chose a reference value of 283.5 s, interpolating between earlier measurements with BATSE (Kretschmar et al. 2021, Fig. 8) and the periods found by Staubert et al. (2004) in INTEGRAL observations carried out in 2003. Uncertainties in the found periods were calculated using Eqs. 3 and B2 in Martin-Carrillo et al. (2012). The periods we found are shown in Table 3. For what concerns Obs 2 and Obs 3, our values agree within the uncertainties with the values determined in the earlier publications. The pulse periods remain stable within uncertainties during the individual observations, since maxima and marked minima of individual pulse cycles align in phase throughout the whole time intervals covered.

For the study of pulse profiles and their behavior during individual pulse cycles, we explored different values of 2ⁿ for the number of phase-bins and settled on 32 as a balance between resolution and typical significance of the count rate in each bin. In the following, we mainly discuss data at the resolution of full pulse cycles. For coherence between different data sets, we used the first occurrence of a marked minimum in the profiles, visible in individual pulse trains, as effective zero point for each observation. As reported in Table 3, these times are offset from the first arrival of each observation (51850.617, 53880.452, and 58606.927 MJD, for Obs 1, Obs 2, and Obs 3, respectively). Ignoring partial pulse cycles at the end of each observation, 189, 434, and 391 full pulse cycles are obtained for the three observations.

A further complication arose from bright flares in Obs 2 and Obs 3 saturating the onboard buffer and disturbing temporal information on timescales of individual pulse cycles (Martínez-Núñez et al. 2014; Diez et al. 2023). For this reason, we excluded the data between T = 76 ks and T = 86 ks in Obs 2, and T = 7.5 ks and T = 11 ks in Obs 3. The final sets of pulse cycles used are listed in Table 3.

3.2. Energy-resolved light curves and hardness ratios

Following Martínez-Núñez et al. (2014), we generated light curves in four energy bands: 1–3, 3–6, 6–8, and 8–10 keV. Poisson uncertainties were calculated using stats.poisson_conf_interval⁵, with “frequentist-confidence” as interval (for further details, see Maxwell 2011) and σ = 1. To correct for the differences in effective areas (Sect. 2.2), we normalized the derived count rates by the factors listed in Table 2.

In order to track changes in the source spectrum as function of time, we used the hardness ratio (HR), defined as

$$\begin{array}{c}\hfill \mathrm{HR}=\frac{H-S}{H+S},\end{array}$$(1)

where H and S denote the flux in the hard and soft band, respectively. In the following, the HR is always calculated between the 8–10 keV (hard) and the 3–6 keV (soft) bands, since this HR has been found to closely follow the evolution of the photoelectric absorption N_H obtained from spectral fitting for Obs 2 (Martínez-Núñez et al. 2014, Figs. 4 and 7) and Obs 3 (Diez et al. 2023, Figs. 5 and 9).

We present the energy-resolved light curves and the HR curves for the three observations in Fig. 2. To contextualize the HR evolution observed in our XMM-Newton data, Fig. 3 compares it to the long-term, orbit-averaged trend derived from MAXI/GSC observations (Abalo et al. 2024). A detailed discussion of this comparison is provided in Sect. 6.2.

Fig. 2. Energy-resolved light curves (top panels) and HRs (Eq. 1) between the 8–10 and 3–6 keV bands (bottom panels), with a time resolution of the respective pulse period, for the three observations. The width of each column is proportional to the duration of the corresponding observation. The light curves are plotted on a logarithmic scale for clarity. The colors represent different energy bands following the legend in the plot. Gray areas mark intervals of bright flares excluded for time-resolved analysis (see Sect. 3.1 and Table 3).

Fig. 3. HR evolution during the three XMM-Newton observations of Vela X-1 compared to the orbit-averaged MAXI/GSC trend (2009–2024), adapted from Abalo et al. (2024). MAXI HRs, based on Crab-like spectra, appear systematically harder and are rescaled for comparison with XMM-Newton, which uses count-based HRs.

Obs 1 shows relatively mild variability by a factor of about 6 in flux. The HR is around a mean value of about −0.5, and shows limited variability with a slight decay at the end of the observation.

Obs 2 is much more variable. The first data points are from a time when the neutron star was coming out of eclipse, a time period partially ignored in Martínez-Núñez et al. (2014) to avoid times close to the radiation belts, but included here since we concentrate on timing instead of spectral analysis. Even outside the excluded flaring period, observed fluxes vary by a factor of more than 360 in the 1–3 keV energy band, or ∼140 in the 8–10 keV range. The HR shows a strong evolution, with both the highest and the lowest values across all data sets reached in this observation. The impact of the correlated changes in absorption are clearly visible in the flux evolution of the light curve at 1–3 keV.

Obs 3 presents an intermediate picture between the two other observations. Flux varies by up to a factor of ∼27 in the 1–3 keV range and ∼6 in the 8–10 keV band. The HR rises markedly in the second half of the observation, which has been associated to the signature of the accretion wake moving across the line of sight (Diez et al. 2023).

3.3. Energy-resolved pulse profiles and hardness ratios

We present in Fig. 4 the energy-resolved pulse profiles averaged over the time intervals taken into account for each observation (Table 3), along with their corresponding HRs. Since the profiles were obtained from the events and not from the light curves, we applied again the factors from Table 2. A comparison between the normalized pulse profiles of the three observations is displayed in the last column, where the normalization was performed by dividing the individual profiles by their mean rate. We applied this method to all the normalized pulse profiles shown in this article.

For the further discussion, we defined two parts of the profile, A and B, associated with the two peaks visible at higher energies, which are generally linked to emission from the two polar regions. Based on the 32 bins used throughout this study, the boundary between these parts is right after the 14^th phase-bin. Motivated by the sub-peaks, visible especially below 3 keV, we further split these intervals, such that A is formed by the two sections A1 (bins 0–7) and A2 (8–13), while B consists of sections B1 (14–19), B2 (20–24), and B3 (25–31). These boundaries are indicated in Fig. 4. This segmentation is similar to the one used by Raubenheimer (1990), but with flipped definition of A and B, and their region C roughly corresponding to our B1.

Fig. 4. Mean pulse profiles of the three observations in the four energy bands used throughout this article and the corresponding HRs between the 8–10 and the 3–6 keV bands. In the rightmost column, we compare the profiles across observations, normalizing the individual profiles by their respective mean count rates. The error bars are smaller than the line width in nearly all of the panels. Dashed lines indicate the five defined phase regions, A1 to B3, as explained in the text.

While generally similar, the profiles shown in Fig. 4 display differences at the level of these components. A2 almost always contains the maximum of the profiles, except for the 1–3 keV profile of Obs 3. In general, the ratio between A1 and A2 varies between different observations. B1, B2, and B3 also show variations in their relative strength, both as a function of energy and between observations.

The last column of Fig. 4 confirms the general stability of the pulse profile shape across many years and different levels of source brightness. The main differences are at low energies, where pulse amplitudes are reduced, especially for Obs 1.

The bottom panels of Fig. 4 show the HRs between the 8–10 and 3–6 keV profiles. For all observations, the HR minimum coincides in pulsed phase with the boundary between A and B, while the maxima coincide with the boundaries between the subpeaks.

4. Pulse profile variations with time and flux

The mean pulse profiles in Fig. 4 average data over extended periods with substantial flux and hardness variability (Fig. 2). If the profile shape depends on the flux, mean shapes will be dominated by the brightest pulse cycles. We next examine how pulse shape evolves with time and varies with flux.

4.1. Variations with time within individual observations

We studied the evolution of the energy-resolved pulse profiles within each observation. For this analysis, we aimed to split the total number of pulse cycles in each observation (Table 3) as evenly as possible, while finding a compromise between the goals of higher time resolution to follow variations in the light curves, and statistical quality, which requires integration over a sufficient number of pulse periods. We settled on a baseline of 13 pulse periods per profile, corresponding to ∼1 h of integration time, with some variations as explained below. This choice was also driven by the will to have profiles of equal integration time before and after the bright flares excluded in the analysis (Sect. 3.1). All of the above led to the following splits per observation. For Obs 1, there are 14 full-length profiles and a final one covering 7 pulse periods. Regarding Obs 2, up to the flare, we have one profile covering 8 pulse periods and 20 full-length profiles; after the flare, 10 full-length profiles. For Obs 3, there are 2 full-length profiles before the flare; after this, 26 full-length profiles and one covering 14 pulse periods. In Fig. 5, we can visualize the changes in profile shapes during the observations, with example profiles at different times marked in the light curves shown above.

Fig. 5. Top panels: 8–10 keV light curves with a time resolution of an individual phase-bin (∼8.86 s). We mark in color the time ranges used to build the time-resolved profiles shown in the lower panels.

The most pronounced pulse profile changes occur in the 1–3 keV band, where pulsations are nearly absent for most of Obs 1 – except near the end when HR drops (Fig. 2) – and similarly in the initial and final time ranges of Obs 2 and Obs 3, respectively, when HR is also high. Profiles at higher energies remain more stable, but still show variations in the relative brightness between parts (A and B) and between sections (A1, A2, B1, B2, and B3). For Obs 1, the amplitude of A slightly decreases with respect to that of B as time goes on. This is even more pronounced in Obs 2, where the amplitude of B surpasses that of A at the end of the observation. For Obs 3, A1 and B3 are relatively brighter in the earlier profiles than in the later ones.

4.2. Variations with flux within and across observations

To investigate a possible relation between the source flux and the observed pulse profile shape, we classified all pulse cycles by their mean count rate in the 8–10 keV energy range. We selected this energy band as it is the least affected by absorption. We chose count rate intervals covering the observed range variability and ensuring that each observation was covered by at least three intervals. This led to the following count rate boundaries between intervals: 4, 8, 15, 20, and 28 counts/s (Fig. 6).

Fig. 6. Distributions of the 8–10 keV mean count rate per pulse cycle. The black dashed lines indicate the boundaries between flux intervals. The pulse profiles obtained for each flux interval and observation are shown in Fig. 7.

The flux-resolved pulse profiles for different observations (Fig. 7) are visibly more similar for higher fluxes. The similarity across observations also appears to increase with energy. Note that the pulse profiles averaged over the brightest pulse cycles – a range only covered within Obs 2 – show a significantly different shape, with A and B being almost equally bright.

Fig. 7. Normalized flux-resolved pulse profiles, obtained as explained in Sect. 4.2. Each column corresponds to a different flux range, given by the count rate C (counts/s) in the 8–10 keV band averaged over the pulse cycle. The number followed by “pc” in the top row is the number of pulse cycles per observation and flux level. Obs 2 covers the whole flux range, while Obs 1 and Obs 3 only have data within 3 and 4 flux intervals, respectively. The distributions of the 8–10 keV mean count rate per pulse cycle and the flux limits are shown in Fig. 6.

5. Variations of individual pulse cycles

At least some of the variations between pulse profiles described in the previous sections will be driven by the known variations between pulse cycles, which are also visible in the light curves of all three observations. In the following, we quantify these variations.

5.1. Total flux variations between consecutive pulse cycles

As a measure for intrinsic variations in X-ray brightness from one pulse cycle to the next, we calculated the ratio of total flux in the 8–10 keV band between consecutive pulse cycles. For Obs 1, the minimum and maximum ratios are 0.61 and 1.75, respectively, while 90% of the values lie in the range 0.78–1.26. For Obs 2, the corresponding ratios are 0.58 and 2.25, and 0.75–1.37, while for Obs 3 they are 0.60 and 1.58, and 0.79–1.28.

5.2. Variable contributions from the individual pulse sections

The fraction of phase covered by A and B is 0.4375 and 0.5625, respectively. However, the fraction of flux is almost perfectly identical on average. Considering all individual pulse cycles of all observations and energy bands, the mean flux fraction in A is 0.498, while the median is 0.499, and for 90% of the pulse cycles this fraction is in the range 0.408–0.587, i.e., a scatter of about 10% from the even distribution.

We also calculated the flux fraction per pulsed phase section (A1, A2, B1, B2, and B3) in each pulse cycle. The mean and standard deviation of these fractions confirm that the contribution of each section to the total emission is fairly constant across all observations and energy bands, with variations of the order of 10–20% of the average. A1 and B1 tend to diminish with energy, while A2 and B2 tend to rise. B3 does not seem to evolve with energy range.

Overall variations may not be representative of fluctuations between consecutive pulse cycles. Therefore, we quantified the variability of each pulsed phase section (X) from one pulse cycle to the next with

$$\begin{array}{c}\hfill {\mathrm{R}}_{X,(j+1)/j}=\frac{{(F_X/F_T)}_{j+1}}{{(F_X/F_T)}_j},\end{array}$$(2)

where j is the pulse count, F_X is the flux in the selected section, and F_T is the total flux in the selected pulse cycle. Considering all energy bands, A2 fluctuates less from one pulse cycle to the next than the rest of the pulsed phase sections, with the geometric standard deviation⁶ of R_A2 between 1.14 and 1.25. On the opposite extreme we find B2 and B3, with the corresponding geometric standard deviations in the ranges 1.17–1.26 and 1.19–1.32, respectively.

5.3. Variations in shape compared to the mean profile

5.3.1. Synthetic light curves scaling the mean profile

Kretschmar et al. (2014) compared the extracted light curves with “ideal” light curves, which were obtained in two steps. They first repeated the mean pulse profile as many times as pulse cycles were in the extracted light curve. The second step was to scale the predicted curve to the extracted one, so that they were comparable regardless of the flux variability.

Following a similar approach, we created “synthetic light curves” (Slc) by repeating the mean pulse profile (PP) for each observation sufficient times to cover the pulse cycles used (Table 3), normalizing this curve, and multiplying by the mean count rate of the real light curve (Rlc) for each pulse cycle

$$\begin{array}{c}\hfill {\mathrm{Slc}}_{i,j}=\frac{{\mathrm{PP}}_i}{\overline{{\mathrm{PP}}_i}}\times {\left(\overline{{\mathrm{Rlc}}_i}\right)}_j,\end{array}$$(3)

where i denotes the phase-bin.

We also tried other scaling approaches, including separate scale factors for the regions A and B, but we finally chose this one as it produced the fewest artifacts when visually inspected. An example comparing a few pulse cycles of real data with a simple repeated profile and the scaled version is shown in Fig. 8.

Fig. 8. Demonstration of the approach to generate the Slc by scaling the repeated PP to match the mean flux of the corresponding pulse cycle of the Rlc. The top panel shows the Rlc and the repeated PP before scaling, while the bottom panel shows the Rlc and the final Slc. See Sect. 5.3.1 for more details.

5.3.2. Quantitative comparison parameters

To measure the degree of similarity between the Rlc and the Slc, we calculated the Pearson Correlation Coefficient (PCC) with the Scipy function stats.pearsonr⁷ (setting 0.9 as the confidence level). Using their definition of PCC, our coefficient per pulse cycle is

$$\begin{array}{c}\hfill {\mathrm{PCC}}_j=\frac{{\sum }_{i=1}^N({\text{Rlc}}_{i,j}-{(\overline{{\mathrm{Rlc}}_i})}_j)\times ({\text{Slc}}_{i,j}-{(\overline{{\mathrm{Slc}}_i})}_j)}{\sqrt{{\displaystyle \sum _{i=1}^N}{({\text{Rlc}}_{i,j}-{(\overline{{\mathrm{Rlc}}_i})}_j)}^2\times {\displaystyle \sum _{i=1}^N}{({\text{Slc}}_{i,j}-{(\overline{{\mathrm{Slc}}_i})}_j)}^2}},\end{array}$$(4)

where N is the number of bins per pulse cycle. PCC can take any value between −1 (perfect negative correlation) and 1 (perfect positive correlation). PCC = 0 means that there is no correlation at all.

To quantify the resemblance between the curves on timescales of the duration of an individual bin, we measured the absolute percentage error (APE)

$$\begin{array}{c}\hfill {\text{APE}}_{i,j}=\frac{|{\text{Rlc}}_{i,j}-{\text{Slc}}_{i,j}|}{{\text{Rlc}}_{i,j}}.\end{array}$$(5)

5.3.3. Distribution of PCC values per pulse cycle

The PCC values for individual pulse cycles range from somewhat anticorrelated (−0.47) to almost perfectly matching (0.97), showing a wide variety of shapes for individual pulse cycles. Examples of good and bad matches between the curves are shown in Fig. 9. To assess the relation between these variations and flux, we plotted the PCC against the mean count rate (Fig. 10).

Fig. 9. Comparison between the Rlc and the Slc for some pulse cycles of Obs 2 in the 3–6 keV range. The pulse cycle number is indicated right after the symbol # in the top panels. The generation of the Slc is explained in Sect. 5.3.1. The similarity between the Rlc and the Slc decreases from left to right, as quantified by PCC (obtained with Eq. 4) and APE (obtained with Eq. 5).

Fig. 10. PCC (obtained with Eq. 4) between the Rlc and the Slc per pulse cycle, plotted against the mean count rate (in logarithmic scale) of the corresponding pulse cycle. The colors represent different energy bands following the legend in the plot. Error bars are not shown for a clearer visualization of the scatter of individual points. The green and red lines mark the interval between perfect correlation and no correlation.

The Rlc and the Slc behave, in general, similarly, regardless of the observation. Nonetheless, Fig. 10 shows that the differences between the curves increase as flux decreases. Due to absorption, the lowest fluxes and highest differences are detected in the 1–3 keV range.

To further quantify the difference in shape between the Rlc and Slc pulse cycles, we calculated the 0.1 and 0.9 quantiles (Q_0.1 and Q_0.9, respectively) of the PCC distributions (Table 4). In this analysis, we sought to exclude pulse cycles without evident pulsations. At the beginning of Obs 2, with the neutron star coming out of the eclipse, the signal of individual pulse cycles is not significantly different from a constant close to zero, even at the highest energies. Making a cut at a mean count rate of 2 counts/s in the 8–10 keV band, we excluded pulse cycles 0–26, 41, and 44 of Obs 2 (eclipse egress started between pulse cycles 21–49), while the other observations were not affected.

The lowest Q_0.1 values for all observations are in the 1–3 keV range. The highest Q_0.1 values in Obs 2 and Obs 3 correspond to the 8–10 keV range, while the highest Q_0.9 values are distributed between the 6–8 keV (Obs 1), 8–10 keV (Obs 2), and 3–6 keV (Obs 3) bands. In summary, the Rlc and the Slc are least similar in the 1–3 keV range, and most similar in the 8–10 keV band. Still, the 3–6 and 6–8 keV values are very close to those of the highest energy band. Comparing the observations, Obs 3 has overall the quantile values closest to 1, while Obs 2 is the furthest away.

5.3.4. PCC-resolved pulse profiles

In order to test the connection between PCC values and systematic shifts in the emission pattern, we chose to create “PCC-resolved pulse profiles”. To do so, we divided the range in PCC values found in the 8–10 keV band (labeled PCC_8 − 10 in the following) across all three observations into five intervals, with PCC_8 − 10 values of 0.6, 0.8, 0.85, and 0.9 as limits (Fig. 11). This choice ensured to have at least some pulse cycles of each observation within the chosen intervals. We then created mean profiles from the individual pulse cycles falling within a certain PCC_8 − 10 range, as shown in Fig. 12.

Fig. 11. Distributions of the 8–10 keV PCC (Eq. 4) between the Rlc and the Slc. The black dashed lines indicate the limits between PCC ranges. The PCC-resolved light curves and pulse profiles are shown in Fig. 12.

Fig. 12. Top panels: 8–10 keV PCC-resolved (see Sect. 5.3.4) light curves with a time resolution of an individual phase-bin (∼8.86 s). The numbers followed by “pc” are the number of pulse cycles per PCC8 − 10 range. First two rows of the bottom panels: normalized PCC-resolved pulse profiles in the 3–6 and 8–10 keV bands. Last row of the bottom panel: HR between the 8–10 and 3–6 keV normalized pulse profiles.

The light curves in Fig. 12 show no evident temporal pattern for PCC_8 − 10 in any observation. The highest PCC_8 − 10 values usually correspond to pulse cycles of high flux; except for pulse cycles ∼230–267 of Obs 2, where during a very active period before the omitted massive flare, individual pulse cycles also deviate strongly from the mean shape. Regarding the PCC-resolved profiles, the amplitude of A1 and A2 increase with respect to that of B2 and B3 (in both energy bands) as PCC_8 − 10 increases, while the amplitude of B1 is much more stable. Note that the HR shape barely varies with PCC_8 − 10. It is also worth mentioning that, in the 8–10 keV range, the PCC-resolved profile from the lowest PCC_8 − 10 range of Obs 1 is significantly different to that of Obs 2 and Obs 3. We visually confirmed that this is due to the small number of pulse cycles (5) combined for Obs 1 in this range.

6. Discussion

6.1. Overall pulse profile stability

The analysis of three observations of Vela X-1 with the same X-ray telescope, XMM-Newton, widely spaced in time (2000, 2006, and 2019) and at different orbital phases and luminosities, confirms again the long-term stability of the general pulse profile shape, and thus emission properties, averaged over many pulse cycles. Smaller, but significant differences, remain in the strength of the visually defined subcomponents of the profiles (Fig. 4 and Sect. 5.2).

6.2. Impact of absorption

Figure 3 places the HR evolution from our three XMM-Newton observations within the broader orbital context of Vela X-1, as illustrated by the long-term averaged MAXI curve from Abalo et al. (2024). The MAXI trend shows systematic modulation across the orbit, with pronounced hardening near eclipse ingress and egress due to increased absorption in dense structures. The XMM-Newton tracks align with this behavior: Obs 1 and the end of Obs 3 coincide with phases of high obscuration, likely caused by wake structures, while the start of Obs 2 shows elevated HRs as the neutron star emerges from eclipse behind the extended stellar atmosphere.

Significant variations are observed in the “soft” pulse profiles within individual observations (Fig. 5), driven by changes in HR and thus N_H (Figs. 2 and 3). Increased absorption and scattering at lower energies lead to large count rate fluctuations and, at times, to a complete loss of the typical pulse shape. This behavior is reflected in the high scatter of individual pulse cycles around the mean profile, particularly in the 1–3 keV band, as quantified by the low PCC values (Fig. 10, Table 4).

The almost flat 1–3 keV pulse profiles in Fig. 5 are also a consequence of high absorption episodes. In Obs 1, the 1–3 keV flux remains persistently low, consistent with dense wake structures crossing the line of sight at this orbital phase (Kaper et al. 1994; Doroshenko et al. 2013; Abalo et al. 2024). As a result, pulsations are poorly defined in this energy range almost throughout the entire observation. A similar suppression is seen at the end of Obs 3, where comparable structures again obscure the neutron star (Malacaria et al. 2016; Diez et al. 2023). At the beginning of Obs 2, with the neutron star emerging from eclipse, photons were absorbed or scattered by the extended atmosphere of the supergiant (Martínez-Núñez et al. 2014). In all three cases, the duration of the suppressed pulsations suggests that large-scale, orbitally modulated structures dominate the absorption, rather than short-lived wind clumps.

6.3. Flux dependence

Flux variations can yield changes in the height of the accretion column, leading to changes in the fraction of photons affected by light bending and, ultimately, to variations in the shape of the pulse profile (Falkner 2018). The flux resolved analysis (Fig. 7) indicates that the overall shape of the pulse profile is similar at all flux levels – excluding the brightest profiles only observed in Obs 2, which are discussed below. Except for these, the profiles appear also more similar across different observations at higher source flux.

In contrast, various accreting X-ray pulsars do show significant variations in their pulse shapes as a function of luminosity – for recent examples see, e.g., Chhotaray et al. (2024), Sharma et al. (2024), Thalhammer et al. (2024), Du et al. (2025). In bright states, such changes are usually associated with the idea of passing a “critical luminosity” (L_crit), above which the emitted radiation is capable of braking the flow above the surface of the neutron star, creating a radiative shock and an extended column (Basko & Sunyaev 1976; Becker et al. 2012; Mushtukov et al. 2015). At first order, L_crit depends on the mass accretion rate and the magnetic field strength of the source. Following Becker et al. (2012), Diez et al. (2022) estimated L_crit for different assumptions on accretion geometry and magnetic field strength, finding ranges from ∼0.13–0.15 × 10³⁷ erg s⁻¹ to ∼5–6 × 10³⁷ erg s⁻¹. Using a competing approach, Mushtukov et al. (2015) stated that L_crit was likely to be around 10³⁷ erg s⁻¹. Diez et al. (2022) obtained luminosities between ∼0.2–1 × 10³⁷ erg s⁻¹ for one NuSTAR observation strictly simultaneous with Obs 3, and argued that Vela X-1 lies between the sub- and super-critical accretion regimes, but that all such estimates have to be treated with caution. Based on the fluxes derived by Martínez-Núñez et al. (2014) for the same data as our brightest data points in Obs 2, we estimated an intrinsic luminosity of about 10³⁷ erg s⁻¹ for these. While the varying absorption does not allow to simply extrapolate this relation to the whole data set, all of the above indicates that the source was close to L_crit during the brightest pulse cycles analyzed in this work.

The brightest flux-resolved pulse profiles are different in that the amplitudes of A and B are almost identical (Fig. 7). Thirty of the thirty-one pulse cycles that make up these profiles correspond to an almost continuous and relatively long time range (pulse cycles ∼218–267) preceding the flare of Obs 2, which was ignored in this work due to the disturbance of temporal information (see Sect. 3.1 for more details). Martínez-Núñez et al. (2014) found that the wind density started to increase during this observation at orbital phase ∼0.185. Considering that they used T₉₀ as time of phase zero, and that we used T_ecl, their mentioned orbital phase corresponds to ∼0.160 (pulse cycle 186) in our data. They claimed that an increase in the emission of the neutron star followed this wind enhancement, probably due to the corresponding gradual accretion of more mass. This also points to the source approaching L_crit during this period, which is reflected in the pulse profile as the change in relative amplitude between A and B and, ultimately, in the low PCC_8 − 10 values obtained for these pulse cycles (Fig. 12).

6.4. Nonperiodic variations on short timescales

In the 8–10 keV range, the majority of pulse cycles are, at most, ∼30% brighter or dimmer than their predecessor, while the most dramatic change was quantified by a factor of 2.25 (Obs 2). All of the above is in agreement with Staubert et al. (1980), who found smooth pulse-to-pulse intensity variation by up to a factor of 2.

Concerning the subcomponents, B1 and B2 are less prone to vary overall, while A1 shows the largest variations. In terms of consecutive pulse cycles, A2 seems to vary less than the rest, and B2 and B3 are more prone to fluctuate. It is interesting to note that the standard deviation of the flux fraction in B2 is one of the lowest and, at the same time, this section fluctuates slightly more than the rest from one pulse cycle to the next.

As found through the PCC-resolved analysis (Sect. 5.3.4, Fig. 12), the minority of pulse cycles with relatively low PCC tend to have B (especially B2 and B3) relatively more prominent, and A somewhat suppressed, compared to mean profiles. At the same time, the HR shape does not change with PCC. If we assume that A and B are related to the two magnetic pole regions, variations in relative brightness imply variations in the relative accretion to the two poles. These abnormalities could be, e.g., sudden and slight increases (or decreases) in the density of the accreted matter.

Overall, pulse-to-pulse variations in Vela X-1 are very stochastic, apparently following no temporal patterns nor trends. This is again in line with Staubert et al. (1980), but also with Frontera et al. (1985), who found “chaotic spiky” pulse shape variations on short timescales in the accreting X-ray pulsar 1A 0535+262.

7. Summary and outlook

We analyzed three long XMM-Newton observations, from different years and orbital phases, to study pulse profiles and pulse-to-pulse variations in the accreting X-ray pulsar Vela X-1. We confirm that the overall shape of the pulse profile is stable, although we detected non-negligible differences on all timescales. The magnetic field and X-ray emission geometry appear very stable, and Vela X-1 seems to be accreting continuously but somewhat erratically from the stellar wind, leading to flux and pulse shape variations. The general structure of Vela X-1’s pulse profile does not heavily depend on the flux. We also observed that, as the brightness increases, the similarity between profiles of different observations also rises. Below ∼3 keV, pulsations can be strongly suppressed by the variable strong absorption, common in Vela X-1 at certain orbital phases. Variations at higher energies seem to be caused by changes in the density of the accretion flow. On short timescales, the pulse shape is also more stable when the source is brighter, and variations are nonperiodic. There is no compelling evidence to suggest that any pulsed phase section varies significantly more than others.

It would be interesting to compare these profiles with new ones, obtained from any future observations of Vela X-1 at different orbital phases. This would allow us to further study pulse profile variations along the binary orbit. More detailed models of the pulsed emission of Vela X-1 at soft and hard X-rays, and accounting for flux variations, could allow for a direct comparison with the observed profiles, and provide additional insights into the causes of the observed variations. Extending this analysis to other accreting X-ray pulsars with complex pulse profiles would help us to better understand the variations presented here.

Acknowledgments

We thank the referee for very positive and constructive comments that identified some remaining gaps in the presented work, and led us to several additions to the text of this publication. VMF acknowledges ESA/ESAC and Universidad Complutense de Madrid (UCM) for the agreement that allows students from UCM to carry out their master’s thesis at ESAC, as this work started within the framework of this agreement. IEM acknowledges support from the grant ANID FONDECYT 11240206. SMN acknowledges support from grant PID2021-122955OB-C41 funded by MICIU/AEI/10.13039/501100011033 and by “ERDF/EU”. This research was supported by the International Space Science Institute (ISSI) in Bern, through ISSI International Team project #495 (Feeding the spinning top). Work on this project profited from activities funded by the Space Science Faculty of the European Space Agency.

References

All Tables

All Figures

	Fig. 1. Top view of the Vela X-1 system during the three XMM-Newton observations analyzed in this article. The red circles represent the periastron (filled) and the apastron (empty). The observer is located facing the system at x = 0 and at minus infinity along the y-axis. The gray zone represents the eclipse range.
In the text

Fig. 2. Energy-resolved light curves (top panels) and HRs (Eq. 1) between the 8–10 and 3–6 keV bands (bottom panels), with a time resolution of the respective pulse period, for the three observations. The width of each column is proportional to the duration of the corresponding observation. The light curves are plotted on a logarithmic scale for clarity. The colors represent different energy bands following the legend in the plot. Gray areas mark intervals of bright flares excluded for time-resolved analysis (see Sect. 3.1 and Table 3).

In the text

	Fig. 3. HR evolution during the three XMM-Newton observations of Vela X-1 compared to the orbit-averaged MAXI/GSC trend (2009–2024), adapted from Abalo et al. (2024). MAXI HRs, based on Crab-like spectra, appear systematically harder and are rescaled for comparison with XMM-Newton, which uses count-based HRs.
In the text

Fig. 4. Mean pulse profiles of the three observations in the four energy bands used throughout this article and the corresponding HRs between the 8–10 and the 3–6 keV bands. In the rightmost column, we compare the profiles across observations, normalizing the individual profiles by their respective mean count rates. The error bars are smaller than the line width in nearly all of the panels. Dashed lines indicate the five defined phase regions, A1 to B3, as explained in the text.

In the text

	Fig. 5. Top panels: 8–10 keV light curves with a time resolution of an individual phase-bin (∼8.86 s). We mark in color the time ranges used to build the time-resolved profiles shown in the lower panels.
In the text

	Fig. 6. Distributions of the 8–10 keV mean count rate per pulse cycle. The black dashed lines indicate the boundaries between flux intervals. The pulse profiles obtained for each flux interval and observation are shown in Fig. 7.
In the text

Fig. 7. Normalized flux-resolved pulse profiles, obtained as explained in Sect. 4.2. Each column corresponds to a different flux range, given by the count rate C (counts/s) in the 8–10 keV band averaged over the pulse cycle. The number followed by “pc” in the top row is the number of pulse cycles per observation and flux level. Obs 2 covers the whole flux range, while Obs 1 and Obs 3 only have data within 3 and 4 flux intervals, respectively. The distributions of the 8–10 keV mean count rate per pulse cycle and the flux limits are shown in Fig. 6.

In the text

	Fig. 8. Demonstration of the approach to generate the Slc by scaling the repeated PP to match the mean flux of the corresponding pulse cycle of the Rlc. The top panel shows the Rlc and the repeated PP before scaling, while the bottom panel shows the Rlc and the final Slc. See Sect. 5.3.1 for more details.
In the text

	Fig. 9. Comparison between the Rlc and the Slc for some pulse cycles of Obs 2 in the 3–6 keV range. The pulse cycle number is indicated right after the symbol # in the top panels. The generation of the Slc is explained in Sect. 5.3.1. The similarity between the Rlc and the Slc decreases from left to right, as quantified by PCC (obtained with Eq. 4) and APE (obtained with Eq. 5).
In the text

Fig. 10. PCC (obtained with Eq. 4) between the Rlc and the Slc per pulse cycle, plotted against the mean count rate (in logarithmic scale) of the corresponding pulse cycle. The colors represent different energy bands following the legend in the plot. Error bars are not shown for a clearer visualization of the scatter of individual points. The green and red lines mark the interval between perfect correlation and no correlation.

In the text

	Fig. 11. Distributions of the 8–10 keV PCC (Eq. 4) between the Rlc and the Slc. The black dashed lines indicate the limits between PCC ranges. The PCC-resolved light curves and pulse profiles are shown in Fig. 12.
In the text

Fig. 12. Top panels: 8–10 keV PCC-resolved (see Sect. 5.3.4) light curves with a time resolution of an individual phase-bin (∼8.86 s). The numbers followed by “pc” are the number of pulse cycles per PCC8 − 10 range. First two rows of the bottom panels: normalized PCC-resolved pulse profiles in the 3–6 and 8–10 keV bands. Last row of the bottom panel: HR between the 8–10 and 3–6 keV normalized pulse profiles.

In the text