© The Authors 2025

1. Introduction

High-mass X-ray binary pulsars are systems that host a strongly magnetized neutron star (NS) that accretes material from a companion star. As the material falls towards the NS, it gets coupled with the magnetic field lines and is channeled into the polar caps. At low accretion rates (Ṁ ≲ 10¹⁷ g s⁻¹), the accreted material is stopped near the NS surface, forming a hotspot rather than an extended column. In this regime, the emission is primarily thermal and originates from the heated polar cap.

At high accretion rates (Ṁ ≳ 10¹⁷ g s⁻¹), the radiation pressure becomes significant, leading to the formation of a radiation-dominated shock above the NS surface. This shock decelerates the inflowing material, creating a column structure where matter settles more gradually onto the NS surface. The column height increases with the accretion rate, influencing the observed X-ray emission pattern (Basko & Sunyaev 1976). The accretion power is converted into radiation via thermal and nonthermal processes such as synchrotron, bremsstrahlung, and the inverse-Comptonization of all of the aforementioned photon fields. Since the rotational axis is usually not aligned with the magnetic one, the X-ray emission will appear to be pulsated to a distant observer.

X-ray spectral analysis is one of the standard tools to infer the properties of X-ray pulsars (XRPs). Besides the phase-averaged spectral analysis (averaging over an observational exposure time much longer than the spin period), we can also study phase-dependent spectral variations in the so-called phase-resolved (or spin-resolved) spectral analysis. To translate photon arrival times into phases (i.e, rotational phases of the NS, defined in the range from 0 to 1, where 0 corresponds to a reference phase such as the minimum of the pulse profile), we would usually adopt a timing solution correcting for the orbital parameters of the system (if known) and, when significantly detected, the spin derivative terms.

The shape of the pulse profile (PP) is the result of the combination of many fundamental properties of the NS (e.g., emission anisotropy and light bending) in a non-linear manner. The majority of the pulsed photons originate from the accretion flow near the magnetic poles of the pulsar. X-ray photons emitted from the surface are redistributed, after interacting with the relativistic particles of the accreted material, via the aforementioned nonthermal radiative processes. This interaction occurs in the presence of strong magnetic and gravitational fields of the XRP. The resulting PPs are thus energy-dependent and, in addition, luminosity-dependent (Wang 1981).

For low luminosities, the photons escape primarily along the magnetic axis (pencil beam). However it has been theorized that beyond a certain luminosity threshold, a radiative shock is formed, causing photons to escape from the sides of the accretion column (fan beam) rather than the upward direction, dramatically changing the PPs (Becker & Wolff 2022, and references therein).

XRPs often exhibit cyclotron-resonant scattering features (CRSFs) in their spectra, appearing as absorption-like features in the 10 − 100 keV range (see Staubert et al. 2019, for a review). These arise near the NS magnetic poles, where electrons in the intense magnetic field are confined to quantized Landau levels. The energy spacing of these levels is given by

$$\begin{array}{c}\hfill E_{\mathrm{cyc}}=\frac{n}{(1+z)}\frac{ħeB}{m_{e}c}\approx \frac{n}{(1+z)}11.6[\mathrm{keV}]\times B_{12},\end{array}$$(1)

where z is the gravitational redshift, n the Landau quantum state and B₁₂ the magnetic field strength in units of 10¹² Gauss. When photons at these energies interact with electrons, they undergo resonant scattering, leading to a reduction in observed flux at E_cyc. The CRSFs (and harmonics) provide the only direct method to estimate the surface magnetic field of a NS, as the line-forming region is believed to be very close to the stellar surface. Observations reveal a complex phenomenology: (i) the profile of the fundamental line is non-Gaussian; (ii) CRSFs vary with pulse phase (e.g., Cen X-3); (iii) their energy shifts with luminosity (both positively and negatively) (Staubert et al. 2019); and (iv) they correlate with the spectral continuum. These behaviors suggest intricate physical processes, with ongoing theoretical efforts seeking to explain their non-linearity and dependence on accretion dynamics (Poutanen et al. 2013; Loudas et al. 2024).

4U 1538–52 is a high-mass X-ray binary pulsar consisting of a NS in orbit with the B0Iab supergiant companion QV Nor. The source is known for its persistent X-ray emission, with notable flaring activity attributed to the close proximity and low eccentricity (e = 0.18) of the binary system (Hemphill et al. 2014). The orbital characteristics of 4U 1538–52 are well-defined, with an orbital period of approximately 3.728 days (Hemphill et al. 2019) and evidence of eclipses lasting about 0.2 days (Rodes-Roca et al. 2009). The orbital inclination is well constrained at 67°(1) (Falanga et al. 2015). The source shows also super-orbital modulation at a period of 14.9130 ± 0.0026 days, consistent with four times the orbital period (Corbet et al. 2021).

The pulsar has a spin period of approximately 530 seconds (Davison et al. 1977). The NS displayed many episodes of torque reversals (Malacaria et al. 2020), with the latest reversal from spin-down to spin-up at the end of 2019, although the general increase in frequency appears to be modulated by many intermittent, months-long spin-down episodes¹, indicative of the complex interactions between the NS and the wind of the companion (Karino 2020; Sharma et al. 2023; Hu et al. 2024). Timing studies have shown that the PPs remain invariant before and after torque reversals, suggesting a stable oblique rotator configuration (Hu et al. 2024).

Early studies of 4U 1538−522 X-ray spectra revealed the presence of cyclotron resonance scattering features (CRSFs) at around 20 keV (Clark et al. 1990). A possible first harmonic was reported by Robba et al. (2001) at 51 keV and by Rodes-Roca et al. (2009) at 47 keV. Hemphill et al. (2016) reported a secular change of the CRSF energy, approximately 1.5 keV over a span of 8.5 years. Varun et al. (2019) reported on the pulse-phase dependence of the spectral parameters using a 50 ks observation with ASTROSAT. They noted a maximum relative shift of the cyclotron line energy of 13% (E_cyc in the20–23 keV range), a variation of the CRSF depth correlated with the pulse shape and the spectral continuum.

In addition, spectral analysis has identified various emission lines, particularly in the iron Kα band. The presence of absorption features, such as one at 2.1 keV, has also been investigated, providing insights into the material surrounding the NS and its magnetic field (Rodes-Roca et al. 2014).

The study of energy dependence of the PPs together with the low energy position of the CRSF, allows us to set constraints on the geometry of the source. Clark et al. (1990) was the first to perform phase-resolved spectroscopy, revealing the variable nature of the CRSF and modeling the energy-dependent PPs using the framework proposed by Meszaros & Nagel (1985b). Later, Cominsky & Moraes (1991) obtained similar results to those of Clark et al. (1990), fitting the PPs with a model similar to that of Parmar et al. (1989) and finding evidence of a magnetic dipolar asymmetry. Further refinements to the accretion geometry were made by Bulik et al. (1992, 1995), who developed detailed models that fully accounted for both relativistic effects and the combined intensities of radiation and magnetic fields. Their findings suggested an asymmetric configuration of the magnetic polar caps, which were neither equal in size nor strictly antipodal.

Ferrigno et al. (2023) developed a pipeline tailored for PP analysis of NuSTAR data, which first showcased the utility of modelling, even phenomenologically, the pulsed fraction spectrum (PFS) to infer PP changes correlated with the CRSFs. Local changes in the PFS were observed at energies corresponding to features in the X-ray spectrum. Specifically, for all XRPs analyzed in that work (Her X-1, 4U 1626−67, Cen X-3, and Cep X-4), a dip in the PFS was detected at the energy of the cyclotron absorption feature. In contrast, the PFS of the XRP V 0332 + 53 (D’Aì et al. 2025) exhibits two distinct and narrow bumps that have been interpreted as indicative of a more complex CRSF profile with extended wings.

In this work, we aim to further investigate the connection between the PP properties and the spectral features of the source 4U 1538–52, using all the available NuSTAR data. The four observations have comparable luminosities, as reported in Table 1, allowing for a coherent analysis of the pulse phase–dependent behavior across epochs. We focused on achieving the best possible modeling of the PFS and find that it differs significantly from previous cases. Notably, we did not observe a dip, nor any narrow emission wing-like features at the cyclotron line energy. Instead, our modeling reveals a broad bump at energies around the spectral energy of the CRSF, highlighting a connection between this feature and the spin-dependent spectral variability of the source. Finally, we explored possible explanations for the diverse phenomenology observed in the PFS of different sources, particularly around CRSF energies.

2. Observations and data analysis

