© The Authors 2025

1. Introduction

High-mass X-ray binaries (HMXRBs) are systems formed by a compact object (i.e., a neutron star or a black hole) orbiting a massive star (the companion). These are excellent laboratories where the process of accretion and the structure of the companion’s stellar wind can be studied (Martínez-Núñez et al. 2017). The study of HMXRBs is important in order to understand the evolution of close binaries, which may become compact-object mergers and eventually sources of gravitational waves and/or short γ-ray bursts. They also provide insight into the behavior of matter at extreme gravitational and magnetic fields. Understanding these processes is fundamental to modern astrophysics and has been the driver of multiple theoretical and observational studies.

The HMXB 4U 0114+65 was discovered by Dower et al. (1977) in the SAS 3 Galactic survey. The companion star is V^* V662 Cas, a luminous BI1a supergiant (Reig et al. 1996). With an orbital period of ∼11.6 d, the neutron star (NS) orbits the companion deeply embedded into its wind at an orbital radius of ≈1.34 − 1.65 R_⋆, thereby offering the opportunity to probe the inner wind of the B1 supergiant star. It shows remarkable temporal and spectral variability across various timescales (Crampton et al. 1985; Hu et al. 2017). In Table 1 a brief summary of the system’s main characteristics is provided.

4U 0114+65 is a non-eclipsing system (Pradhan et al. 2015) and hosts an X-ray pulsar with a ∼2.6 h spin period, which makes it one of the slowest X-ray pulsars ever known. In order to explain such a short spin period, Li & van den Heuvel (1999) and Sanjurjo-Ferrín et al. (2017) suggested that this object could have been born as a magnetar. Long periods of X-ray pulsars might point to high magnetic fields of accreting NSs (Ikhsanov & Beskrovnaya 2013). Other slowly rotating NSs were thought to be magnetars, similar to the case of the slowly rotating NSs found in the supernova remnant RCW103 (Rea et al. 2016; D’Aì et al. 2016).

Magnetars are NSs that possess very strong magnetic fields, B ≃ 10¹⁴ − 10¹⁵ G (Duncan & Thompson 1992) and they can be present in HMXRBs (Bozzo et al. 2008). Alternatively, Koenigsberger et al. (2006) proposed that the observed pulse in this system could be explained by the accretion of a structured wind produced by tidally driven oscillations in the B-supergiant star photosphere induced by the closely orbiting NS.

The pulse period of 4U 0114+65 has been changing quickly over the years. It was initially estimated to be around 2.78 hours (Finley et al. 1992). Subsequent observations reported shorter periods: in 2000, Hall et al. (2000) found a period of approximately 2.73 hours, followed by (Bonning & Falanga 2005) who measured it to be about 2.67 hours. In 2006, Farrell et al. (2006) obtained a period of roughly 2.65 hours. Wang (2011) noted a further decrease in the spin period, from approximately 2.67 hours to 2.63 hours between 2003 and 2008, with a corresponding increase in the NS spin rate of approximately 1.09 × 10⁻⁶ s⁻¹; in Sanjurjo-Ferrín et al. (2017) a period of 2.6 h was reported. The pulse period is long enough to allow the detailed study of intrapulse morphology.

Furthermore, the system’s long-term behavior exhibits super-orbital modulations and torque reversals, suggesting complex interactions between the NS and the accreting material (Farrell et al. 2006; Sood et al. 2006). Hu et al. (2017) proposed that the changing behavior of the NS spin period could be caused by the effect of a transient disk.

In Sanjurjo-Ferrín et al. (2017) the spin behavior (duration and fast evolution) of 4U 0114+65 was explained within the framework of the quasi-spherical settling accretion theory proposed by Shakura et al. (2012). In this mode of wind accretion, typically present in sources exhibiting moderate X-ray luminosity (below ≃4 × 10³⁶ ergs⁻¹), a convective quasi-spherical shell forms above the NS magnetosphere. The inflow rate of plasma from this shell through the magnetosphere is governed by the cooling mechanisms acting on the hot plasma, primarily driven by Compton and radiative energy losses.

Relative recent observations with NuSTAR and XMM-Newton have shed new light on the system’s activity, revealing features such as a corotating interaction region (CIR), whose presence is commonly invoked to explain the periodic variability observed in the UV spectral lines that form in stellar winds of B-type supergiants (Lobel & Blomme 2008) as well as the modulation of X-ray emission from O stars (Oskinova et al. 2001; Massa et al. 2014), off-states, and potential signatures of cyclotron resonant scattering features (CRSFs) (Abdallah et al. 2023), which could be at odds with the magnetar nature of the system (see Fig. 1).

Fig. 1. Pole-on sketch of the system showing the relative sizes of the companion star and the NS orbit (NS size not to scale, the observer would be situated in a perpendicular plane with respect to the orbit). The light curves are displayed parallel to the orbit. Green represents the 9 uninterrupted Chandra light curves from the dataset analyzed in this paper scaled to the count rate. Red represents a former XMM-Newton light curve taken in 2015 scaled accordingly to fit the figure.

In this paper we present a study of 4U 0114+65 based on 200 ks of Chandra grating observations. This data allowed us to infer the characteristics of the accretion flow along several phases of the spin and orbital period.

