© The Authors

1. Introduction

Although most Be X-ray binaries (BeXRBs) are transient sources composed of a neutron star (NS) in a wide and eccentric orbit (with e>0.3) around the Be star, a growing class of BeXRBs is that of the persistent Be/NS binaries (Reig & Roche 1999), which are characterized by a rather constant and low luminosity (L_X∼10³⁴⁻³⁵ erg s⁻¹) and a long spin period (P_spin>100 s). These properties suggest that the NS orbit around the Be star is rather wide (P_orb≳30 days) and nearly circular (e<0.2, Pfahl et al. 2002). Therefore, it is continuously accreting material from the low–density regions of the fast polar wind of the companion star.

The first sample of these sources was composed of the four sources X Persei (which is considered the prototypical persistent BeXRB), RX J0146.9+6121, RX J1037.5-5647, and RX J0440.9+4431. In recent years, the observations performed with the most sensitive X-ray observatories led to the identification of other confirmed or candidate sources of this type, such as Swift J045106.8-694803 (Bartlett et al. 2013), CXOU J225355.1+624336 (La Palombara et al. 2021), XTE J1906+090 (Sguera et al. 2023), and 4XMM J182531.5-144036 (Mason et al. 2024). This suggests that this class of X-Persei–like persistent BeXRBs could constitute a significant part of the population of Galactic unidentified low-luminosity X-ray sources.

An additional candidate member of the persistent BeXRBs is the Galactic Be pulsar 4U 0728-25, which is a poorly studied X-ray source. It was discovered with Uhuru (Forman et al. 1978) and observed also with Ariel V (Warwick et al. 1981) and HEAO-1 (Wood et al. 1984). The optical counterpart, V441 Pup, is at an estimated distance (based on Gaia EDR3 data) $d=7.6_{-0.8}^{+1.0}$ kpc (Bailer-Jones et al. 2021). According to Negueruela et al. (1996), it is classified as an O8–9Ve star, while for Maíz Apellániz et al. (2016) it is a star of O5Ve type; very likely, the reason for this disagreement between the two spectral classifications is that only Negueruela et al. (1996) corrected for the He infilling emission. Subsequent RXTE observations led to the discovery of pulsed emission (P_spin = 103.2±0.1 s) and to the detection of the orbital period (P_orb≃ 34.5 days; Corbet & Peele 1997). The source monitoring performed between 2005 and 2016 with the BAT instrument of the Swift Neil Gehrels Observatory confirmed the orbital period and revealed an almost steady emission with no brightening (Corbet et al. 2016). The source was also included in the catalogues of persistent X-ray sources detected with BeppoSAX WFC (Capitanio et al. 2011) and INTEGRAL IBIS (Bird et al. 2016). The latest observation of 4U 0728-25 was performed in 2016 with AstroSat, when the source was in a non-flaring persistent state at a flux level f_X (0.4–7 keV) = 6.6×10⁻¹¹ erg cm⁻² s⁻¹. It provided a new measure of the pulse period (P_spin = 103.144±0.001 s) and revealed an energy-dependent pulse profile and a broad (with width σ = 1.1 keV) iron line at E≃6.3 keV in the energy spectrum (Roy et al. 2020).

The almost constant X-ray flux of 4U 0728-25 over long timescales, and its long orbital and spin periods, strongly suggest that also this source belongs to the class of the X-Persei–like persistent BeXRBs. On the other hand, the observed X-ray flux and estimated distance imply a rather high luminosity (L_X∼ a few 10³⁵ erg s⁻¹) compared with the other sources of this type. Therefore, we observed it with XMM–Newton to investigate its timing and spectral properties.

2. Observation and data reduction

The XMM–Newton observation of 4U 0728-25 was performed on April 11, 2023 (MJD 60045) and its duration was ≃33.5 ks. The source was observed with both the pn (Strüder et al. 2001) and the two MOS (Turner et al. 2001) cameras of the EPIC instrument; in addition, also the Reflection Grating Spectrometer (RGS) was used (den Herder et al. 2001). Due to the rather high source flux, which implies a high count rate (CR), for the EPIC cameras the Small Window operating mode was used, in order to prevent photon pileup. The configuration of the observing instruments is reported in Table 1, where the net exposure time of the EPIC cameras does not include the readout dead time, which contributes for 2.9% and 29.9% for the MOS and pn cameras, respectively.