We used data from the NuSTAR satellite (Harrison et al. 2013) which collects X-ray events through two independent and similar detectors, the focal plane module A (FPMA) and B (FPMB). NuSTAR has observed 4U 1538–52 four times to date (see the observations log in Table 1). The first observation in August 2016 (30201028002) shows the system during a full eclipse passage. The second (30401025002) and third observations (30602024002), taken in May 2019 and in February 2021, respectively, show the source in the persistent state. In the final observation (30602024004), which occurred just a few days after the third, an eclipse was observed again during the second half of the observation. Hereafter, we will refer to the single observations using the label shortcuts 302, 304, 306-2, and 306-4.

We processed the data using the NuSTAR Data Analysis Software (nustardas, v. 1.9.7, as in HEASoft v. 6.31), and the most recent NuSTAR calibration files (CALDB v20240104). For higher level scientific products, we used our dedicated toolkit NUSTARPIPELINE². After re-processing the data to clean level 2 data stage, we defined a circular source region centered on the most accurate X-ray coordinates of the source, with a radius of 2 arc-minutes, and a background region of the same size in an area free of any contaminating sources for both the FPMA and FPMB detector images.

We excluded, for every single observation, time intervals where the source exhibited clear signs of dips, flares, and noticeable spectral changes, including the eclipse events in the observations 302 and 306-4. To do this, we defined the hardness ratio between the source light curves within the 3–7 keV and 7–30 keV energy bands, using 0.1-second time bins and then examined the time evolution. They were subsequently re-binned to maintain a minimum signal-to-noise ratio (S/N) of 10. If no major variations were detected, we focused on the flux changes in the re-binned light curves, constructing a histogram of the overall 3–30 keV rates. Assuming a log-normal distribution for the histogram, we fit a Gaussian in logarithmic space and filtered out intervals that deviated by more than 5σ from it (see Appendix D).

2.1. Timing and pulse profile analysis

For the timing and the PPs analysis, we used our pipeline, tailored for NuSTAR data, described in detail in Ferrigno et al. (2023). Here, we summarize the main steps of the analysis for every ObsID. First, the arrival times of each event are converted to the Solar System barycenter frame. We perform a Lomb-Scargle search to identify coherent signals at the light curve of the energy range of 3–7 keV, usually sufficient to find pulsations in case of a bright XRP, and define the preliminary spin frequency (P_spin). After verifying this value for consistency with the values reported in the literature, we refined the spin frequency by performing an epoch folding search within a small frequency interval around ±5% of the preliminary P_spin. We report the refined frequencies in Table 1. We also correct the arrival times for the binary motion of the NS by using the orbital parameters³ available from the Fermi Gamma-ray Burst Monitor (GBM) online database (Malacaria et al. 2020).

To study the PPs properties as a function of energy, we generated multiple energy–phase matrices with different numbers of phase bins for each observation, folding the data using the best period spin value. We computed separate source and background matrices for the two detectors, FPMA and FPMB, and then summed the resulting matrices. The energy bin spacing is set to match the intrinsic energy resolution of the FPM. For each energy bin, we calculated the S/N of the pulse while accounting for background events and summed the adjacent bins to reach a minimum S/N of 15 across the NuSTAR energy band. For a quick-look visualization of these matrices, we produced color-intensity heat maps of the energy-resolved profiles (Fig. 1).

Fig. 1. Energy-phase map of the observation 304, generated with 32 phase bins and a minimum S/N of 15: normalized PPs at equally spaced logarithmic energy intervals. Each pulse in the bin was normalized by subtracting its mean and dividing by the standard deviation.

We did not attempt to find a unique timing solution to phase-connect the pulses of the different observations, and, as no statistical evidence was found, no spin derivative term was taken into account for the single observation. However, the PPs of the different observations are very similar, regardless of the energy range from which they are extracted. We therefore assumed that the source phase-dependent properties also remain similar, when the profiles are correctly aligned.

To align the PPs of the four different ObsIDs, we took as a reference PP the ObsID 304, which is the one with the highest statistics. We extracted the energy-averaged 3–70 keV PP using 64 phase bins. For this PP, the zero phase corresponds to the center of the bin with the minimum normalized rate value. To align the phases of the other PPs, we calculated the sum of the squared differences (SSD) between the PP to be aligned and the reference one. For each of the 64 phase bins, we computed the squared difference between the respective values of the two PPs and sum these differences over all bins. We then shifted the profile by one phase bin and repeated the process, resulting in 64 SSD values for each observation. We visually identified the phase bin corresponding to the minimum value and refined the phase measurement value further: we modeled the points around the minimum with a constant plus a Gaussian line to reduce statistical fluctuations. The best-fitting centroid value of the Gaussian provides the phase shift, and by using this value, we could determine the new time of reference for the aligned profiles (see Fig. 2). The PPs of the different observations aligned according to this method, are shown in Fig. 3.

Fig. 2. Alignment of the energy-phase matrices. The plot shows the SSD values calculated for each phase shift between the PP of a given observation and the reference PP (ObsID 304). The minimum SSD value corresponds to the best alignment, and the points around this minimum are fitted with a Gaussian to refine the phase shift determination and account for statistical fluctuations. The centroid of the Gaussian provides the precise phase shift used for alignment.

Fig. 3. Aligned PPs in the 3–70 keV band (errors not shown for clarity’s sake). Two periods are shown for better clarity.

We pass from a qualitative description of pulse-energy dependence to a more quantitative analysis using the PFSs, following the methods described in Ferrigno et al. (2023). In brief, to calculate the pulsed fraction (PF) value of a single energy bin, we applied a fast Fourier transform (FFT) on each energy-resolved pulsed profile, defined by the re-binned energy-phase matrix. The PF value is derived using a harmonic decomposition, given by the equation

$$\begin{array}{c}\hfill {\mathrm{PF}}_{\mathrm{FFT}}=\frac{\sqrt{2{\displaystyle \sum _{k=1}^{N_{\mathrm{harm}}}}{|A_k|}^2}}{|A_0|},\end{array}$$(2)

where |A₀| is the average value of the PP (the zero term of the FFT transform), and |A_k| represents the amplitude of the k-th harmonic. We extended the number of harmonics necessary to describe the pulse with a statistical acceptance level of at least 10%, which typically corresponds to three to five harmonics. Even though other definitions exist for calculating PF values, such as the min-max method, which relies on the relative difference between the maximum and minimum flux of the PP, Ferrigno et al. (2023) demonstrated that the FFT method is the most robust and least biased with respect to the number of phase bins used.

The uncertainty for each PF point is assigned with a bootstrap method. We generated 1000 faked profiles, assuming Poisson statistics for the counts registered in each phase bin, we calculated the PF value for the generic faked profile, and from this Monte-Carlo sample, we derived the standard deviation and assumed it as the uncertainty at 1-σ confidence level.

We then modeled the broad energy dependence of the PFS using a phenomenological polynomial function. To keep the polynomial order as low as possible, we split the full dataset into a low-energy and a high-energy band. The splitting energy was treated as a free parameter within the 10–15 keV range and was chosen as the energy corresponding to a zero in the first derivative of an interpolating function (Ferrigno et al. 2023). The polynomial order is the minimum one to ensure a p-value greater than 0.05. We used Gaussian profiles to model local features in both the low- and high-energy bands.

We use the lmfit package (v1.3.2) for the PFS fitting, with parameter uncertainties estimated via Markov chain Monte Carlo (MCMC) sampling using the emcee package (v3.1.4). The algorithm from Goodman & Weare (2010) was employed with 50 walkers, 600 burn-in steps, and 6000 total steps. The best-fit parameters were derived from the 50th percentile, with uncertainties from the 16th and 84th percentiles. Optimal parameters, such as an S/N of 15 and 32 phase bins, were chosen for this source to balance spectral resolution and S/N. Lastly, we compute the correlation and lag spectra for each observation using the methods described in Ferrigno et al. (2023). We used the 3–70 keV and the 35–70 keV as the PPs of reference, but excluding the energy bin of the PP we want to cross-correlate when carrying out this template PP computation.

2.2. Spectral analysis

We generated the ancillary response and redistribution matrix files using numkarf and numkrmf, respectively. The source spectra were processed using the ftgrouppha tool⁴, following the optimization algorithm described by Kaastra & Bleeker (2016) and maintaining a minimum S/N of 5 for ObsID 304 and a minimum S/N of 3 for the ObsID 306-2. We only analyze the spectra of ObsIDs 304 and 306-2, as these observations provide the longest time windows when the source is in a persistent state (no eclipses), thereby giving sufficient statistics for also conducting a constraining phase-resolved analysis.