2. Observations and analysis

The dataset analyzed consists of nine different Chandra observations spread along the NS orbit, taken at different epochs. The Chandra telescope, launched by NASA in 1999, employs a system of four nested grazing-incidence mirrors coated with iridium. These mirrors focus X-rays onto detectors with an angular resolution of 0.5 arcseconds. Chandra operates over an energy range of 0.1 − 10 keV using two primary instruments: the High Resolution Camera (HRC) for imaging and the Advanced CCD Imaging Spectrometer (ACIS) for detailed spectroscopy. Positioned in a highly elliptical orbit, it minimizes Earth’s interference and enables long-duration observations (Weisskopf et al. 2000). For this analysis the ACIS data were utilized.

The log of the observations, the corresponding orbit, and the order within each orbit are collected in Table 2.

The spectra and response (arf and rmf) files were generated using standard procedures with CIAO software v4.15, CalDB 4.15, (Fruscione et al. 2006). The first dispersion orders (m = ±1) of HEG and MEG were extracted and combined as described in Sect. 4. The data were considered useful when in the 0.2–10 keV range.

The spectra were analyzed and modeled with the ISIS¹ package. The emission lines were identified using the ATOMDB² data base. The orbital modulations present within this dataset were analyzed using the XRAYBINARYORBIT Python package (Sanjurjo-Ferrín et al. 2024)³.

3. Timing analysis

Within this dataset there are nine different uninterrupted light curves spread across four different orbits. Each light curve is referred to as the orbit in which it occurred followed by a number in chronological ascending order. As an example, the first light curve is Orb 1-1 (see Table 2 and Fig. A.1). The orbital phases were calculated following the ephemeris provided by Hu et al. (2017) (see Eq. (2) in their paper). The light curves were extracted in 10 s bins and then further re-binned to reach a signal-to-noise ratio of 7 in each bin. Thus, the final light curves have different time resolutions in the range 30–200 s. The pulse period of the NS can be easily discerned by the eye. To characterize the spin period is not straightforward however (see below).

To facilitate the analysis of the light curve and perform a spectral analysis, we segmented each observation into peaks and valleys. This segmentation was achieved by re-binning the light curve using a weighted mean approach. A valley was identified at each minimum, while a transition to a peak was considered when the count recovery exceeded three times the median difference in counts between bins for each light curve (see Fig. 2 for an illustrative example).

Fig. 2. Example of peak and valley division on Orb 1-1. Upper panel: High-energy (3–10 keV, red) and low-energy (0.2–3 keV, blue) light curves. The vertical lines indicate the divisions between the peaks and valleys. Lower panel: Color ratio (CR), defined as the high-energy light curve divided by the low-energy light curve.

In Table 3 we provide comprehensive data on the light curves, including their duration, phase, NS period if detected, pulse fraction, duration of valleys and peaks, and weighted mean color ratio (CR).

3.1. Spin period

The slow pulsation of the X-ray source can be appreciated at first sight. However the pulse is not consistently detected in every observation, as evidenced by the Lomb-Scargle periodogram (Scargle 1982). Specifically, during Orbit 14, the pulse is absent in the first two observations, although this might be caused by their duration being close to the NS spin period. In observations 5 and 6 the pulse is weakly detected. Furthermore, in Obs 32-1 the pulse appears irregular and interrupted, which leads to its nondetection in the periodogram.

The pulse fraction is a measure of the amplitude of variability in a periodic or pulsating signal. It can be calculated as

$$\text{PF}=\frac{I_{\text{max}}-I_{\text{min}}}{I_{\text{max}}+I_{\text{min}}},$$

where I_max and I_min are the maximum and minimum intensity observed during the pulse, respectively.

Within this source, the pulse fraction varies considerably across observations (Table 3), ranging from 0.6 to 0.2. There is a strong positive correlation (r = 0.9) between the evolution of the pulse fraction through the orbit and a sinusoidal orbital evolution.

The source exhibits secular spin-up with periods of reversal, as documented in Hu et al. (2017) and Abdallah et al. (2023). To accurately calculate the spin period, we utilized a Lomb-Scargle periodogram on the observations that exhibit a pulse fraction higher than 0.3, where the pulse could be detected: Orbit 1-1, Orbit 14-4, and Orbit 31-1. The weighted mean period, along with its error, was determined to be 9050 ± 100 seconds. The source is thus continuing the spin-up trend previously reported.

3.2. NS pulse shape

A remarkable characteristic of this source is that the pulses are different from each other. Moreover, Orb 1-1, Orb 14-3, and Orb 32-1 show a decreasing evolution of the peak count rate, both in the high- and low-energy light curves.

To further investigate this issue we selected all the peaks within this dataset lasting more than 4000 s and calculated the weighted mean of the count rates for the high-energy light curve (whc, 3–10 keV), the low-energy light curve (wlc, 0.3–3 keV), and the hardness ratio (HR), defined as $\frac{\mathrm{whc}-\mathrm{wlc}}{\mathrm{whc}+\mathrm{wlc}}$, and used a k-means algorithm to classify them. k-means (MacQueen 1967) is a popular clustering algorithm used in machine learning to partition a dataset into a predetermined number of clusters. It works by iteratively assigning data points to the nearest cluster centroid and then updating the centroids based on the mean of the points assigned to each cluster.