We used the version 21 of the Science Analysis System (SAS) to process the collected data. First of all, we looked for the possible presence of periods characterized by a high rate of soft protons (SP), which for EPIC could affect the results of the data analysis. We verified that no SP flares occured during the observation: the rate of SPs was negligible along the whole exposure and, then, we used the whole exposure time.

In the case of the EPIC data, we considered events with pattern in the range 0–12 (corresponding to up to 4 contiguous pixels) for the MOS data and in the range 0–4 (i.e. mono- and bi-pixels) for the pn data. In all EPIC cameras the source position was very near the window edge. Therefore, we were forced to select source events only from a circular region, centered at the source position, with a small radius of 30 arcsec. In the case of the pn camera we had the same constraint also for the selection of the background events, for which we considered a circular extraction region of the same size; instead for both the MOS cameras we considered a background extraction radius of 120 arcsec. Finally, we verified also that the data were not affected by photon pileup: in fact, although the count rate was rather high (Table 1), it was well below the threshold for significant pile-up estimated by Jethwa et al. (2015).

Our analysis revealed that the RGS data do not show the presence of any significant emission or absorption feature and provided no further information compared with the EPIC ones. Therefore, in the following we will not discuss them.

3. Timing analysis

The arrival times of the EPIC events were reported to the Solar System barycentre by using the SAS tool BARYCENTER. In our timing analysis, we considered both the full energy range 0.15–12 keV and two sub-ranges, one soft (S) between 0.15 and 2.5 keV and one hard (H) between 2.5 and 12 keV: in this way, we obtained a similar source count number in both ranges. For each of the three cameras, we accumulated a light curve in each of these three energy ranges; then, we used the SAS tool EPICLCCORR to properly take into account both the background and the extraction region. In this way we obtained a net source CR in the full energy range of ≃2.8 counts s⁻¹ and ≃1.0 counts s⁻¹ for, respectively, the pn camera and each of the two MOS cameras. Afterwards, for each energy range we summed the light curves of the three cameras to obtain the overall EPIC light curves. These curves are shown in Fig. 1, toghether with the hardness ratio (HR) between the H and the S light curves (HR = H/S). The average CR is ≃2.5 and ≃2.3 counts s⁻¹ for the S and the H range, respectively.

Fig. 1. Light curves of 4U 0728-25, corrected for the background and the extraction region, in the S (top panel) and H (second panel) energy ranges, together with their ratio (third panel); the binning time is 100 s/bin. In the bottom panel the curve in the full energy range is reported, with a binning time of 50 s/bin.

The light curves reported in Fig. 1 clearly show that 4U 0728-25 is characterized by a very large flux variability on short time scales: in both the partial ranges, where the bin time is 100 s, the CR varies by a factor up to ∼3 between two consecutive bins. This is confirmed by the large values of the CR standard deviation reported in Table 1. This finding is further confirmed in the full energy range, where the binning time is 50 s only: also in this case, indeed, we can observe a comparable variability factor between two nearby bins. On the other hand, the source is rather stable on the long time scales, since all the large short-scale variations occur always around the same average CR value. Even the HR is rather constant: it shows variations of ∼50% at most and it is not correlated with the source CR.

Once corrected for the Solar System barycentre, the datasets of the three instruments were merged into a unique event list to increase the count statistics. The best-fit period, obtained with the phase-fitting technique, is P = 103.301(5) s. Based on this period, we built the folded light curve in each of the three energy ranges, together with the folded HR between the H and the S curves. They are reported in Fig. 2. In both energy ranges, the pulse profile is characterized by a single broad peak, which extends for almost half of the pulse period. Both the rise and the decay of the CR are rather smooth, but in the H range the peak is slightly delayed compared with the S range. The HR indicates some spectral variation during the rising part of the pulse.