The spectral analysis results from these NuSTAR observations have been previously reported in the literature (Hemphill et al. 2019; Sharma et al. 2023; Hu et al. 2024; Tamang et al. 2024). Therefore, rather than attempting to determine the best-fit model for the entire spectrum, our primary focus is on determining the CRSF parameters while ensuring they remain uncorrelated with the continuum parameters. To achieve this, we performed a Bayesian fit on the re-binned spectra in the 4–50 keV energy range, determined by the quality of the data.

We first performed a fit using a standard minimization algorithm, then explored the neighboring parameter space using a MCMC method. Specifically, we employed the Goodman & Weare algorithm (Goodman & Weare 2010) with 200 walkers, a burn-in phase of 20 000 steps, and a total length of 200 000 steps. We used the W statistic, which is analogous to the Cash statistic (C-stat) when a background dataset is provided, assuming Poisson distributions for both source and background spectra (see Sect. 7 in Arnaud et al. 2011). To ensure the robustness of our results, we visually inspected the fit statistic distribution and parameter corner plots to verify chain convergence and confirm that the cyclotron line parameters remain uncorrelated with those of the continuum. All spectral fitting is performed using xspec.

To account for systematic differences between the FPMA and FPMB instruments of NuSTAR, we included a multiplicative constant in the model. This constant was fixed at 1 for the FPMA data, while it was left as a free parameter for the FPMB data to allow for calibration adjustments between the two detectors. In our fits, the constant for the FPMB remained close to 1, indicating minimal cross-calibration discrepancies between the two instruments and ensuring a consistent flux scale across both modules.

3. Results of the timing and pulse profile analysis

The PPs of all four observations display similar behavior, as illustrated in Fig. 3. The 3.0–70.0 keV energy-averaged PP for all ObsIDs reveals an asymmetric double-peaked shape. An exception is observed in the first peak of ObsID 304, which appears more structured, while the peak’s width remains similar across the other three observations. To analyze the energy dependence of the PPs, we generated energy-resolved PPs in the 3–10, 10–19, 19–24, 24–30, 30–40, and 40–70 keV energy ranges (see Fig. 4).

Fig. 4. Energy-selected PPs for all examined observations.

The amplitude of the secondary peak is strongly energy-dependent: it is aptly detected below 19 keV, nearly undetectable in the 19–24 keV range, reappears in the 24–30 keV range, and then vanishes at higher energies. This is consistent with the results reported by Clark et al. (1990). Notably, the 19–24 keV energy range corresponds to the CRSF energy, whose centroid energy is approximately 21 keV. This trend is consistent across all four observations, even though in the 302 observation the reappearance of the second peak is barely detectable. No significant phase shift of the secondary peak is observed. In contrast, the amplitude of the first peak shows less energy dependence, although its position shifts to earlier phases as the energy increases. These broad-band, energy-resolved PPs clearly demonstrate the strong dependence of the PPs on the spectral shape.

3.1. Pulsed fraction spectral fit results

We carried out a phenomenological fit of the PFS following the methodology described in Section 2.1, dividing the fitted range into a soft and a hard band split at around 11 keV. In all four observations, the PFS exhibit a similar trend. The PF values begin at approximately 0.3 in the low-energy band (around 3 keV) and gradually increase up to 0.8 at 19.5−20 keV. Subsequently, the PF values decrease, only to rise again after 30 keV, reaching again values close to or higher than 0.8 in the highest energy band (above 40 keV). Notably, there is a clear dip in the PF values at 6.4 keV, the energy corresponding to the Kα iron line (Rodes-Roca et al. 2014). This dip is clearly detected in the observations with the most robust statistics, 304 and 306-2, and was modeled using a negative Gaussian profile, while in the other two observations, 302 an 306-4, there are no significant residuals in this energy range.

To model the PFS, we explored two different phenomenological approaches: fitting a polynomial combined with a positive broad Gaussian centered at approximately 20 keV, and a negative broad Gaussian centered around 25 keV. In all four PFS examined, the model incorporating the positive Gaussian consistently provided a better fit than the one with the negative Gaussian (see Fig. 5 and Appendix A). The corresponding fit results are summarized in Table 2. Across all fits, the centroid energy of the positive Gaussian lies between 19.5 and 20.5 keV, approximately 1 keV below the cyclotron line energy in the corresponding spectra (Hu et al. 2024). The Gaussian width is found to be between 2.2 and 2.9 keV, which is comparable to the observed width of the CRSF, of around 3 keV, thus suggesting a close link between this feature in the PFS and the CRSF. Given this result, we next sought to investigate the physical origin of this bump in the PFS and its potential connection to the cyclotron line.

Fig. 5. Summary plot for ObsID 304 using 32 phase bins and a S/N 15. Panel a: PF (green points) along with its best-fit model (solid lines); polynomial fits are also included for comparison. Panel b: residuals from the fit. Panels c–f: Phases and normalized amplitudes of the first (A1, ϕ1) and second harmonic (A2, ϕ2), respectively. The vertical shaded regions denote the energy ranges and widths of the Gaussian functions fitted to the PF. Panel g: Selection of normalized PPs at equally spaced logarithmic energy intervals, horizontally offset for clarity. Each pulse in the bin was normalized by subtracting its mean and dividing by the standard deviation. Panel h: Color-map display of the normalized PPs as a function of energy, with 20 evenly spaced contour lines drawn for reference. Panel i: CC between the PP at each energy band and the overall average profile. Panel j: corresponding phase lags. The colored vertical regions are consistent with those in panels d–f.

3.2. Harmonic decomposition: Energy dependence of the amplitudes

The contribution of the first two harmonics of the Fourier decomposition dominate the 3–70 keV energy averaged PP in all the four observations, as shown in Fig. 6 with the truncated decomposition of ObsID 304. The energy-dependent behavior of the first two harmonics of the FFT is similar across all four observations. The amplitude of the first harmonic (A₁), computed as the Fourier amplitude normalized to the mean count rate of the PP, follows the general trend of the PFS, showing increased values around 20 keV, followed by a decrease and then an increase at energies above 30 keV. In addition, a decrease is observed at the iron line energy in the two observations with the highest statistics, 304-2 and 306. The second harmonic amplitude (A₂) also shows a increase, peaking at energies just below 20 keV, which is lower than the broad feature in the PFS. However, this rise is more gradual and less pronounced compared to the first harmonic. At energies corresponding to the broad feature in the PFS, the second harmonic shows a slight dip for a few keV before briefly rising again (see Fig. 5 and Appendix A).

Fig. 6. Total 3–70 keV PP of ObsID 304 and the first two Fourier harmonics from its FFT decomposition.

3.3. Cross-correlation and lag spectra

We computed the energy-resolved cross-correlation (CC) spectra for all four observations to evaluate the energy-dependent coherence of the PPs with respect to reference PPs. In the summary plots (see Fig. 5 and Appendix A), we show the CC and lag spectra using the 3−70 keV averaged PP as template.

The behavior across all observations exhibits similar patterns. Specifically, the CC shows a steady increase up to approximately 10 keV, followed by a decrease until 20 keV, which aligns with the broad feature observed in the PFS. Beyond 20 keV, the values rise again until around 30 keV. After this point, in observations 306-2 and 306-4, the data points become too sparse to discern any clear trend. In the other two observations, the CC values decrease.

In the lag spectra at the cyclotron line energy, we observe that global trend of the lag values drastically changes, from a relatively flat trend to a steep decline of the lag values towards earlier phases.

We also found the CC spectrum to be particularly informative when using the PP of energies above 35 keV as a template, as shown in Fig. 7. In this case, the CC spectrum more clearly highlights the strong correlation between the PP extracted around the CRSF energy range, where the bump on the PFS is observed, and the highest energy band.

Fig. 7. CC spectrum extracted during ObsID 304 using the PP from the 35–70 keV energy range as template.

4. Results of spectral analysis

To model the broadband continuum emission we used the newhcut component, first introduced by Burderi et al. (2000). This model represents a power-law with a high-energy cutoff, where a smoothing function in the form of a polynomial, is applied near the cutoff energy. The width of the polynomial, which determines the energy range over which the transition from the power-law to the cutoff occurs, is typically fixed at 5 keV. Preliminary fit results indicated that this component alone was insufficient to fit the broadband 4–50 keV spectrum; therefore we added to the model a soft blackbody component (bbody), which is often required in modeling the spectra of bright accreting XRPs as also done in the literature for this source in particular Sharma et al. (2023). Since our analysis is limited to energies above 4 keV, we fixed the spectral parameters that are less constrained due to the lack of the softer spectral coverage, namely, the black-body temperature at 1 keV and the value of the equivalent hydrogen absorption column (modeled using the tbabs component in xspec), N_H = 1 × 10²² cm⁻², both values taken from the fits in Sharma et al. (2023).