The best result obtained was for four groups with a silhouette score of 0.52. The silhouette score (Rousseeuw 1987) measures how similar an object is to its own cluster compared to other clusters, with values ranging from −1 to 1. A score of −1 indicates no separation between clusters, while a score of 1 suggests a perfect separation. A score of 0.52 indicates a reasonably good clustering result, suggesting that the data points are appropriately grouped into distinct clusters. Based on their characteristics, we labeled each of the groups as Bright, Intermediate, Faint, or Absorbed (Fig. A.1 and Table 4).

In the Orb 32-1 light curve (Fig. 3), the first peak lasts 13.7 ks, approximately ×1.5 the duration of the rest of the other peaks (∼9.0 ± 0.7 ks on average). The subsequent peaks recover the ∼9 ks duration. However, a dramatic change occurs. The first valley is placed in the NS corresponding valley phase (0.9 − 0.1), but after the first 13 ks peak the NS phase is shifted by 0.5, so that now the peaks are where the valleys should be.

Fig. 3. Orb 32-1 light curve. The green vertical lines divide the peaks from the valleys. The NS phase 0.9 − 0.1 is shadowed in blue.

In Sanjurjo-Ferrín et al. (2017) a low-luminosity episode, making one of the NS spin valleys longer than expected, was interpreted as a substantial departure from the spherical accretion symmetry. The Bondi radius is an idealization, which could change by a factor of 2 or 3 in different directions because the stellar wind velocity changes on the characteristic scale of 10¹¹ cm, which is comparable to R_B (Table 1). The stellar wind could then be momentarily captured in a nonspherical way. Since the magnetospheric radius and the Bondi radius differ only by a factor of a few, the captured matter can reach the magnetosphere nonspherically, thereby enabling a switching between the magnetic poles of the NS. Our hypothesis is that a pole switching could cause a longer than normal peak and an immediate recovery to the usual values, maintaining, from now on, a 0.5 NS phase shift.

To analyze the pulse profile, we folded the peaks identified within the k-means classification for each of the four clusters (see Fig. 4). Prior to folding, the peaks were normalized to ensure that each peak contributed equally to the average pulse profile. The folded pulse for the Bright group displays a broad maximum, but is skewed to the left, with higher intensity observed at lower NS spin phases. The folded pulse for the Intermediate group indicates a subtle double-peaked structure, with the left peak appearing more pronounced than the right. In both cases, the maximum is situated around NS phase 0.4. In contrast, the Faint and Absorbed groups exhibit irregular profiles with large error bars, precluding detailed analysis. Bulik et al. (2003) noted the scarcity of double-peaked accretion-powered pulsars with respect to what is expected from theoretical models, attributing this scarcity to the likely alignment of the magnetic and rotational axes, which would increase the probability of single-peaked pulse profiles. When compared to the pulse profiles calculated using XMM-Newton data by Sanjurjo-Ferrín et al. (2017) (see their Fig. 4), the pulses in this dataset share characteristics with those observed in the first part of their observation (T1).

Fig. 4. Normalized NS spin pulse shapes. The averaged pulses of the four distinct orbits included in this dataset are shown in the top panels, while the averaged pulses for the k-means classifications are displayed in the bottom panels. To compute the average pulse shape, all pulses were normalized beforehand to ensure that each contributed equally to the final result.

3.3. Spikes

The very slow cadence of the NS pulsations allows us to appreciate a great deal of interesting details within each peak that reveals the complexity of the accretion processes. In this section we want to focus on the flickering variability exhibited by this source, especially within the peaks where some spikes can be observed.

There are two main challenges in the spike detection. On the one hand, we have to distinguish them from random noise. On the other hand, we need to retain a sufficient time resolution to obtain reliable results. Unfortunately, within this dataset the low-energy light curve is very suppressed, so we focused solely on the spikes detected in the high-energy light curve.

To distinguish real spikes from the noise, we conducted a comparative analysis using a semi-supervised dip-detection algorithm based on anomaly detection with autoencoders. Autoencoders are known for their effectiveness in anomaly detection tasks as they excel at learning and reconstructing normal data patterns (Bank et al. 2020). By leveraging the reconstruction error, they offer a robust mechanism to identify data points that deviate from the norm.

For our study we generated synthetic light curves by de-trending the original data and creating a random time series with the same statistical properties as the original light curve. The spikes detected in these synthetic light curves were then used as a training set for our autoencoder model. Spikes observed in the real light curve that deviated significantly from those detected in the synthetic data were identified and preserved as genuine spikes likely caused by astrophysical phenomena. Additionally, we calculated the signal-to-noise ratio (S/N) based on the mean uncertainty of the signal around each spike, and spikes with an S/N greater than five were also accepted.

We define the prominence of each spike as their increase in count rate with respect to their baseline divided by the baseline. Thus, the prominence represents the increase in percentage of the count rate (see Fig. 5).