Fig. 2. Light curves of 4U 0728-25 in the full (bottom panel), S (top panel) and H (second panel) energy ranges, together with the HR (third panel), folded at the pulse period P = 103.301 s. The phase ranges selected for the phase-resolved spectral analysis are delimited with vertical dashed lines.

The dependence of the pulse profile on the photon energy is confirmed by Fig. 3, where we report the folded light curves in four narrower energy ranges: 0.15–1.5, 1.5–2.5, 2.5–4, and 4–12 keV. As in the previous case, all of them show a single broad peak, with no specific features. The flux variability shows a clear dependence with the photon energy: the average pulsed fraction (PF), defined as PF = (CR_max–CR_min)/(2 × CR_average), gradually increases from ≃15% at E < 1.5 keV to ≃23% at E > 4 keV.

Fig. 3. Pulse profile of 4U 0728-25 in the energy ranges 0.15–1.5, 1.5–2.5, 2.5–4, and 4–12 keV.

4. EPIC time-averaged spectral analysis

Although the source flux of 4U 0728-25 was highly variable on short timescales during the XMM–Newton observation, over the whole observation it showed neither an increasing nor a decreasing trend. In addition, the HR showed a limited variability (Fig. 1), thus indicating that the spectral properties were almost constant. For these reasons we performed the spectral analysis of the source over the whole exposure, with no time selection. To this aim, we obtained a unique source spectrum for each EPIC camera by selecting events from the same extraction regions already considered for the light curves. We used the SAS tool SPECGROUP to perform the rebinning of each spectrum with a minimum significance of 3σ for each energy bin. Then, the spectral channels consistent with zero (after the background subtraction) were removed. The response matrices and ancillary files were produced using, respectively, the SAS tasks RMFGEN and ARFGEN. The spectral analysis was performed using version 12.12.1 of XSPEC (Arnaud 1996); to this aim we considered the energy range between 0.2 and 12 keV. In all cases the spectral uncertainties were calculated at the 90% confidence level for each interesting parameter and we considered a source distance d = 7.6 kpc. The spectral fitting was performed on the three EPIC spectra simultaneously, since there was no evidence that the separate fits provided inconsistent results. This was done by introducing a free normalization parameter among the three cameras, in order to properly take into account the possible uncertainties in the instrumental responses: more precisely, once the normalisation factor of the pn spectrum was fixed to 1, we obtained relative normalisation factors 0.938 ± 0.008 and 1.091 ± 0.008 for, respectively, the MOS1 and the MOS2 spectra. In our analysis we used the elemental abundances reported by Wilms et al. (2000), the photoelectric absorption cross-sections provided by Verner et al. (1996), and the absorption model TBABS defined in XSPEC. The estimated fluxes were calculated with the XSPEC tool CFLUX.

Our first attempt was to describe the EPIC continuum spectrum with a simple power-law (PL) model that provides a phenomenological description of the emission due to the matter accretion from the companion star onto the neutron star: with it the best-fit values of the photon index and interstellar absorption were, respectively, $\mathrm{\Gamma }=1.{33}_{-0.1}^{+0.2}$ and $N_{\mathrm{H}}=(6.9_{-0.2}^{+0.1})\times {10}^{21}$ cm⁻², and the corresponding reduced chi-squared was ${\chi }_{\nu }^2$ = 1.215 for 2508 degrees of freedom (d.o.f.). Starting from this model, we obtained a significantly better result in two alternative ways: either with the introduction of a partial covering fraction absorption (TBPCF in XSPEC), which describes the absorption due to an inhomogeneous absorber medium, or by replacing the PL model with a cutoff power law (CPL). In both cases the reduced chi-squared decreased to ${\chi }_{\nu }^2$ = 1.063, with 2506 and 2507 degrees of freedom (d.o.f.), respectively, and the observed source flux is f_X≃1.7×10⁻¹¹ erg cm⁻² s⁻¹. We tested also models with a single thermal component, but in all cases we obtained a worse fit.