We find that the CRSF shape is statistically better modeled using a Lorentzian profile (cyclabs in xspec) rather than a Gaussian profile (gabs in xspec). The latter was more susceptible to spectral degeneracy between parameters governing the line shape and the continuum in some of the phase-resolved spectra. The iron line complex was modeled using a single Gaussian emission line (gauss in xspec).

As reported in Sect. 1, a CRSF harmonic has been reported for this source by several authors (Robba et al. 2001; Hemphill et al. 2019), although we did not find strong evidence of it in the two examined phase-averaged spectra. Although the addition of an absorption feature at the spectral model, at around 50 keV does improve the fit, its energy position is poorly constrained while its width takes large values greater than 10 keV. This suggests that the feature likely models residuals in the continuum rather than representing an absorption feature at these energies. The fits are shown in Fig. 8 for ObsID 304 and in Fig. B.1 for ObsID 306-2. The results of these phase-averaged fits are reported in Table 3.

Fig. 8. (a) Data and best-fitting model for the phase-averaged spectrum of ObsID 304. (b) Residuals of the best-fitting model when the depth of the cyclotron line is set to zero. (c) Residuals of the best-fitting model.

To investigate the phase-dependent spectral variations of the cyclotron line, we performed an independent phase-resolved spectral analysis for the two observations, 304 and 306-2. We found that a satisfactory compromise between spectral resolution and statistics is obtained by using eight evenly spaced phase bins. This also ensured a good coverage for energies up to 40 keV for each phase-selected spectrum. The phase bins are defined as shown in Fig. 9 and will be referenced accordingly from now on.

Fig. 9. Definition and naming of the phase bins used in the analysis.

We used the same spectral model of the phase-averaged analysis to fit the phase-resolved spectra. The variability of the best-fitting parameters is shown in Fig. B.4. The plot of all best-fit models from the phase-resolved analysis (normalized for clarity to unity at 10 keV) can be seen at Fig. 10, which illustrates the variability of both the continuum and the cyclotron line.

Fig. 10. Phase-dependent spectral variability for ObsID 304: each curve represents one of the best-fitting model from the eight phase-selected spectra. Models are normalized to unity at the reference energy of 10 keV. The blue and red curves correspond to spectra extracted in the phase intervals of the primary and secondary peak, respectively. The black spectrum shows a very shallow cyclotron line in the primary falling flank. The gray spectra correspond to the primary rising flank and secondary falling flank.

The continuum parameters exhibit significant phase variability, with similar trends seen for both observations. The photon index is higher in the primary rising and secondary falling flanks, whereas the folding energy decreases before starting to rise again in the secondary falling flank. The normalization of the blackbody component closely follows the shape of the energy-averaged PP.

The only phase bin with significantly different parameter values in the two observations is the primary falling flank. For instance, in observation 304, the flux of the blackbody component is the highest among all phase bins, while in observation 306-2, it is the lowest and constrained only by an upper limit. Conversely, the normalization of the cutoff power law shows the opposite trend. Additionally, satisfactory fits are generally obtained when the temperature of the blackbody component is fixed at 1 keV, while for just this phase bin we had to set it to 1.5 keV in order to achieve a good fit. These findings suggest stronger degeneracy and inter-correlation of the parameters for this specific phase bin. Regarding the CRSF in both observations, it is only marginally detected in this particular phase bin, whereas it is detected in all the other bin intervals.

Although the Lorentzian profile provides a better fit for the CRSF shape compared to the Gaussian absorption profile, it remains poorly constrained in certain phase-resolved spectra. To restrict the parameter space degeneracy, we chose to freeze the width of the cyclotron line to the phase-averaged value when the fit tended to converge to nonphysical, very broad values, such as in the primary rising flank, in the secondary falling flank, and (just in the case of ObsID 306-2) in the primary peak A.

In the primary falling flank, we found that the CRSF is only marginally detected. In ObsID 304, although the fit could constrain a value for the line width, it had a large uncertainty. In contrast, in ObsID 306-2, where the line depth is even lower, the line width had to be frozen to the value determined for the same phase bin in the ObsID 304 spectrum to stabilize the fit.

The general variability of the CRSF between the primary and secondary peaks is arguably the most intriguing observational fact, particularly the modulation of the line’s depth. This is clearly shown in Fig. 10, where the best-fitting models of the phase bins of the primary and secondary peak are colored in blue and red respectively, to better highlight the different shapes.

The CRSF energy, E_cyc, exhibits significant phase-dependent shifts. The highest values, reaching up to ≈22 keV, are observed during the primary peak phases in both observations, while the lowest values, ranging between ≈19.1 and ≈19.5 keV, occur in the secondary falling flank, approximately 3 keV lower than in the primary peak. The line depth is significantly greater in the secondary peak than in the primary peak. The difference in energy and depth values between the two peaks, D_cyc, is illustrated in Fig. 11, where it can be seen that the values corresponding to the two pulse peaks tend to cluster in opposite regions of the diagram.

Fig. 11. CRSF energy versus line depth values for the phase bins corresponding to the primary and secondary peaks (0.250–0.500 and 0.625–0.875).

5. Discussion

The relation between the CRSF and the energy resolved PPs of 4U 1538–52 has been extensively studied in the past. Clark et al. (1990) was the first to not only detect the CRSF in the X-ray spectrum but also to highlight its phase-dependent variability using data from Ginga, along with an early attempt to physically model the PP shape changes with energy. A more detailed analysis was done by Bulik et al. (1995), where they showcased the impact of general relativistic effects on the phase-dependent spectral features of 4U 1538–52. By employing model calculations of inhomogeneous magnetized NS atmospheres and fitting Ginga data, they determined geometric and magnetic parameters of the system. Their results suggested that the polar caps are unequal and non-antipodal, indicating either a misaligned or distorted magnetic field structure. Additionally, the study found that relativistic effects significantly influence the observed flux and PPs, with evidence for a bent or off-center magnetic axis. A similar conclusion was reached by Burderi et al. (2000) for Cen X-3, reinforcing the idea that deviations from a simple dipolar geometry may be a common feature in accreting pulsars.

In this work, we analyze the PPs of 4U 1538–52 using data from NuSTAR, which provides improved sensitivity and broadband coverage in the hard X-ray band. With its high spectral resolution, NuSTAR allows for a more detailed investigation of the PP morphology and the spectral phase-dependent behavior of the CRSF. We examined the persistent emission of the source (i.e. outside of the eclipse) across the four NuSTAR observations. No significant variation in the accretion regime was noted (the average source flux varied by a factor of less than two among the observations and the spectral shape remained similar) and the analysis of pulse-phase dependent behavior across epochs produced consistent results. We focused on how the CRSF characterizes and modifies the energy-resolved PPs, by extracting the analytical PPs elements from the energy-phase matrix from the available NuSTAR observations. In particular, we emphasize the importance of modeling the PFS to gain insights into the spin dependence of the CRSF parameters. We found consistent results for all the examined observations, except for eclipse intervals where the low S/N prevented a detailed investigation. By generating energy-resolved PPs across different energy bands, we found that the secondary peak of the PP nearly disappears at the CRSF energy and then reappears at higher energies, consistent with the findings of Clark et al. (1990). In the PFS, we observed a significant increase of the PF values just before the CRSF energy, followed by downtrend after the CRSF energy. We modeled the PFS phenomenologically using a polynomial, a negative Gaussian to the iron line energy and a positive, moderately broad, Gaussian in the energy range of the CRSF. The best-fitting parameters of the Gaussian resulted in a width ranging from 2.2 to 2.9 keV across the four analyzed ObsIDs. Such values are comparable to the width of the CRSF found spectrally. The centroid energy, around 20 keV, is 1 keV lower than the spectral value (but we note that the spectral centroid energy systematically depends on the functional shape used to fit the CRSF, namely Lorentzian profiles tend to downshift measured energies with respect to Gaussian profiles). Furthermore, when using the energy-resolved PP above 35 keV as a template for constructing the CC spectrum, we observed a maximum at the CRSF energy, consistent with the disappearance of the secondary peak at these energies. Additionally, phase-resolved spectroscopy revealed that the cyclotron line depth differs significantly between the primary and secondary peak: in the phase interval of the secondary peak, it is much greater than in the phase interval of the primary peak. The centroid energy also differs, ranging from around 20.7 to nearly 22 keV at the primary peak and from around 20.4 to 20.7 keV at the secondary peak. The cutoff and e-fold energies of the continuum model also vary significantly between the two peaks, as the spectrum of the secondary peak drops off more rapidly above 35 keV compared to that of the primary peak, consistent with the PP analysis.