Fig. 5. Detected spikes at their respective orbital phases. The radial axis represents the prominence. The size of the bubbles is proportional to their duration. The color gradient is used to visually support the prominence gradient.

After obtaining the real spikes, we applied a Gaussian mixture clustering algorithm. Unlike the k-means algorithm, which assumes clusters of spherical shape, Gaussian mixture provides a more flexible clustering approach by modeling the data as a mixture of multiple Gaussian distributions, allowing for elliptical cluster shapes and different cluster sizes (Reynolds 2009). This approach led to better results in this particular case, achieving a silhouette score greater than 0.6, which indicates the presence of two well-separated clusters. A higher silhouette score (between 0.7–0.8) was achieved by clustering each observation’s spikes individually, pointing to a difference in the spike durations and prominences through the different observations (see Fig. 6).

Fig. 6. Prominence vs. duration of spikes detected in this dataset. The blue squares and red triangles represent the two separated clusters obtained by a Gaussian mixture algorithm. The leftmost panel provides the results obtained for all the observations combined, while the three panels on the right display the results for the three individual observations that produced the better results in number of spikes detected, and thus in clustering.

Two distinct groups of spikes are observed within this dataset, clearly distinguished by their duration and prominence. We refer to them as short and long spikes based on their duration. The short spikes, within each individual observation, exhibit a narrow range of durations (see the three right panels in Fig. 6). Notably, the duration of these short spikes is inversely proportional to the brightness of the source for each observation, as measured by the weighted mean number of counts. In contrast, no specific trend is observed for the long-duration spikes (see Fig. 7). Similar flaring activity has been attributed to the development of magnetospheric instability in NSs in other HMXRBs, such as supergiant fast X-ray transients (Sidoli et al. 2019).

Fig. 7. Average spike duration for each of the two clusters (with the upper panel displaying the short spikes and the lower panel displaying the long spikes) as predicted by the Gaussian mixture algorithm for each observation. The x-axis represents the brightness of each observation, measured as the error-weighted count average (whc).

4. Spectral analysis

We used two different criteria to perform the spectral analysis. On the one hand, we analyzed the combined spectra for each of the four orbits in this dataset. To achieve this we isolated and combined the peaks within each observation.

These spectra are Orb 1, Orb 14, Orb 31, and Orb 32. However, the limited count rate prevented us from conducting the same exercise with each orbit for the valleys. Thus, all the valleys were combined into one single spectrum, which we refer to as valley combined. This classification is referred to as orbital spectra. The best fit for each spectrum is represented in Fig. B.1, and the best fit parameters are collected in Table 5 (see below for the description of the model used).

On the other hand, we combined the spectra based on the peak morphology grouped as suggested by the k-means algorithm (see Sect. 3.2 and Table 4). We refer to them as Bright, Intermediate, Faint, and Absorbed. This classification is referred to as k-means spectra. The best fit for each spectrum is represented in Fig. 8 and the best fit parameters are collected in Table 6.

Fig. 8. Upper row (left to right): Spectra, model, and χ2 values for Bright, Intermediate, and Faint. Lower row (left to right): Spectra, model, and χ2 values for Absorbed, and all spectra displayed for comparison.

The continuum has been described using the bulk motion Comptonization (bmc) model. bmc is an analytical model that describes the Comptonization of soft seed photons by matter undergoing relativistic bulk motion (Titarchuk et al. 1997). The model parameters are the characteristic blackbody temperature of the soft photon source, a spectral energy index (α), and an illumination parameter characterizing the fractional illumination of the bulk motion flow by the thermal photon source.

The absorption is modeled by the Tuebingen-Boulder interstellar medium (ISM) absorption model (Tbnew). This model calculates the cross section for X-ray absorption by the ISM as the sum of the cross sections due to the gas-phase ISM, the grain-phase ISM, and the molecules in the ISM (Wilms et al. 2000). The detected emission lines were modeled as Gaussians. Their parameters are collected in Table C.1.

With the aim of maintaining consistency and comparability, in our first attempt we employed the same model described in Sanjurjo-Ferrín et al. (2017). However, for the present dataset the partial covering was not needed, and a soft excess could be observed at low energies. It was successfully fitted with an absorbed bbody. The temperature and absorption of this soft excess were 0.04 − 0.07 keV and a 0.5 − 0.8 × 10²² cm⁻², respectively, in every case. Thus, we decided to fix the temperature to the more frequent value (0.06 keV) and fix the absorption to the interstellar value toward the source (0.8 × 10²² cm⁻²).

The complete model is described by the equation

$$\begin{array}{cc}\hfill F(E)& =exp(-N_H\sigma (E))\times (\mathrm{BMC}(E)+{\displaystyle \sum _{i=1}^6}(G_i))\hfill \\ \hfill & +exp(-N_H^{\mathrm{ISM}}\sigma (E))\times \mathrm{BBody},\hfill \end{array}$$(1)

where G_i represent the Gaussians describing the emission lines. The orbital and k-mean averages were both fit with the same model.

The radius of the source emitting the seed soft photons, which are subsequently Comptonized, can be estimated assuming that the source is radiating as a blackbody of area $\pi R_W^2$ (Torrejón et al. 2004)