We obtained a slightly better fit (${\chi }_{\nu }^2$/d.o.f. = 1.058/2505) with a CPL plus a component due to collisionally ionised gas (APEC model in XSPEC) that accounts for the thermal emission from a possible optically thin gas around the accreting NS. An even better fit (${\chi }_{\nu }^2$/d.o.f. = 1.050/2506) is obtained with a PL plus a black body (BB), where the latter component describes the thermal emission of a possible optically thick plasma on the NS surface, a scenario that is often observed in this type of sources. In both cases the non-thermal component (either the PL or the CPL) dominates over the thermal one (either the BB or the APEC). The F-test statistics value is very high (≃200) for the PL + BB model, which implies a very negligible probability (<10⁻⁸⁰) that the additional BB component is an artifact. This is confirmed by the spectrum reported in Fig. 4, where we show that the data-model residuals are much larger in the case of the fit with a PL or a BB model, compared with the PL + BB model. For the CPL + APEC model, instead, the F-test statistics value is much lower (≃7.7) and, then, the probability that the APEC component is spurious is just below 10⁻³. Apart from these two models, with other two-component models it was not possible to obtain equally good fits. In particular, the PL + APEC model resulted in a significantly worse fit (${\chi }_{\nu }^2$/d.o.f. = 1.182/2506), while the best-fit CPL + BB model was fully equivalent to the PL + BB model, since E_cut>100 keV.

Fig. 4. Time-averaged spectrum of 4U 0728-25, where the black, red, and green symbols represent the pn, MOS1, and MOS2 data, respectively. Upper panel: superposition of the best-fitting PL + BB model, for the pn spectrum only, with the EPIC spectra. Middle panels: residuals (in units of σ) between data and model in the case of the single PL and of the single BB. Bottom panel: residuals in the case of the PL + BB.

We noted that the fit of the continuum spectrum with each of the previous models left weak residuals at some specific energies, such as at ≃0.55 keV. We verified that it was possible to describe these features with a Gaussian model, but further investigations of these features with Monte Carlo simulations, performed using the XSPEC routine SIMFTES, estimated a high probability that these features are not significant. Therefore, in the following analyses we will ignore them. Moreover, we looked for the possible presence of emission lines due to iron K_α transition, at energies between 6 and 7 keV and with widths in the range 0–0.5 keV (in steps of 0.1 keV). We did not find any feature of this type, with an upper limit of ≲100 eV (at 90% c.l.) on its equivalent width.

The best-fitting parameters we obtained by using the four previous models are reported in Table 2. All models are characterozed by a moderate interstellar absorption, since N_H ∼(5−7)×10²¹ cm⁻², and a rather hard non-thermal component, since Γ≃0.7 in the case of the CPL and 1.5–1.8 in the case of the PL; moreover, the cut-off energy of the CPL model is rather low. The unabsorbed flux of the non-thermal component is always f_X∼(2−3)×10⁻¹¹ erg cm⁻² s⁻¹ and, in the cases of the two-component models, contributes most of the source flux and dominates over the thermal component at all energies. Both the additional thermal components have a rather high temperature but, while in the case of the BB it provides a significant fraction of the total source flux, in the case of the APEC its contribution is only marginal.

5. EPIC phase-resolved spectral analysis

Fig. 2 reveals that the HR is variable along the pulse phase, with a primary and a secondary peak just after and before the flux maximum, respectively. Since this implies that the source spectrum varies depending on the spin phase, we performed a specific spectral analysis in two distinct phase ranges: for each EPIC camera, we accumulated a hard spectrum (H) in the phase range Δϕ = 0.675−0.925 (where HR > 0.95) and a soft spectrum (S) in the phase range Δϕ = 0.975−1.275 (where HR < 0.9). In this way, the two spectra are characterized by very similar number of source counts, in spite of the different length of the two phase intervals (0.25 and 0.3 for, respectively, spectrum H and S).