Fig. 12. PFS of observation 304, calculated using 8 phase bins and an S/N of 5 for consistency with the phase-resolved spectral analysis, overplotted with the primary (blue) and secondary (red) peak spectra.

5.1. Explaining the dips and bumps at the CRSF energy in the pulsed fraction spectra

The combined analysis of the PFS, PPs, and phase-resolved spectroscopy provides clear hints that the bump in the PFS of 4U 1538–52 arises from the strong spin-phase variability of the CRSF. This is illustrated in Fig. 12, where we present thebest-fitting models of the spectra of the primary (0.250–0.500) and secondary (0.625–0.875) peaks, over-plotted with the PFS. The PFS is calculated using eight phase bins and a S/N of five to maintain consistency with the phase-resolved spectroscopy.

A key question that arises is why the PFS at the cyclotron energy deviates from the trends reported in the literature for other sources. Unlike the dips noted by Ferrigno et al. (2023) or the two narrow bumps observed adjacent to the cyclotron line in V0332+53, as described by D’Aì et al. (2025) – interpreted as cyclotron line wings – 4U 1538–52 does not follow either of these patterns. For this discussion, we aim to recreate the conditions that lead to a clear bump, or a dip, through a simple phase-dependent modeling of the energy spectra.

We first attempted to reproduce the PFS of 4U 1538–52, with a simple calculation. We assumed that we have two spectra at different flux levels, described by a simple power-law and an energy cutoff, and we insert an absorption feature on these two spectra, as shown in Fig. 13 (left). To roughly mimic the observed behavior of 4U 1538–52, we inserted an absorption feature at 20 keV (width: 3 keV) into the lower-flux (secondary peak) spectrum, and another absorption feature at 22 keV (width: 4 keV) into the higher-flux (primary peak) spectrum. Additionally, the primary peak spectrum is assigned a higher energy cutoff (12 keV) compared to the secondary peak spectrum (8 keV). We simulated multiple cases by varying the depth of the absorption feature in the primary peak.

Fig. 13. (a) Simulated spectra for the primary and secondary peaks, modeled by a power-law continuum with an exponential energy cutoff. The cutoff energy is set at 12 keV for the primary peak spectrum (blue) and 8 keV for the secondary peak spectrum (red). The absorption line in the secondary peak spectrum is centered at 20 keV with a width of 3 keV, while the absorption line in the primary peak spectrum is at 22 keV with a width of 4 keV. (b) Simulated PFS of the spectra shown in (a) using the CV proxy. The red and blue vertical lines highlight the energies of the absorption feature in the secondary and primary peaks, at 20 and 22 keV, respectively. (c) Simulated spectra for the primary and secondary peaks, where the energies of the absorption features have been exchanged between the two. (d) Simulated PFS of the spectra shown in (c) using the CV proxy.

To quantify the PFS for each case, we adopted the coefficient of variation (CV) as a proxy for the value of the PF value of each energy bin, defined as CV = σ/μ, where σ and μ are the standard deviation and the mean value of the points in an energy bin. The CV is thus roughly equivalent to the root-mean-square (RMS) calculation method, which, as demonstrated by Ferrigno et al. (2023), is analogous to the FFT approach. For the purposes of this exercise, we do not take into account model uncertainties. The results of this simple exercise are illustrated in Fig. 13.

In Figs. 13a and b, we observe that when the absorption line in the primary peak is very shallow, a distinct bump appears around the line energy, very similar to what we see on the PFS of 4U 1538–52. However, as the depth of the line increases across the different cases (in Fig. 13, from darker to lighter blue colors), the bump becomes less pronounced, its peak shifts to lower energies, and a dip emerges at higher energies. As the absorption line in the spectrum of the primary peak becomes deeper than that of the spectrum in the secondary peak, the dip becomes more prominent and the bump becomes progressively suppressed. On the other hand, when the energy of the cyclotron line in the spectrum of the primary peak is lower than that of the secondary peak spectrum, the dip and bump would appear in opposite locations, with the dip occurring at lower energies and the bump at higher energies, as shown in Figs. 13c and d. Therefore, modeling the PFS phenomenologically with only negative or positive Gaussian features may work only as a good first-order estimation, as the observed profile can get more complicated.

Additionally, we observe that the higher cutoff energy in the primary peak spectrum relative to the secondary peak spectrum leads to an upward trend in the PFS. This is because thespectrum of the secondary peak falls off more rapidly at high energies, thus causing the secondary peak to disappear in the PP, which increases the spectral variability between the two spectra at higher energies, and therefore the PF values.

In the case of 4U 1538–52, the fit described in Section 3 performs reasonably well, with reduced χ² values ranging from 1 to 1.5 across all four ObsIDs. We attempted alternative fitting strategies for the PFS, including a linear or second-degree polynomial continuum combined with a positive Gaussian near ≈20 keV and an additional negative Gaussian at ≈25 keV. However, these models resulted in poorer fit statistics, indicating that the additional negative Gaussian did not lead to a better description of the feature.

It is important to note that the above exercise, based on only two spectra (primary and secondary peaks), is a simplification. The total spin-dependent spectral variability of an XRP cannot be reduced to just two spectra, and the resulting PFS is more complex than what was captured in our simplified two-component case. Nevertheless, this exercise is still effective in illustrating how phase-dependent features such as bumps and dips emerge in the PFS, and in explaining why the energy at which they appear may not exactly coincide with the spectral energy of the cyclotron line. This framework provides a useful reference, especially for relatively straightforward cases such as 4U 1538–52.

Another simple case is Her X-1, one of the sources analyzed by Ferrigno et al. (2023). In this source, a pronounced dip is observed in the PFS near the CRSF energy. Phase-resolved spectroscopy by Fürst et al. (2013) revealed that the line’s depth varies strongly with pulse phase. In the off-pulse phases, the CRSF is shallow and centered at ∼35 keV, while at the pulse peak, its depth increases by a factor of about three, and its centroid shifts to ∼38 keV. This behavior is opposite to what is observed in 4U 1538–52, where the CRSF is deeper at low flux and shallower at high flux, namely, the correlation between spin-dependent flux and line depth is reversed.

Based on this simple exercise, we would thus expect the PFS of Her X-1 to show a more prominent dip rather than a bump, and for the dip to appear at higher energies than the spectral value of the CRSF. This is indeed consistent with the findings of Ferrigno et al. (2023), who report a PFS dip at 40.44 ± 0.15 keV, compared to the CRSF centroid of 37.4 ± 0.2 keV measured by Fürst et al. (2013) using the same dataset.

Consequently, the PFS effectively serves as an indicator of underlying spectral variability and provides additional constraints on the CRSF, which might otherwise be missed due to spectral degeneracy between the continuum and CRSF parameters. A systematic study of the different behaviors observed in the PFS and the phase-resolved CRSF spectral variability goes beyond the scope of the present work, though we have already identified and demonstrated the key drivers of the apparent diversity in PF trends around the cyclotron line energies. Additionally, we have shown that the upward trend of the PFS can be attributed to the lower high-cut energy in the spectrum of the secondary peak compared to that of the primary–or, equivalently, to the disappearance of the secondary peak at higher energies.

5.2. Physical modeling simulations

To understand the possible origin of the reported phase-energy behavior of the flux of the source, we perform physical modeling of emission from magnetic poles of a NS and qualitative comparison with observations. Here, we focus on general trends in energy-dependent behavior of PPs, variability of the fundamental cyclotron resonance, and the PFS around these energies. Specifically, we aim to address the formation of the two separate peaks in the low-energy (or energy-integrate) PPs, the reduction of the secondary peak near the observed CRSF centroid energy, the general appearance of the fundamental CRSF between the main and the secondary peak, and the resulting bump in the PFS.