$$\begin{array}{c}\hfill R_W=0.6\sqrt{L_{34}}{(kT)}^{-2}[\mathrm{km}],\end{array}$$(2)

where L₃₄ is the luminosity of the bmc component in units of 10³⁴ ergs s⁻¹ and kT is the temperature of the seed soft photons. In Tables 5 and 6 we show the best fit parameters for each spectrum.

The Orb 32 spectra and the Absorbed spectra both exhibit a high absorption column. These spectra where taken in orbital phases close to ϕ ∼ 0. In addition to the six lines detected in all the other spectra, they exhibit numerous emission lines that are not observed elsewhere in this dataset.

Particularly intriguing is the presence of a bump near ∼2 keV of these two same spectra. To fit our spectra, depending on the number of counts and the quality of the data, we applied a minimum S/N ratio of 5 − 15 per bin. Upon closer examination, reducing the data binning to a minimum signal-to-noise ratio of 2.5 per bin revealed that the bump is formed by several emission lines. This change in the signal-to-noise ratio for the binning in the 0.6 − 2.2 keV energy range corresponds to a significant increase in resolution, from R = 4 − 9 to R = 25 − 34 in terms of energy. This enhancement is well within the HEG and MEG capabilities of the Chandra ACIS instrument (see Weisskopf et al. (2000) for detailed information). The result is shown in Fig. 9. To fit this section of the spectra, we thawed the bbody parameters and add Gaussians to fit the detected emission lines. Interestingly, the bbody parameters were compatible with those obtained in the rest of the spectra. Thus, we finally set the temperature and absorption column to 0.06 keV and 0.8 × 10²² cm⁻², respectively, and let the Gaussians vary freely. These line parameters are collected in Table D.1. The emission lines detected on the high-energy continuum are reported in Table E.1 (see Eqs. (3) and (1) for a comparison):

$$\begin{array}{cc}\hfill F(E)& =exp(-N_H\sigma (E))\times (\mathrm{BMC}(E)+{\displaystyle \sum _{i=1}^{6+4}}(G_i))\hfill \\ \hfill & +exp(-N_H^{\mathrm{ISM}}\sigma (E))\times (\mathrm{BBody}\times +{\displaystyle \sum _{i=1}^8}(G_i)).\hfill \end{array}$$(3)

Fig. 9. Emission lines observed in the low-energy range within the spectra Orb 32 and the absorbed spectra.

5. Discussion

5.1. The circumstellar absorption

The absorption column traversed by the radiation emitted by the compact object on its way to the observer depends on the wind density it crosses in its path toward the observer. This path depends on the orbital phase, the orbital separation, eccentricity, inclination, and argument of periapsis. To estimate the density we used the wind structure predicted by the CAK models, adopting an acceleration parameter beta β = 0.8 (Castor et al. 1975).

Since the system has an inclination and eccentricity different from 0, the circumstellar absorption column varies through the orbit, being much higher in phases close to ϕ = 0. This explains the high absorption observed in Orb 32 or in the Absorbed spectra as they were taken during and around ϕ = 0 (Tables 5 and 6).

By comparing the evolution of N_H through the orbit with the expected theoretical values we can constrain the inclination and mass-loss rate of the system. In doing so, we assume that the stellar wind is un-ionized and spherically distributed. Even when this is an oversimplification of the real situation, this analysis can provide insight into the system configuration.

Utilizing the N_H collected in Table 5 for the orbital spectra Orb 1, Orb 31, and Orb 32 and the orbital parameters summarized in Table 1, we obtain an estimation on the mass-loss rate and inclination of the system. Orb 14 comprises several different orbital phases, and thus this data is not suitable for this calculation; the same problem holds for the K-means spectra.

The complexity of the model prevents us from utilizing a classical least-squares approach; hence, the data were fitted using a particle swarm optimization (PSO) algorithm. In PSO, each particle represents a potential solution to an optimization problem. The particles move through the solution space to find the optimal parameter configuration (Bonyadi & Michalewicz 2017).

The best fit (r = 0.9) was achieved for an inclination of 40 ± 3° and a lower limit for the mass-loss rate of 8.6 ± 1.7 × 10⁻⁷ M_⊙ yr⁻¹; however, it must be taken into account that the companion’s wind could be partly ionized. This result is in agreement with a long-term study of the optical spectra (Reig et al. 2016) where a strong Hα emission was found, sometimes showing a P Cygni profile, and therefore pointing to a substantial wind density (Ṁ ∼ 10⁻⁶ − 10⁻⁷ M_⊙ yr⁻¹).

Fig. 10. Representation of the system. The dotted circle represents the NS orbit and the white circle the companion star. The red arrows represent the path that the emitting radiation transverses in its path to the observer. This path is different in length and wind density depending on the orbital phase in which the radiation was emitted, provided that the system has an inclination and eccentricity different from zero. The mass-loss rate and inclination could be inferred by the differences found within the NH in three different orbital phases, corresponding to Orb 1-1, Orb 31-1, and Orb 32-1 (see Tables 5 for the NH values and 3 for the orbital phases).