Our first aim was to verify if the best-fitting models of the time-averaged spectrum can provide a satisfactory description also of the two phase-resolved spectra. Therefore, we performed an independent fit of the two spectra using the four models described in Section 4. The results obtained with this approach are reported in Table 3. We found that both spectra can be described with each model and that in all cases the source flux decreases from (1.88 ± 0.05) × 10⁻¹¹ erg cm⁻² s⁻¹ for spectrum H to (1.34 ± 0.04) × 10⁻¹¹ erg cm⁻² s⁻¹ for spectrum S, corresponding to a reduction of ≃28%. However, for each model there is a significant difference of the best-fit parameter values between the two spectra: this holds for the interstellar absorption, the photon index of the non-thermal component and, in the case of the PL + BB model, also for the temperature of the thermal component. In the case of the PL + BB model the relative contribution of the two components to the total flux does not vary. For the CPL + APEC model, instead, it is not possible to make any comparison, due to the large uncertainties in the parameters of the APEC component.

The previous results show that the spectral shape varies along the pulse phase. In the following, we try to investigate whether, in the PL + BB model, the variation can be ascribed to only one of the two components. In this case, we found that it was not possible to describe the two spectra by allowing a variation of only the PL flux (Null hypothesis probability ≃0.02), while we obtained a good fit with a common PL component and variable BB flux (Null hypothesis probability ≃0.22). We also investigated two alternative scenarios regarding the variability of the two components:

• On the one hand, we forced common values between the two spectra for N_H, Γ, and kT_BB, allowing only the normalization of the two components to vary (Case 1): this implies that the spectral variablity along the pulse phase is attributed to a difference in the relative contribution of the two components.

• On the other hand, we forced a common BB component by allowing only the PL parameters to vary (Case 2): in this case the spectral variability is fully attributed to a change in the non-thermal component.

The results obtained with this approach are reported in Table 4. They show that the two cases are fully comparable from the point of view of the fit quality, since they have almost equal ${\chi }_{\nu }^2$: in Case 1 the flux variability is almost fully due to a variation of the BB flux, since the PL flux is comparable in the two spectra; in Case 2 the constant BB component needs that not only the flux but also the photon index of the PL component is significantly variable between the two spectra.

6. Discussion

After its discovery in the ’70s, 4U 0728-25 was detected several times by various X-ray telescopes. In the upper panel of Fig. 5 we provide an overview of the X-ray luminosity levels (in the energy range 2–10 keV), which are estimated based on the flux detected by each telescope and assuming a source distance of 7.6 kpc. The first three values are based on the flux value reported in overall catalogues of X-ray sources detected with Uhuru (Forman et al. 1978) and Ariel (Warwick et al. 1981; Steiner et al. 1984): since these catalogues were obtained with observations performed over a long time span, the reported luminosity values should be considered as average values (as indicated by the horizontal error bars). For these observations, and also for the detections obtained with HEAO-1 (Steiner et al. 1984; Wood et al. 1984), we estimated the 2–10 keV flux based on the published CR values and the appropriate count-rate-to-flux conversion factors, assuming a Crab-like PL spectrum; instead, in the case of the RXTE and AstroSat observations, we used the reported best-fit models to estimate the source flux. The figure shows that the first four detections revealed 4U 0728-25 at a luminosity level in the range (5–6) × 10³⁵ erg s⁻¹, while in the observations 5–8 the luminosity level remained always in the range (1–4) × 10³⁵ erg s⁻¹; moreover, the last observation performed with XMM–Newton in 2023 revealed 4U 0728-25 at the lowest flux ever, with L_X≃9×10³⁴ erg s⁻¹. The reported values show that, over a time interval of about 50 years, 4U 0728-25 was characterized by a remarkable stability, since its flux varied less than a factor 10; the corresponding luminosities were always in the range ∼10³⁵⁻³⁶ erg s⁻¹.

Fig. 5. Long-term evolution of the luminosity in the energy range 0.2–12 keV (upper panel) and of the pulse period (lower panel) of 4U 0728-25. Telescope and epoch: (1) Uhuru catalogue: 12/12/1970–18/03/1973 (Forman et al. 1978); (2) Ariel V: 01/11/1974–30/04/1977 (Steiner et al. 1984); (3) Ariel V 3A catalogue: 15/10/1974–14/03/1980 (Warwick et al. 1981); (4) HEAO1: 19–23/10/1977 (Steiner et al. 1984); (5) HEAO1: 15/08/1977–15/02/1978 (Wood et al. 1984); (6) RXTE: 05/05/1997 (Corbet & Peele 1997); (7) AstroSat: 06–07/05/2016 (Roy et al. 2020); (8) XMM–Newton: 11/04/2023 (this work).