We follow the same setup as presented in Meszaros & Nagel (1985a), performing polarized radiative transfer calculations in a homogeneous, strongly magnetized slab-like atmosphere. The atmosphere is self-emitting, with magnetic Compton scattering and bremsstrahlung, including the first cyclotron resonance and the first-order relativistic corrections. For these simulations, we use the FINRAD code (Sokolova-Lapa et al. 2021, 2023). Within this model, the set of parameters determining the emission is the following: the cyclotron energy, corresponding to the magnetic field, E_cyc, the electron plasma temperature, kT_e, the electron number density, n_e, and the total Thomson optical depth, τ_T. We obtain similar emission patterns to those described in Meszaros & Nagel (1985b): a slab of the depth τ_T ∼ 20 is transparent close to the magnetic field direction (θ ≲ 10°) for photons of soft energies, 2 − 6 keV, but has much higher opacity in this direction for photons of higher energies. At the same time the polarization-averaged differential flux dF = I(θ, E) cos θ falls down towards 90° more rapidly in the vicinity of the cyclotron energy than in the low-energy continuum.

The predicted dF allows us to set the emission from each surface element of two spot-like regions on the surface of the NS, associated with its magnetic poles. To obtain the projection of the emission onto the plane of a remote observer, we use the ray-tracing code by Falkner et al. (2024), which calculates photon trajectories from a rotating NS in the Schwarzschild metric. The observed flux depends on the locations of the poles with respect to the rotation axis, Θ₁ and Θ₂, the azimuthal separation of the poles, ΔΦ, and the inclination of the rotation axis to the observer’s line of sight, i. These parameters are referred to in the following as “geometric” ones. Here, we consider only the standard values for the NS mass, M_NS = 1.4 M_⊙ and radius, R_NS = 12 km, resulting in the surface gravitational redshift of 1 + z ≈ 1.24.

The preliminary study showed that for a fixed E_cyc, the main driver of the observed phase and energy variability is geometric parameters. We fix E_cyc = 25 keV, which results in observed redshifted cyclotron energy of ≈20 keV. In addition, we set the other parameters for the internal emission: kT_e = 6 keV, n_e = 10²³ cm⁻³, and τ_T = 20. Within the framework of the adopted model, the energy-integrated PPs, corresponding to the observed two clearly separated peaks of different heights (see, e.g., Fig. 3), occur only in a localized region of the parameter space, requiring i ≳ 45° and 45° ≲Θ₁ ≲ 90° simultaneously. At the same time, only geometries with |i − Θ₁|≲25° and Θ₁ ≲ 85° result in a significant reduction of the secondary peak in PPs obtained near the fundamental CRSF energy. Consistent with previous studies (Bulik et al. 1992, 1995), we conclude that the positions of the two peaks in the PPs of 4U 1538–52 require a slight asymmetry in the pole locations of ≈10 − 15°. Here, we qualitatively reproduce some of the observables and illustrate them with a showcase below. A more detailed exploration of the parameter space will be presented in future work focused on the physical model.

Our simulations indicate that, given a small size of the poles on the NS surface (with a radius of about a few hundred meters), the geometric parameters and the internal emission pattern, I(θ), determine the observed spectral flux variations. The angle, under which the observer sees each of the poles at a certain rotation phase, is fully determined by the geometry and the compactness of the NS. For a selected geometry, which, according to our simulations, is appropriate for 4U 1538–52, i = 67° (fixed to the inclination of the system), Θ₁ = 80°, Θ₂ = 105° and ΔΦ = 168°, we show the phase dependency of the emission angle in Fig. 15a and the adopted geometry in Fig. 14. This variation is a crucial component for understanding the observed cyclotron line behavior. The scattering cross sections exhibit strong anisotropy dominated by angle-dependent Doppler broadening near the cyclotron resonance (see, e.g., Schwarm et al. 2017). As a result, the closer the emission angle is to the magnetic field direction, along which the thermal motion of electrons is not restricted, the shallower the line profile observed (Meszaros & Nagel 1985b).

Fig. 14. Geometry adopted for 4U 1538–52, showing the NS, magnetic poles at Θ1 = 80° and Θ2 = 105°, azimuthal separation ΔΦ = 168°, and inclination to the observer’s line of sight i = 67°.

Fig. 15. Phase-dependent emission angle and various observables modeled for a selected geometry: i = 67°, Θ1 = 80°, Θ2 = 105° and ΔΦ = 168°. The mass and radius of the NS are MNS = 1.4 M⊙ and RNS = 12 km, the size of spot-like poles is r0 = 300 m. The emission from each emitting surface element is calculated for: Ecyc = 25 keV, kTe = 6 keV, ne = 1023 cm−3 and τT = 20. (a) Emission angle, visible for a remote observer at each rotation phase. The angle is calculated for a center of a pole on the NS surface in the Schwarzschild metric (blue dashed-dotted line corresponds to pole 1, red dashed line – to pole 2). The shadowed region shows when a pole is obscured by the NS. (b) Energy-integrated (5 − 80 keV) pulsed profile (black thick line) for the same geometry, with individual pole contributions. Light-gray line shows the PP in the energy range 19 − 24 keV, corresponding to the core of the fundamental CRSF in the observer’s rest frame. (c) Modeled energy flux near the fundamental line. The blue and red lines are obtained for the central regions of the primary (“peak 1”) and secondary (“peak 2”) peaks in the PPs; the respective ranges of phases are marked by shadow regions in panel b. The black dotted line marks the redshifted value for the cyclotron energy. (d) PFS in the same energy range as in panel c.

Figure 15b shows PPs for the chosen geometry. As indicated by the dashed-dotted blue and dashed red lines, each of two peaks in the profile is fully dominated by the contribution from one of the poles. From Fig. 15a, it is clear that near the first peak, we observe photons emitted from pole 1 at an angle of ≈10° to the magnetic field, while near the secondary peak, photons are emitted from pole 2 at ≈30°. This difference in the emission angle results in a significantly more shallow profile of the fundamental CRSF observed at phases near the main peak (Fig. 15c). As was shown in Sect. 5.1, this fact leads to a pronounced bump in the PFS, which we show in Fig. 15d. The PF peaks at slightly lower energies than the geometric center of the observed lines, indicating the largest difference between the line fluxes. Due to the complex redistribution occurring during the formation of broad lines near the continuum cutoff, one can expect some variation of the centroid energies of the lines seen under different angles. Specifically, a simple fit with a Gaussian absorption in the two spectra of Fig. 15c gives a centroid energy of 24.5 ± 0.2 keV for the spectrum of the primary peak, and 22.7 ± 0.4 keV for that of the secondary peak. The bump appears at lower energies, consistent with the expectations discussed in Section 5.1.

A similar physical model for the emission was previously adopted to describe spectra (Bulik et al. 1992) and PPs (Bulik et al. 1995) of 4U 1538–52. Bulik et al. (1992) included the fully relativistic Compton scattering cross sections and obtained observed spectral fluxes for different rotation phases. They considered a more complex shape of emission regions – spherically curved polar caps, but omitted the general relativistic effects on photon trajectories. The latter effect was included in Bulik et al. (1995) who fitted PPs of the source at different energies. The emission was calculated from an inhomogeneous atmosphere, but with the same non-relativistic cross sections as used in the current work. We would like to note that in the qualitative study performed here, the light-bending effect plays a more crucial role in shaping the general behavior of observables than the relativistic effects in the scattering cross sections. The latter is not expected to dramatically alter the emission pattern from magnetized plasma but does noticeably affect the width of the cyclotron resonances, generally making them narrower and deeper (Alexander et al. 1989; Bulik et al. 1992). Light bending plays a crucial role in the observability of emitting regions on (and near) the neutron surface, causing radiation beaming even from an isotropically emitting spot and significantly increasing the contribution from the back side of the star (Pechenick et al. 1983). In general, it leads to a decrease of the PF, except in extreme cases of high compactness or extended accretion columns, where photons bent through large angles create a bright ring around the NS, resulting in a spike in the light curve (Pechenick et al. 1983; Falkner 2018; Markozov & Mushtukov 2024). For the geometry adopted here (see Fig. 14), our simulations show that although the overall shape of the PPs and PFS is preserved in the flat space-time projection, accounting for light bending changes the pulse fraction across the whole spectrum (on average by ∼40%) and makes the minima in PPs significantly shallower.