A visualization of the system is shown in Fig. 10. Due to the system’s nonzero eccentricity, the density of the stellar wind encountered by the NS varies through the orbit. This variation arises because the distance between the compact object and the companion star changes as the object moves along its elliptical trajectory. Interestingly, the observed pulse fraction varies in a manner very similar to the density through the orbit expected for the best fit parameters (see Fig. 11), providing further validation of the model.

Fig. 11. Density variation encountered by the NS as it travels around the eccentric orbit (red continuous line). The measured pulse fraction (Table 3) might be correlated with the density of the NS environment.

5.2. Soft excess

An intriguing characteristic of many HMXB systems is that their spectra show an excess of emission below 0.5 keV, known as the soft excess. Its origin is still not clear and might differ between sources, having several different contributions (Hess et al. 1993; Hickox et al. 2004).

At first glance, the Absorbed and Orb 32 spectra exhibited an excess around 2 keV, initially interpreted as either a blackbody component or a reflection of the continuum. However, upon closer inspection, testing different data binning, a set of emission lines became discernible. One possibility then is that part of the soft excess is produced by unresolved emission lines generated in the stellar wind. This bump is overwhelmed by the NS emission (bmc component) in less absorbed spectra.

In our case, we successfully fitted this excess emission with an absorbed low-temperature bbody (kT = 0.06 keV and $N_{\mathrm{H}}^{22}=0.8$). A straightforward candidate for this emission is the companion star. However, the luminosity expected for a B companion star should be ∼×10³² erg s⁻¹, while the unabsorbed luminosity for this component turns out to be ∼×10³ higher (see Tables 6 and 5).

Throughout the orbital spectra, the radius of this blackbody and its normalization are very similar, but in the K-means classification a progression can be observed, being enhanced during the Bright spectra and more suppressed in the Faint spectra. This is consistent with the scenario in which (at least part of) the soft excess are unresolved emission lines in the stellar wind, enhanced by a higher X-ray illumination from the compact object.

5.3. Fe Kα line

The Fe Kα emission line in 4U 0114+65 is much weaker than normally found in HMXBs. The companion in 4U 0114+65 (B1Ia) is the coolest among the supergiant X-ray binaries (Krtička et al. 2015). With a T_eff = 24 kK, 4U 0114+65 is in the middle of the bistability jump (see lower panel in Fig. 3 presented by Vink et al. 1999). This jump is caused by the sudden change in the ionization balance of Fe (III, IV), which is a major contributor to the acceleration of the wind. Sanjurjo-Ferrín et al. (2017) suggest that this could have a strong impact on the efficiency of the mechanism to form and/or destroy clumps, and conclude that 4U 0114+65 have the typical thick wind of a supergiant star, but with a much lower degree of clumping.

Fig. 12. Ionization map in the plane of the orbit: Shown (from left to right) are the Bright, Intermediate, Faint, and Absorbed spectra. The delimited regions represent the Fe Kα region, i.e., with an ionization parameter lower than 100 and higher than 25. We note that the orbital phase is ignored as the spectra are not phase-resolved. For visualization purposes, the orbital phase ϕ = 0.12 is used.

In order to have a strong Fe Kα line, the ionization parameter $\xi =L_{\mathrm{X}}/n(r_{\mathrm{X}})r_{\mathrm{X}}^2$, where n(r_X) is the local particle density and r_X is the distance from the X-ray source (the NS), must be ≤10². If we assume the former calculated mass-loss rate of 8.6 ± 1.7 × 10⁻⁷ M_⊙ yr⁻¹, the distance to the NS at which ξ equals to 100 is 0.4, 0.3, 0.1, and 0.2 R_⋆ for the Bright, Intermediate, Faint, and Absorbed spectra, respectively. This indicates that, within the Bright spectra, the X-ray source is powerful enough to illuminate the gas at a distance of at least 0.4 R_⋆ from the compact object.

An ionization map is shown in Fig. 12, in which the Fe Kα formation area (i.e., ξ between 25 and 100) is represented by the delimited region inside the map. It is large in the Bright spectra and very compact in the Faint spectra. The area size displays a positive trend with the Fe Kα intensity (Fig. 13). We note that, in this particular case, the orbital phase is ignored because the spectra are classified using k-means clustering and are not phase-resolved. For visualization purposes, the orbital phase ϕ = 0.12 is used.

Fig. 13. Fe Kα formation region (ionization parameter in the range 50 − 100) vs. the Fe Kα line area.

5.4. NS pulse to pulse variability: A change of accretion regime

The observed pulse-to-pulse luminosity variations and the spikes present within the light curve of 4U 0114+65 can be explained within the framework of the theory of quasi-spherical settling accretion onto slowly rotating magnetized NS (see Shakura et al. 2012) and its further development in subsequent works (Shakura et al. 2014). A comprehensive overview of this theory and its applications can be found in the review articles by Shakura et al. (2015) and Shakura & Postnov (2017).

In sources with moderate X-ray luminosity (below ∼4 × 10³⁴ erg s⁻¹), a hot convective quasi-spherical shell forms above the NS magnetosphere. The matter enters the magnetosphere due to Rayleigh-Taylor instability (RTI) at the magnetospheric boundary. This accretion processes in HMXRBs has been known and studied for a long time (see Arons & Lea 1976 and Elsner & Lamb 1977). The settling velocity in this regime is approximately two times smaller than the free-fall velocity, and the RTI is mediated by Compton and radiative plasma cooling.