The pulse period of 103.301 ± 0.005 s measured with XMM–Newton in 2023 is ≃0.15% longer than that of 103.144 ± 0.001 s measured with AstroSat in 2016 (Fig 5, lower panel) and cannot be attributed to the NS orbit around the companion Be star. In fact, assuming a mass M_Be∼10 M_⊙, from the orbital period P_orb≃34.5 days we can estimate an average orbital velocity v_orb≃140 km s⁻¹, which implies a variation of the spin period ΔP_spin=P_spin×v_orb/c≃0.05 s. Since this value is about three times lower than the dP_spin between the XMM–Newton and AstroSat observations, the observed P_spin increase is most likely due to the NS spin down, which occured at an average rate Ṗ = 7.2 × 10⁻¹⁰ s s⁻¹ during the seven years which separates the two observations. For comparison with other persistent BeXRBs, we note that this value is a factor ∼10 lower than that estimated for RX J0440.9+4431 during its quiescence phase (La Palombara et al. 2012) and very similar to that measured for CXOU J225355.1+624336 (La Palombara et al. 2021). The increase in the spin period suggests that the matter transfer from the Be star onto the NS occurs via wind accretion, without any formation of an accretion disk.

The pulse shape of 4U 0728-25 observed with XMM–Newton is characterized by a single broad peak at all energies, which extends for almost half of the pulse period. It is very similar to that revealed with RXTE and AstroSat, which also observed a single-peak profile; moreover, the same type of pulse shape is typical of several other persistent BeXRBs, as in the case of RX J0146.9+6121 (La Palombara & Mereghetti 2006) and RX J1037.5-5647 (La Palombara et al. 2009). The HR is almost correlated (with a little delay) with the CR and varies of ≃±10% around its average value along the pulse phase: this finding is rather common in the binary pulsars and reveals a limited energy dependence of the pulse profile. We also observed a slight increase in the pulsed fraction with the photon energy but, compared with other persistent BeXRBs (La Palombara et al. 2024), in the case of 4U 0728-25 the PF is much lower.

Since a simple absorbed PL model provided a rather poor fit of the source spectrum, we used various more elaborated models to describe it: on one hand, we considered two models with a single non-thermal component, either a partially covered PL or a CPL; on the other hand, we introduced a two-component model composed of a simple PL or CPL plus a thermal component, either a BB or an APEC model. All these models provide a good description of the source spectrum and none of them is clearly favoured compared with the other ones. For all models the estimated hydrogen column density is N_H ∼(5−7)×10²¹ cm⁻², which is consistent with the values of interstellar absorption estimated with RXTE and AstroSat. For comparison, we note that the color excess of the optical counterpart is E(B–V) = 0.75 (Negueruela et al. 1996), which implies an extinction A_V = 3.1×E(B−V) = 2.325; since the average relation between optical extinction and X–ray absorption is N_H=(2.87±0.12)×10⁻²² A_V cm⁻² (Foight et al. 2016), this implieas that N_H=(6.7±0.3)×10²¹ cm⁻², a value fully consistent with our result.

In the two spectral models with a PL component, the photon index Γ is between 1.5 and 1.8, in agreement with the index values obtained in the cases of the RXTE and AstroSat spectra. Instead, with the two models with the CPL component we obtained best-fit values for Γ and E_cut that are significantly lower than for the corresponding parameters obtained with AstroSat. These differences are very likely due to the limited energy range covered by XMM–Newton (0.2–12 keV). Therefore, it is possible that the power-law and cutoff components are not well constrained, particularly in the two-component models. Finally, regarding the partially absorbed PL model, we note that the inhomogeneous absorber medium described with the TBPCF component is very likely the clumpy polar wind coming from the Be companion star, which crosses the line of sight to the X-ray source (Grinberg et al. 2020).