Although earlier works incorporated more complex physics – such as relativistic cross sections in Bulik et al. (1992) and self-consistent atmospheric structure calculations in Bulik et al. (1995) – the otherwise similar models to the one presented were generally unable to describe the data of 4U 1538–52 in a statistically satisfactory manner without introducing additional components of ambiguous physical interpretation. For example, even with relativistic cross-sections, the shape of the cyclotron line in Bulik et al. (1992) deviated significantly from the observed one, prompting the authors to adopt a mixed, arbitrary weighted contributions from emitting surface elements with different magnetic fields, covering a wide range of cyclotron energies, E_cyc = 16 − 26 keV to obtain the total flux. Bulik et al. (1995) also encountered the difficulty of simultaneously describing a deep shape of the cyclotron line and the phase modulation of the flux. They introduced an additional “nonmagnetic” component summed with the magnetic one. The latter determined the location of the cyclotron line, while the nonmagnetic component played role of a continuum contribution, reducing the line depth. Taking into account these facts, together with the high quality of the data presented in the current work and similar state of the physical modeling, we do not aim to perform phase-dependent spectral fitting, which would inevitably require the introduction of additional components. Instead, we argue that the general trends observed in the PPs, spectral flux, and PFS of the source can be qualitatively described and well understood within the framework of a physical model based on thermal resonant Comptonization in the slab-like emission regions and its relativistic projection onto the observer’splane.

5.3. Hints of a cyclotron harmonic in a single phase bin

Our phase-resolved spectral analysis of ObsID 304 revealed residuals around 40 keV in the phase bin of the primary peak A (0.375−0.500), shown in Fig. B.2 of the Appendix B. Modeling these residuals with an additional absorption feature using a Lorentzian profile, whose width was fixed to that of the first harmonic, as it could not be otherwise constrained, improves the reduced chi-squared from 243.9/234 to 217/232. The best-fit parameters indicate a line energy at 40.0^+1.3−1.1 keV with a depth of 0.74^+0.28−0.26, corresponding to a detection significance of approximately 2.6σ, which is suggestive but not yet statistically compelling.

Given that the energy of the absorption feature is almost exactly twice the fundamental cyclotron line energy, one possible origin could be the harmonic of the cyclotron line. This specific phase bin appears to be unique among the four NuSTAR observations: as discussed in Section 3, the first peak in observation 304 exhibits a more structured profile, showing a higher count rate at a phase of approximately 0.4. In contrast, the other three observations do not display this increase in count rate at the same phase.

We argue that the higher statistics at this phase might facilitate the detection of this potential harmonic line, whereas for the same phase bin of the other observations, the non-detection is mainly due to lack of statistics. A harmonic CRSF has been reported for this source before (Robba et al. 2001; Rodes-Roca et al. 2009), though at significantly higher energies, suggesting that the feature detected here may have a different origin or correspond to a different accretionregime.

Finally, it should be kept in mind that the formation of the second harmonic may not mirror that of the fundamental line. As shown by Schwarm (2017) (Fig. 5), the first and second harmonics reach their maximum depths at different angles relative to the magnetic field, which can be attributed to differences in redistribution and the ratio of E_cycⁿ to the plasma temperature. It is therefore possible that the second harmonic emerges more clearly at certain viewing angles or magnetic configurations, which could be met only in specific pulse phases–such as in this particular bin of ObsID 304.

6. Conclusions

We analyzed all available NuSTAR data of 4U 1538–52 to investigate the energy-resolved PP variations and the phase-dependent spectral variability, with a particular focus on the fundamental CRSF. Through decomposing the PPs into Fourier harmonics and computing the PFS, CC, and lag spectra, we identified key spectral features associated with the CRSF. Our timing and PP analysis revealed that the secondary peak of the PP disappears at the CRSF energies, which correspond to where a broad bump is observed in the PFS. Phase-resolved spectroscopy further shows that the CRSF is significantly deeper during the phase bins of the secondary peak, consistent with the disappearance of the secondary peak at theseenergies.

We demonstrated that broad bumps and dips observed in the PFS around the CRSF energy are a natural consequence of strong phase-dependent variations in the CRSF parameters. In the case of 4U 1538–52, the observed bump arises because the CRSF is much deeper in the low-flux phase bins, while the high-flux bins show a much shallower feature. Conversely, in sources like Her X-1, where dips are seen in the PFS, the CRSF is shallower during the off-pulse (low-flux) phases and deeper during the pulse peak. Through a simple exercise, we also showed that the energy of these features in the PFS does not necessarily coincide with the spectral energy of the CRSF when the line energy itself varies with spin phase. Specifically, if the CRSF energy increases with flux, as seen in both Her X-1 and 4U 1538–52, then the location of the dip or bump in the PFS will shift accordingly: if the line is deeper at high flux, a dip will appear at higher energies; if it is shallower, a bump will appear at lower energies. This naturally explains, for instance, why the dip in the PFS of Her X-1 is observed at a higher energy than the CRSF centroid, while the bump in the PFS of 4U 1538–52 appears at a lower energy. Furthermore, variations in the cutoff energy as a function of the spin-dependent flux can account for the overall upward or downward trends observed in the PFS. These findings underscore the importance of the PFS as a diagnostic tool for CRSFs, providing an independent and direct way to probe their spectral spin-dependent behavior. The PFS can offer immediate information on the relationship between the CRSF strength as a function of the spin-dependent flux, enabling a more direct interpretation of the underlying accretion geometry and magnetic fieldstructure.

To provide a qualitative description to the observed energy- and phase-dependent phenomena, we also performed a modeling of emission from strongly magnetized plasma and studied the effect of geometry – the magnetic pole locations and the observer’s line of sight – on PPs, spectra, and energy dependency of the PF. Our model indicates the high observer’s inclination and magnetic angle, with a difference between the two angles of |i − Θ₁|≲25° based on a significant reduction of the secondary peak in the PP. Within the frame of the model, the behavior of the CRSF at the two peaks associated with the two poles, is explained by observing photons emitted under different angles to the magnetic field at these phases. This results in a shallow line with a slightly higher centroid energy at the primary peak and a deeper line with a lower energy at the secondary peak. In turn, this leads to the discussed variation of the PFS at energies near the fundamental CRSF. The behavior of the CRSF described here favors line formation under thermal Comptonization. In contrast, the angular dependence of line parameters formed under the influence of bulk Comptonization is more complex and is associated with significantly larger shifts in centroid line energies (see, e.g., Schwarm 2017, their Chapter 5). We found a similar asymmetry in the location of the two poles as reported before, namely, ∼10 − 15°. Both in terms of previous studies (Bulik et al. 1992, 1995) and the offset from a pure dipole, indicate that the conditions at the two poles, that is, the characteristic magnetic field strength and the plasma temperature, might differ slightly. We would like to note that to describe the qualitative behavior observables, the two identical poles assumed in this work were sufficient. However, for the future studies aiming at the more careful and statistically acceptable description of the data, this and other effects will likely need to beconsidered.

Data availability

The data of the spectral and timing analysis are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/700/A70. An online service that reproduces our current results of timing analysis is available on the Renku-lab platform of the Swiss Science Data Centre at https://renkulab.io/projects/carlo.ferrigno/ppanda-light/

Acknowledgments

This research was supported by the International Space Science Institute (ISSI) in Bern, through the ISSI Working Group project Disentangling Pulse Profiles of (Accreting) Neutron Stars and by the European Space Agency (ESA) through the Archival Research Visitor Program. AD, EA, GC, VLP acknowledge funding from the Italian Space Agency, contract ASI/INAF n.I/004/11/4. CP, MDS, AD, FP, GR acknowledge support from SEAWIND grant funded by the European Union – Next Generation EU, Mission 4 Component 1 CUP C53D23001330006. FP, AD, MDS, GR and CP aknowledge the INAF GO/GTO grant number 1.05.23.05.12 for the project OBIWAN (Observing high B-fIeld Whispers from Accreting Neutron stars). EA aknowledges support from the INAF MINI-GRANTS number 1.05.23.04.04 for the project PPANDA (Pulse Profiles of Accreting Neutron Stars deeply analyzed). AA acknowledges financial support from ASI-INAF Accordo Attuativo HERMES Pathfinder operazioni n. 2022-25-HH.0 ESL acknowledges support from Deutsche Forschungsgemeinschaft grant WI 1860/11-2. We made use of Heasoft and NASA archives for the NuSTAR data. We developed our own timing code for the epoch folding, orbital correction, building of time-phase and energy-phase matrices. This code is based partly on available Python packages such as: astropy (Astropy Collaboration 2013), lmfit (Newville et al. 2024), matplotlib (Hunter 2007), emcee (Foreman-Mackey et al. 2013), stingray (Huppenkothen et al. 2019), corner (Foreman-Mackey 2016), scipy (Virtanen et al. 2020).