We hypothesize that a long spike duration t_B, reflects the variability in the density and velocity of the gravitationally captured stellar wind and should scale with the free-fall time from the Bondi radius:

$$\begin{array}{c}\hfill t_B=\frac{R_B^{3/2}}{{(\mathrm{GM})}^{1/2}}\sim 2\times {10}^3\mathrm{s}.\end{array}$$(4)

The spread in this timescale can be attributed to variations in the stellar wind velocity, which influences the Bondi radius and, consequently, the duration of t_B. For instance, a 10% variation in the wind velocity results in approximately a 30% variation in t_B as the Bondi radius is inversely proportional to the square of the wind velocity.

A shorter timescale of variability, where t_A represent short spike durations, is produced by the RTI near the magnetosphere. The growth rate of the RTI scales as

$$\begin{array}{c}\hfill \sigma =\frac{1}{t_A}=\sqrt{\mathit{Agk}},\end{array}$$(5)

where A is the Atwood number (assumed to be close to 1 in our case), $g=\frac{\mathrm{GM}}{R_A^2}$ is the gravitational acceleration at the magnetospheric boundary, and k = 2π/λ is the RTI wave number. The most unstable modes scale with R_A, giving

$$\begin{array}{c}\hfill t_A\sim \frac{R_A^{3/2}}{{(\mathrm{GM})}^{1/2}}·\end{array}$$(6)

Plugging in our system’s Alfvén radius, a 1.4 M_⊙ NS mass, and for A ∼ 1, we get t_A ∼ 150 s, in agreement with the observed spike duration (see Fig. 6). Since the magnetospheric radius is inversely proportional to the luminosity of the source, a negative correlation is expected between t_A and the brightness, and this what is observed (Fig. 7, upper panel), supporting RTI in the magnetosphere as the most likely origin of these short-duration spikes.

The light curve exhibits a cyclic pattern where bright pulses are followed by intermediate and faint ones; the bright luminosity is up to 14 times higher than the faint luminosity. Eventually, the intensity of the pulse recovers and the cycle begins again (see Fig. 14).

Fig. 14. Light curves showing the decreasing pulse intensity trend (Orb 1-1, Orb 14-3, and Orb 31-1).

In this scenario, matter accumulates within the magnetospheric boundary until the pressure triggers more efficient Compton cooling and a faster RTI, leading to increased accreted matter flow, and thus the source brightness. As the settling velocity is approximately two times smaller than the free-fall velocity, the accumulated matter is gradually depleted, reducing the accretion rate, and thus the luminosity until sufficient matter accumulates again, triggering another phase of the efficient Compton cooling. Considering the Atwood number $A\sim 1=\frac{{\rho }_1-{\rho }_2}{{\rho }_1+{\rho }_2}$, where ρ₂ is the density inside the magnetosphere and ρ₁ represents the density outside the magnetosphere, and assuming ρ₂ remains approximately constant with ρ₂ ≪ ρ₁, an increase in ρ₁ (corresponding to greater matter accumulation above the Alfvén radius) would result in a shorter t_A, and vice versa. This behavior aligns with our observations (shorter spikes in brighter spectra, corresponding with high accumulation of matter and longer spikes in faint spectra, when the matter above the magnetosphere is depleted).

The spectral parameters derived form the spectral model support this scenario. The local absorption column N_H is about two times larger in the Bright spectra than in the Faint spectra, suggesting an accumulation of circumstellar matter (see Table 6). The α parameter from the bmc component changes significantly, from ∼0.1 in the case of the Bright spectra, indicating a very efficient Comptonization mechanism, compared to ∼4 for the Faint spectra.

6. Summary and conclusions

In this analysis of the 4U 0114+65 system, based on 200 ks of Chandra observations, we examine the variability in the accretion process of one of the slowest pulsars known. The key results are the following:

• Circumstellar absorption. The system’s inclination angle with respect to the observer and the companion’s mass-loss rate were successfully constrained by fitting the variation of N_H through the orbit using a PSO algorithm. The result suggests an inclination of 40 ± 3° and a lower limit for the mass-loss rate of 8.6 ± 1.6 × 10⁻⁷ M_⊙ yr⁻¹, in alignment with previous studies.

• Soft excess. The soft excess present in the X-ray spectrum, especially around 2 keV, is at least partly formed by unresolved emission lines. This feature could only be revealed in the more absorbed spectra, when the emission produced by the NS at low-energy ranges is partly masked by the circumstellar absorption.

• Fe Kα line emission site. This emission line, generally strong within HMXRBs, is unusually weak in this source. We hypothesize that 4U 0114+65 has the typical thick wind of a supergiant star, but with a much lower degree of clumping. In order to have a strong Fe Kα line, the ionization parameter ξ must be ≤10². This ionization parameter is reached at distances from the NS of 0.4, 0.3, 0.1, and 0.2 R_⋆ for the Bright, Intermediate, Faint, and Absorbed spectra, respectively. This indicates that, within the Bright spectra, the X-ray source is powerful enough to illuminate the gas at distances of at least 0.3 R_⋆ resulting in a higher Fe Kα line flux.