Compared with a simple PL model, in the PL + BB model the additional thermal component improves significantly the spectral fit; it contributes for ≃24% of the total X-ray flux and its temperature and radius are kT_BB≃1.5 keV and R_BB≃240 m, respectively. Therefore, for 4U 0728-25 the properties of the detected BB excess are very similar to those observed in several other persistent BeXRBs: for most of these sources, in fact, this component has a temperature in the range 1–2 keV and a size of few hundred meters, and contributes for ∼20–40% to the total X-ray flux (La Palombara et al. 2024). In the CPL + APEC model, instead, the additional thermal component provides only a marginal contribution to the total flux; moreover, the fit improvement obtained with the addition of the APEC component to the simple CPL is much lower than in the case of the BB addition to a simple PL.

In the EPIC spectrum of 4U 0728-25 we found no evidence of any emission or absorption feature, so we cannot confirm the presence of the broad Fe K_α line detected in the AstroSat spectrum. On the contrary, our analysis provided an upper limit of ≲100 eV on the equivalent width of a possible Fe line. It is very interesting to note that this result is fully aligned with what obtained for all the currently known persistent BeXRBs, since in none of them this type of feature has been observed up to now, and a comparable upper limit on its equivalent width has been estimated (La Palombara et al. 2024).

We performed a phase-resolved spectral analysis in order to investigate the spectral dependence on the pulse phase. To this aim, we compared the spectra accumulated in two different phase intervals, corresponding to the maximum (H) and the minimum (S) of the HR. The independent fit of these two spectra showed that the source flux decreases of ≃30% from phase H to phase S and that there is also a clear variability of the spectral parameters. For the PL + BB model, we found that the flux of both components is clearly variable and that the BB is significant at 99% c.l. in both spectra. Instead this is not true for the CPL + APEC model, since in both spectra the additional APEC component cannot be constrained. This result, together with the marginal APEC contribution to the time-average flux, disfavours the real presence of this type of component in the source spectrum. We also performed a simultaneous fit of the two phase-resolved spectra, in order to investigate in deeper detail the spectral variability. For the PL + BB model we verified that the spectral variability can be well described with a variation of only the BB flux, while this is not true in the case of a variation of only the PL flux; moreover, we verified that the two spectra are compatible with a constant BB component only if both the PL flux and photon index vary. Therefore, based on this results, the presence of a variable BB component is very likely.

The size of the BB emission region and the strong clues that it is variable suggest that it originates on the surface of the accreting NS, in agreement with the scenario proposed by Hickox et al. (2004) to explain the thermal excess in the low-luminosity HMXRBs. This hypothesis has been confirmed in most of the persistent BeXRBs, where it has been proven that the radius of the observed BB excess is consistent with the estimated radius of the accretion column onto the NS polar caps (La Palombara et al. 2024). In order to verify if this hypothesis holds also in the case of 4U 0728-25, we calculated the column radius R_col and compared it with R_BB. To this aim, for the accreting NS we assumed a mass M_NS = 1.4 M_⊙, a radius R_NS = 10⁶ cm, and a magnetic field B_NS = 10¹² G. Based on the estimated luminosity L_X(0.2–12 keV) ≃1.3×10³⁵ erg s⁻¹ and on the relation Ṁ = L_XR_NS/(GM_NS), we derived an accretion rate Ṁ ≃ 7 × 10¹⁴ g s⁻¹ and a magnetospheric radius R_m≃8.6×10⁸ cm (Campana et al. 1998). Finally, we used the relation R_col∼R_NS (R_NS/R_m)^0.5 (Hickox et al. 2004) to obtain R_col≃340 m. This value is only ∼40% larger than the estimated BB radius R_BB = 240 m, which implies that the BB size is comparable with that of the NS polar caps. Therefore, it is very likely that they are the source of the observed BB component, and this hypothesis is further strengthened by the BB variability along the pulse phase.

7. Conclusions