References

Appendix A: Summary plots of the pulse profile properties of 302, 306-2 and 306-4

Fig. A.1. PP properties for 4U 1538-522 in ObsID 302 (see the caption of Fig.5), for 32 phase bins and a S/N of 15.

Fig. A.2. PP properties for 4U 1538-522 in ObsID 306 (see caption of Fig.5), for 32 phase bins and a S/N of 15.

Fig. A.3. PP properties for 4U 1538-522 in ObsID 306-4 (see caption of Fig.5), for 32 phase bins and a S/N of 15.

Appendix B: Spectral analysis plots

Fig. B.1. Observation 306-2 PA

B.1. Phase-resolved spectral analysis of 304

Fig. B.2. Observation 304: A grid of eight figures in four rows and two columns, showing phase bins 0.000-0.125 to 0.875-1.000.

B.2. Phase-resolved spectral analysis of 306-2

Fig. B.3. Observation 306: A grid of eight figures in four rows and two columns, showing phase bins 0.000-0.125 to 0.875-1.000.

Fig. B.4. Phase-resolved parameters

Appendix B: Tables of spectral analysis

Appendix D: Good time intervals

Fig. D.1. Filtering of light curves using a log-normal distribution of counts in the 3–70 keV band. Left panels: Total light curve of each observation. The horizontal lines mark the adopted count rate limits, while the pink vertical bands indicate time intervals that were excluded from further analysis. Right panels: Logarithmic histograms of count rates, with horizontal lines denoting the same adopted limits, where good time intervals were defined as those corresponding to count rates within the ±5σ of the central peak of the log-normal distribution.

All Tables

All Figures

	Fig. 1. Energy-phase map of the observation 304, generated with 32 phase bins and a minimum S/N of 15: normalized PPs at equally spaced logarithmic energy intervals. Each pulse in the bin was normalized by subtracting its mean and dividing by the standard deviation.
In the text

Fig. 2. Alignment of the energy-phase matrices. The plot shows the SSD values calculated for each phase shift between the PP of a given observation and the reference PP (ObsID 304). The minimum SSD value corresponds to the best alignment, and the points around this minimum are fitted with a Gaussian to refine the phase shift determination and account for statistical fluctuations. The centroid of the Gaussian provides the precise phase shift used for alignment.

In the text

	Fig. 3. Aligned PPs in the 3–70 keV band (errors not shown for clarity’s sake). Two periods are shown for better clarity.
In the text

	Fig. 4. Energy-selected PPs for all examined observations.
In the text

Fig. 5. Summary plot for ObsID 304 using 32 phase bins and a S/N 15. Panel a: PF (green points) along with its best-fit model (solid lines); polynomial fits are also included for comparison. Panel b: residuals from the fit. Panels c–f: Phases and normalized amplitudes of the first (A1, ϕ1) and second harmonic (A2, ϕ2), respectively. The vertical shaded regions denote the energy ranges and widths of the Gaussian functions fitted to the PF. Panel g: Selection of normalized PPs at equally spaced logarithmic energy intervals, horizontally offset for clarity. Each pulse in the bin was normalized by subtracting its mean and dividing by the standard deviation. Panel h: Color-map display of the normalized PPs as a function of energy, with 20 evenly spaced contour lines drawn for reference. Panel i: CC between the PP at each energy band and the overall average profile. Panel j: corresponding phase lags. The colored vertical regions are consistent with those in panels d–f.

In the text

	Fig. 6. Total 3–70 keV PP of ObsID 304 and the first two Fourier harmonics from its FFT decomposition.
In the text

	Fig. 7. CC spectrum extracted during ObsID 304 using the PP from the 35–70 keV energy range as template.
In the text

	Fig. 8. (a) Data and best-fitting model for the phase-averaged spectrum of ObsID 304. (b) Residuals of the best-fitting model when the depth of the cyclotron line is set to zero. (c) Residuals of the best-fitting model.
In the text

	Fig. 9. Definition and naming of the phase bins used in the analysis.
In the text

Fig. 10. Phase-dependent spectral variability for ObsID 304: each curve represents one of the best-fitting model from the eight phase-selected spectra. Models are normalized to unity at the reference energy of 10 keV. The blue and red curves correspond to spectra extracted in the phase intervals of the primary and secondary peak, respectively. The black spectrum shows a very shallow cyclotron line in the primary falling flank. The gray spectra correspond to the primary rising flank and secondary falling flank.

In the text

	Fig. 11. CRSF energy versus line depth values for the phase bins corresponding to the primary and secondary peaks (0.250–0.500 and 0.625–0.875).
In the text

	Fig. 12. PFS of observation 304, calculated using 8 phase bins and an S/N of 5 for consistency with the phase-resolved spectral analysis, overplotted with the primary (blue) and secondary (red) peak spectra.
In the text

Fig. 13. (a) Simulated spectra for the primary and secondary peaks, modeled by a power-law continuum with an exponential energy cutoff. The cutoff energy is set at 12 keV for the primary peak spectrum (blue) and 8 keV for the secondary peak spectrum (red). The absorption line in the secondary peak spectrum is centered at 20 keV with a width of 3 keV, while the absorption line in the primary peak spectrum is at 22 keV with a width of 4 keV. (b) Simulated PFS of the spectra shown in (a) using the CV proxy. The red and blue vertical lines highlight the energies of the absorption feature in the secondary and primary peaks, at 20 and 22 keV, respectively. (c) Simulated spectra for the primary and secondary peaks, where the energies of the absorption features have been exchanged between the two. (d) Simulated PFS of the spectra shown in (c) using the CV proxy.

In the text

	Fig. 14. Geometry adopted for 4U 1538–52, showing the NS, magnetic poles at Θ1 = 80° and Θ2 = 105°, azimuthal separation ΔΦ = 168°, and inclination to the observer’s line of sight i = 67°.
In the text

Fig. 15. Phase-dependent emission angle and various observables modeled for a selected geometry: i = 67°, Θ1 = 80°, Θ2 = 105° and ΔΦ = 168°. The mass and radius of the NS are MNS = 1.4 M⊙ and RNS = 12 km, the size of spot-like poles is r0 = 300 m. The emission from each emitting surface element is calculated for: Ecyc = 25 keV, kTe = 6 keV, ne = 1023 cm−3 and τT = 20. (a) Emission angle, visible for a remote observer at each rotation phase. The angle is calculated for a center of a pole on the NS surface in the Schwarzschild metric (blue dashed-dotted line corresponds to pole 1, red dashed line – to pole 2). The shadowed region shows when a pole is obscured by the NS. (b) Energy-integrated (5 − 80 keV) pulsed profile (black thick line) for the same geometry, with individual pole contributions. Light-gray line shows the PP in the energy range 19 − 24 keV, corresponding to the core of the fundamental CRSF in the observer’s rest frame. (c) Modeled energy flux near the fundamental line. The blue and red lines are obtained for the central regions of the primary (“peak 1”) and secondary (“peak 2”) peaks in the PPs; the respective ranges of phases are marked by shadow regions in panel b. The black dotted line marks the redshifted value for the cyclotron energy. (d) PFS in the same energy range as in panel c.

In the text

	Fig. A.1. PP properties for 4U 1538-522 in ObsID 302 (see the caption of Fig.5), for 32 phase bins and a S/N of 15.
In the text

	Fig. A.2. PP properties for 4U 1538-522 in ObsID 306 (see caption of Fig.5), for 32 phase bins and a S/N of 15.
In the text

	Fig. A.3. PP properties for 4U 1538-522 in ObsID 306-4 (see caption of Fig.5), for 32 phase bins and a S/N of 15.
In the text

	Fig. B.1. Observation 306-2 PA
In the text

	Fig. B.2. Observation 304: A grid of eight figures in four rows and two columns, showing phase bins 0.000-0.125 to 0.875-1.000.
In the text

	Fig. B.3. Observation 306: A grid of eight figures in four rows and two columns, showing phase bins 0.000-0.125 to 0.875-1.000.
In the text

	Fig. B.4. Phase-resolved parameters
In the text

Fig. D.1. Filtering of light curves using a log-normal distribution of counts in the 3–70 keV band. Left panels: Total light curve of each observation. The horizontal lines mark the adopted count rate limits, while the pink vertical bands indicate time intervals that were excluded from further analysis. Right panels: Logarithmic histograms of count rates, with horizontal lines denoting the same adopted limits, where good time intervals were defined as those corresponding to count rates within the ±5σ of the central peak of the log-normal distribution.

In the text