• Accretion regime cycle. The pulse-to-pulse variability and spikes in the light curve can be explained within the theory of quasi-spherical settling accretion. Longer-duration spikes correspond to variations close to the Bondi radius, while shorter-duration spikes are likely caused by RTI near the magnetosphere. As the settling velocity is approximately two times lower than the free-fall velocity, matter accumulates within the magnetospheric boundary until the pressure triggers efficient Compton cooling, leading to increased source brightness. The accumulated matter is gradually depleted, reducing the accretion rate, and thus the luminosity, until sufficient matter accumulates again, triggering another phase of efficient Compton cooling.

Data availability

The data analyzed in this study can be found in the Chandra archive under the observation identification numbers 23432, 24479, 24480, 24481, 24482, 24483, 16177, 26178, and 16180. The Appendices can be found on Zenodo: https://zenodo.org/records/14645709.

Acknowledgments

This research has been funded by the ASFAE/2022/02 project from the Generalitat Valenciana. N. Schulz and M. Nowak were supported by NASA Chandra grants GO3-24018A and GO3-24018B, respectively. We acknowledge the constructive criticism of the referee whose comments improved the content of the paper.

References

All Tables

All Figures

Fig. 1. Pole-on sketch of the system showing the relative sizes of the companion star and the NS orbit (NS size not to scale, the observer would be situated in a perpendicular plane with respect to the orbit). The light curves are displayed parallel to the orbit. Green represents the 9 uninterrupted Chandra light curves from the dataset analyzed in this paper scaled to the count rate. Red represents a former XMM-Newton light curve taken in 2015 scaled accordingly to fit the figure.

In the text

	Fig. 2. Example of peak and valley division on Orb 1-1. Upper panel: High-energy (3–10 keV, red) and low-energy (0.2–3 keV, blue) light curves. The vertical lines indicate the divisions between the peaks and valleys. Lower panel: Color ratio (CR), defined as the high-energy light curve divided by the low-energy light curve.
In the text

	Fig. 3. Orb 32-1 light curve. The green vertical lines divide the peaks from the valleys. The NS phase 0.9 − 0.1 is shadowed in blue.
In the text

	Fig. 4. Normalized NS spin pulse shapes. The averaged pulses of the four distinct orbits included in this dataset are shown in the top panels, while the averaged pulses for the k-means classifications are displayed in the bottom panels. To compute the average pulse shape, all pulses were normalized beforehand to ensure that each contributed equally to the final result.
In the text

	Fig. 5. Detected spikes at their respective orbital phases. The radial axis represents the prominence. The size of the bubbles is proportional to their duration. The color gradient is used to visually support the prominence gradient.
In the text

Fig. 6. Prominence vs. duration of spikes detected in this dataset. The blue squares and red triangles represent the two separated clusters obtained by a Gaussian mixture algorithm. The leftmost panel provides the results obtained for all the observations combined, while the three panels on the right display the results for the three individual observations that produced the better results in number of spikes detected, and thus in clustering.

In the text

	Fig. 7. Average spike duration for each of the two clusters (with the upper panel displaying the short spikes and the lower panel displaying the long spikes) as predicted by the Gaussian mixture algorithm for each observation. The x-axis represents the brightness of each observation, measured as the error-weighted count average (whc).
In the text

	Fig. 8. Upper row (left to right): Spectra, model, and χ2 values for Bright, Intermediate, and Faint. Lower row (left to right): Spectra, model, and χ2 values for Absorbed, and all spectra displayed for comparison.
In the text

	Fig. 9. Emission lines observed in the low-energy range within the spectra Orb 32 and the absorbed spectra.
In the text

Fig. 10. Representation of the system. The dotted circle represents the NS orbit and the white circle the companion star. The red arrows represent the path that the emitting radiation transverses in its path to the observer. This path is different in length and wind density depending on the orbital phase in which the radiation was emitted, provided that the system has an inclination and eccentricity different from zero. The mass-loss rate and inclination could be inferred by the differences found within the NH in three different orbital phases, corresponding to Orb 1-1, Orb 31-1, and Orb 32-1 (see Tables 5 for the NH values and 3 for the orbital phases).

In the text

	Fig. 11. Density variation encountered by the NS as it travels around the eccentric orbit (red continuous line). The measured pulse fraction (Table 3) might be correlated with the density of the NS environment.
In the text

Fig. 12. Ionization map in the plane of the orbit: Shown (from left to right) are the Bright, Intermediate, Faint, and Absorbed spectra. The delimited regions represent the Fe Kα region, i.e., with an ionization parameter lower than 100 and higher than 25. We note that the orbital phase is ignored as the spectra are not phase-resolved. For visualization purposes, the orbital phase ϕ = 0.12 is used.

In the text

	Fig. 13. Fe Kα formation region (ionization parameter in the range 50 − 100) vs. the Fe Kα line area.
In the text

	Fig. 14. Light curves showing the decreasing pulse intensity trend (Orb 1-1, Orb 14-3, and Orb 31-1).
In the text