We have described the results obtained by analizing a ≃33 ks XMM–Newton observation of the BeXRB 4U 0728-25 performed in April 2023. The source was detected at a luminosity level L_X≃9×10³⁴ erg s⁻¹ in the energy range 2–10 keV: this is the lowest level ever observed since the source discovery in the ‘70 s, and it is also about 50% lower than the minimum luminosity previously observed. On the other hand, since this value is less than one order magnitude lower than the maximum luminosity observed in almost 50 years, this XMM–Newton observation confirmed the limited flux variability of 4U 0728-25 over long time scales.

We measured a pulse period P_spin = 103.301±0.005 s, which implies a spin down Ṗ = 7.2 × 10⁻¹⁰ s s⁻¹ as compared to the previous measurement in 2016. This moderate spin down is coherent with the low rate of matter transfer, via wind accretion, from the Be companion star onto the NS. The pulse shape shows a simple single-peak profile, with a moderate pulsed fraction at all energies.

The spectral analysis revealed the presence of a flux excess above the main power-law component, which can be described with a BB model with high temperature (kT_BB≃1.5 keV) and small emission radius (R_BB≃240 m); it contributes for ∼25% to the total flux and has a size comparable with that of the NS polar caps. Instead we found no evidence of emission features due to narrow (σ<0.5 keV) Fe K_α lines, with an upper limit of ≃100 eV on the equivalent width of this type of feature.

We also performed a phase-resolved spectral analysis, which proved the spectral variability along the pulse phase. In this way we proved that the BB component is always present, and we also obtained strong clues that it is variable along the pulse phase. Therefore, we attributed its origin to the NS polar caps.

In summary, we verified that the timing and spectral properties of 4U 0728-25 are fully consistent with those of most of the persistent BeXRBs. For this reason, we can consider it as a confirmed member of this class of sources.

Acknowledgments

NLP, LS, and PE acknowledge funding from INAF through the grants “Bando Ricerca Fondamentale INAF 2023” and “Bando Astrofisica Fondamentale 2024”. PE acknowledges support from the PRIN MUR SEAWIND (2022Y2T94C), funded by the European Union – Next Generation EU, Mission 4 Component 1 CUP C53D23001330006, and the INAF Grant BLOSSOM.

References

All Tables

All Figures

	Fig. 1. Light curves of 4U 0728-25, corrected for the background and the extraction region, in the S (top panel) and H (second panel) energy ranges, together with their ratio (third panel); the binning time is 100 s/bin. In the bottom panel the curve in the full energy range is reported, with a binning time of 50 s/bin.
In the text

	Fig. 2. Light curves of 4U 0728-25 in the full (bottom panel), S (top panel) and H (second panel) energy ranges, together with the HR (third panel), folded at the pulse period P = 103.301 s. The phase ranges selected for the phase-resolved spectral analysis are delimited with vertical dashed lines.
In the text

	Fig. 3. Pulse profile of 4U 0728-25 in the energy ranges 0.15–1.5, 1.5–2.5, 2.5–4, and 4–12 keV.
In the text

Fig. 4. Time-averaged spectrum of 4U 0728-25, where the black, red, and green symbols represent the pn, MOS1, and MOS2 data, respectively. Upper panel: superposition of the best-fitting PL + BB model, for the pn spectrum only, with the EPIC spectra. Middle panels: residuals (in units of σ) between data and model in the case of the single PL and of the single BB. Bottom panel: residuals in the case of the PL + BB.

In the text

Fig. 5. Long-term evolution of the luminosity in the energy range 0.2–12 keV (upper panel) and of the pulse period (lower panel) of 4U 0728-25. Telescope and epoch: (1) Uhuru catalogue: 12/12/1970–18/03/1973 (Forman et al. 1978); (2) Ariel V: 01/11/1974–30/04/1977 (Steiner et al. 1984); (3) Ariel V 3A catalogue: 15/10/1974–14/03/1980 (Warwick et al. 1981); (4) HEAO1: 19–23/10/1977 (Steiner et al. 1984); (5) HEAO1: 15/08/1977–15/02/1978 (Wood et al. 1984); (6) RXTE: 05/05/1997 (Corbet & Peele 1997); (7) AstroSat: 06–07/05/2016 (Roy et al. 2020); (8) XMM–Newton: 11/04/2023 (this work).

In the text