Probing the emission mechanism and nature of the pulsating compact object in the X-ray binary SAX J1324.4-6200

Recently, there has been renewed interest in the Be X-ray binary (Be/XRB) SAX J1324.4-6200 because of its spatial coincidence with a gamma-ray source detected by Fermi/LAT. To explore more thoroughly its properties, new observations were carried out in 2023 by NuSTAR, XMM-Newton, and Swift, jointly covering the energy range 0.2-79 keV. The X-ray spectrum of SAX J1324.4-6200 fits well with an absorbed power law with a high energy cut-off. We measured a NuSTAR spin period of 175.8127 +/- 0.0036 s and an XMM-Newton spin period of 175.862 +/- 0.025 s. All the available spin period measurements of SAX J1324.4-6200, spanning 29 years, are correlated with time, resulting in a remarkably stable spin-down of dP/dt=(6.09 +/- 0.06)*1E-9 s/s. If SAX J1324.4-6200 hosts an accretion powered pulsar, accretion torque models indicate a surface magnetic field of ~1E12-1E13 G. The X-ray properties emerging from our analysis strenghten the hypothesis that SAX J1324.4-6200 belongs to the small group of persistent Be/XRBs. We also performed radio observations with the Parkes Murriyang telescope, to search for radio pulsations. However, no radio pulsations compatible with the rotational ephemeris of SAX J1324.4-6200 were detected. We rule out the hypothesis that SAX J1324.4-6200 is a gamma-ray binary where the emission is produced by interactions between the pulsar and the companion winds. Other models commonly used to account for the production of gamma-rays in accreting pulsars cannot reproduce the bright emission from SAX J1324.4-6200. We examined other mechanisms for the gamma-ray emission and noted that there is a ~0.5% chance probability that an unknown extragalactic AGN observed through the Galactic plane may coincidentally fall within the Fermi/LAT error circle of the source and be the responsible of the gamma-ray emission. [Abridged]


Introduction
SAX J1324.4−6200 was discovered in 1997 with BeppoSAX (Angelini et al. 1998).Then, it was observed by ASCA (Angelini et al. 1998;Lin et al. 2002), Swift (Mereghetti et al. 2008), Chandra, and XMM-Newton (Kaur et al. 2009).It was never observed in X-rays at energies above 10 keV until 2023.SAX J1324.4−6200 was classified as a likely Be/X-ray binary (Be/XRBs, a sub-class of high mass X-ray binaries; for a review, see Reig 2011) at a distance in the range 1.5-15 kpc, based on two photometric analyses (Mereghetti et al. 2008;Kaur et al. 2009).No Gaia DR3 counterpart is catalogued within the Chandra 3σ error region (Gaia Collaboration et al. 2023).SAX J1324.4−6200hosts a pulsar with a spin period of P spin ≈ 173 s (Angelini et al. 1998).Since its discovery, it has shown a constant spin down with rate of Ṗ ≈ 6 × 10 −9 s s −1 (Kaur et al. 2009).A tentative detection of an orbital period at ≈ 27 hr was reported by Lin et al. (2002) in an ASCA observation.However, some doubts about this periodicity were raised by the same authors, who pointed out that the detection was based on only two cycles of the period, and the folded lightcurve had a suspicious not smoothed shape (see, also, the discussion in Mereghetti et al. 2008).The 1 − 10 keV flux has not shown so far large variability, always being ≈ 0.5 − 1 × 10 −11 erg cm −2 s −1 .The X-ray spectrum has been described by an absorbed power law, with N H ≈ 5 × 10 22 cm −2 and photon index ∼ 1.25 (Angelini et al. 1998;Lin et al. 2002;Mereghetti et al. 2008;Kaur et al. 2009).
The steady and relatively low X-ray luminosity of SAX J1324.4−6200(≈ 10 33−35 erg s −1 , assuming 1.5 d 15 kpc) is compatible with those displayed by the small subgroup of persistent Be/XRB pulsars (Pfahl et al. 2002), whose most famous member is X Persei.The relatively faint and stable X-ray emission of these sources is thought to be due to the wide and circular neutron star (NS) orbit (P orb ≫ 30 d, e < 0.2) and accretion of the NS from the low-density polar wind of the companion star (see, e.g., Mereghetti et al. 2008;Reig & Roche 1999;La Palombara & Mereghetti 2007;La Palombara et al. 2021, and references therein).
The nature of SAX J1324.4−6200 as an accreting NS was thought to be settled until Harvey et al. (2022) reported the discovery of a persistent γ−ray emission over 12.5 years of Fermi/LAT data from a region consistent with the position of SAX J1324.4−6200.In addition, they found evidence of variability in the γ-ray source, with relatively higher γ-ray emission over an 18-month period in 2018 and 2019.They ruled out that this emission was caused by other already catalogued γ-ray sources in the vicinity of SAX J1324.4−6200.There are no γ−ray pulsars in the third Fermi/LAT catalog associated with SAX J1324.4−6200(Smith et al. 2023), nor other X-ray binaries (Avakyan et al. 2023;Neumann et al. 2023;Fortin et al. 2023Fortin et al. , 2024)), nor other bright X-ray sources (search performed in Vizier 1 ).The source has no associated persistent radio emission (Harvey et al. 2022).Taking into account only the photons detected during the 2018-2019 bright event, Harvey et al. (2022) found that the γ-ray source has an offset from SAX J1324.4−6200 of ∼ 0.07 • , well within the 95 per cent containment radius (∼ 0.19 • ).They concluded that the 18month excess is spatially coincident with SAX J1324.4−6200.The γ-ray source was bright enough to allow spectral analysis.Harvey et al. (2022) found that the best fitting model is a power law with spectral index −2.43,and that the 100 MeV−500 GeV flux is ∼ 2.98 × 10 −6 MeV cm −2 s −1 .Harvey et al. (2022) pointed out that if the γ-ray emission is associated with SAX J1324.4−6200, it could be produced by the collision between the winds from the pulsar and the companion star, similar to what has been proposed for some γ-ray binaries (see, e.g., Dubus 2013;Paredes & Bordas 2019;Dubus et al. 2017;Chernyakova & Malyshev 2020, and references therein).This hypothesis was supported by the γ-ray luminosity and spectral index measured by Harvey et al. (2022), which are consistent with those of other γ-ray binaries with pulsar (see, e.g., Dubus 2013).This scenario casts serious doubts about accretion as the mechanism to explain the X-ray emission from SAX J1324.4−6200, as it was previously accepted.
In light of these recent findings, we gathered more information about SAX J1324.4−6200, to understand its nature.
Here we report on the analysis of X-ray observations of SAX J1324.4−6200spanning the 0.2-79 keV band and carried out in 2023 with NuSTAR, XMM-Newton, and Swift.To gain further insights, an observation was also performed with the Parkes Murriyang telescope in search for radio pulsations.Then, we discuss the nature of SAX J1324.4−6200 using the results of our Xray and radio data analysis together with the other information available for this source.

NuSTAR
The Nuclear Spectroscopic Telescope Array (NuSTAR) satellite is equipped with two identical co-aligned telescopes with focal plane modules FPMA and FPMB.Both operate in the 3−79 keV energy band (Harrison et al. 2013).NuSTAR observed SAX J1324.4−6200 from 1 to 3 July 2023, for a net exposure time of about 62 ks (see Table 1).Data were reduced using NUSTARDAS v2.1.2,which is part of HEASOFT v6.32.1, and the calibration files distributed with the CALDB v20230802 (Madsen et al. 2022).We extracted the source events from a circular region centred on it and with a radius of 87 ′′ and 75 ′′ for 1 https://vizier.cfa.harvard.edu/viz-bin/VizieRFPMA and FPMB, respectively.These radii were calculated to have the maximum signal-to-noise ratio, with the constraint to keep the extraction area within the same detector (Det-0).For the background, we extracted the events from circular regions located on the same detector of the source (Det-0) but in a zone of the focal plane free from the emission of SAX J1324.4−6200.Event times were corrected from satellite frame to the solar system barycenter using the ephemeris JPLEPH.200using the barycorr task.The spectra were rebinned using the optimal binning method by Kaastra & Bleeker (2016) and to have at least 25 counts per bin to enable the use of the χ 2 as fit stastistic.This variation of the Kaastra & Bleeker (2016) technique is implemented in the HEASOFT tool ftgrouppha.

XMM-Newton
The X-ray Multi-Mirror Mission (XMM-Newton) hosts the European Photon Imaging Camera (EPIC).It comprises the pn, Metal Oxide Semi-conductor 1 and 2 (MOS1 and MOS2) CCDs, operating in the 0.2−12 keV energy band (Jansen et al. 2001;Strüder et al. 2001;Turner et al. 2001).XMM-Newton observed SAX J1324.4−6200 on 12 July 2023 (see Table 1).We reduced the data using the XMM-Newton Science Analysis System (SAS, v21.0.0), with the latest calibration files available in the XMM-Newton calibration database (CCF).Calibrated event lists for the pn, MOS1, and MOS2 were obtained from the raw data with the SAS tasks epproc and emproc.For the pn, we used single-and double-pixel events, while for the MOS data we used single-to quadruple-pixel events.We excluded time intervals where the background was too high for a meaningful spectral analysis using the standard procedures 2 .The net exposure times obtained for the XMM-Newton observation are reported in Table 1.Source events were extracted using circular regions centred on the best known position of SAX J1324.4−6200.The radii of these extraction regions were r pn = 36.45′′ , r MOS1 = 50 ′′ , and r MOS2 = 51 ′′ .These radii were calculated with the SAS task eregionanalyse to have the maximum signal-to-noise ratio.Background events were accumulated for each of the three cameras using extraction regions in the same CCD where SAX J1324.4−6200 is located and not contaminated by its emission.The effective area of the pn, MOS1, and MOS2 were corrected in accordance with the CCF Release Note XMM-CCF-REL-3883 to improve the alignment with the NuSTAR spectra.Event times were corrected from satellite frame to the solar system barycenter using the ephemeris DE200 using the SAS task barycen.The spectra were rebinned with the tool ftgrouppha using the optimal binning method and to have at least 25 counts per bin to enable the use of the χ 2 as fit stastistic.

Swfit/XRT
SAX J1324.4−6200 was observed by the X-ray Telescope (XRT, 0.3-10 keV) on board of The Neil Gehrels Swift Observatory on 1 July 2023, for about 5.4 ks (see Table 1).The data were processed using the standard software (HEASOFT v6.31) and calibration (CALDB v20230725).Swift/XRT data were filtered with the task xrtpipeline (v0.13.7).No pile-up correction was necessary.We extracted the source and background events for the spectral analysis with xselect.For the source, we used a cir- cular extraction region centered on the source with radius of 18 pixels, while the background was extracted from an annular region centered on the same position, with inner and outer radii of 60 and 100 pixels.The instrumental channels were combined to include at least 20 photons per bin using the grppha task.The bandwidth was split into 1-MHz wide channels, which were 8-bit sampled every 1.024 millisecond, and only total intensity information was recorded.
The data were folded with the dspsr 4 (van Straten & Bailes 2011) pulsar package using the rotational ephemeris derived in this paper from the X-ray observations (see §3.1) and with a dispersion measure (DM) of 370 pc cm −3 .This initial DM was derived as the median between the minimum and maximum values obtained using the 1.5-15 kpc distance range mentioned above and the NE2001 and YMW16 models for the distribution of the electrons in the interstellar medium (Cordes & Lazio 2002;Yao et al. 2017).Given the extremely long period of the pulsar, the dispersive delay is negligible even in the worst case of a true DM being double the maximum predicted by the electron density models (1600 pc/cm 3 , in which case the dispersive delay across the band would be about 14 seconds, i.e. less than a 1/10 of the pulsar period).Notwithstanding, we searched the data over a DM range from 0 to 1600 pc/cm 3 to maximise the signal-tonoise ratio of a possible pulsed emission.This DM search was done jointly with a spin period search spanning ±350 ms around the nominal period predicted at the epoch of the radio observations using the X-ray ephemeris.This was carried out with pdmp from the software package psrchive 5 (Hotan et al. 2004;van Straten et al. 2012).

Timing and variability analysis
The NuSTAR data (modules A and B) and the XMM-Newton pn data were used for the timing analysis.The source and background events were extracted in the same regions adopted for the spectral extraction.We searched for the spin period in the 0.2 − 12 keV events of pn and in the 3 − 60 keV combined events of modules A and B of NuSTAR, using a Rayleigh test Z 2 with three harmonics (see, e.g., Buccheri et al. 1983).We refined the measurement of the detected signal using the phase-fitting method (see, e.g., Dall'Osso et al. 2003).
Using this technique, we measured a NuSTAR spin period of P spin,NuSTAR = 175.8127± 0.0036 s and an XMM-Newton spin period of P spin,XMM = 175.862± 0.025 s (see Table 2).The XMM-Newton period is consistent with the NuSTAR measurement at ∼ 1.97σ confidence level.We fitted all the previous and our new measured values of the spin period (Angelini et al. 1998;Lin et al. 2002;Mereghetti et al. 2008;Kaur et al. 2009) with the linear relation P(t) = P 0 + Ṗt.The data are highly correlated: the Pearson's linear coefficient is r = 0.9998, and the null hypothesis probability is 1.7 × 10 −9 (see Fig. 3).We obtained a spin period derivative of Ṗ = 6.09 ± 0.06 × 10 −9 s s −1 .The ASCA (Lin et al. 2002) and our recent NuSTAR spin period measurement are outliers from the general trend shown in Fig. 3.They could indicate the sporadic occurrence of small fluctuations of the accretion torque rate, or they might be caused by systematic effects in the spin determination.To produce the lightcurves, the events were rebinned with a time resolution of 0.2 s for pn and 0.5 s for NuSTAR.The final lightcurves were then created by subtracting the background.We folded the NuSTAR and pn lightcurves using their respective spin period.The pulse profiles, for different energy bands, are shown in Figures 1 and 2. We selected these energy bands to have a similar number of average counts.For each pulse profile, we calculated the pulsed fraction using the formula: where R max and R max are the maximum and minimum value of the rates (c/s) of the folded lightcurves.The measured values of p f are in Table 3.We investigated the short-term spectral variability using the hardness ratios (HRs).The HRs were defined as counts in the hard band divided by counts in the soft band, H/S , where the energy bands are: S = 0.2 − 5.2 keV and H = 5.2 − 12 keV for XMM-Newton, and S = 3 − 9 keV and H = 9 − 60 keV for NuSTAR.The bin times were defined as three times the spin period.Figure 4 does not show any significant variability in hardness within each observation.

Spectral analysis
Swift/XRT and NuSTAR data were collected during overlapping periods, while the XMM-Newton observation was made after 11 days.To investigate possible significant spectral variability between the two different data collection periods, and to verify that all the spectra can be fitted simultaneously, we proceeded as follows.First, we fit simultaneously Swift/XRT and NuSTAR data and afterwards XMM-Newton data (pn, MOS1, and MOS2).In the first case (XRT+NuSTAR), we obtained an acceptable fit with an absorbed power law with a high energy cutoff (po*highecut in XSPEC6 ).To model the photoelectric absorption, we used the tbabs model in XSPEC, and we set the abundances to those of the interstellar medium (wilm in XSPEC; Wilms et al. 2000).Renormalisation constant factors were included in the spectral fitting to account for intercalibration uncertainties between instruments.Due to the limited energy coverage of the XMM-Newton spectra, they were insensitive to the high energy cutoff parameters.To investigate possible spectral variability in the soft part of the XMM-Newton spectrum with respect to XRT+NuSTAR data, we froze E c and E f to the values found from the XRT+NuSTAR bestfit.The results are shown in Table 4 and Fig. 5.A close look at the best fit parameters in Table 4 shows that N H and Γ from the XRT+NuSTAR fit are slightly different compared to those from XMM-Newton.However, a comparison of the contour plots of Γ versus N H (Fig. 8), shows that these parameters are consistent between the two fits within 2σ.The reduced χ 2  of the XRT+NuSTAR and XMM-Newton cases are suspiciously high (XRT+NuSTAR: χ 2 red = 1.21;XMM-Newton: χ 2 red = 1.15).Therefore, we rebinned the spectrum for the purpose of visual inspection only, to improve the clarity of the residuals panel.This rebinning, obtained using setplot rebin in XSPEC, does not affect the results of the spectral fitting, which are obtained using the original binned spectrum.Using this approach, the new residual panels (see Fig. 6) show a wave-like structure between ∼ 3 keV and ∼ 20 keV, with a minimum at ∼ 7 keV.Given the general agreement between the best fit parameters obtained from these two fits, we decided to fit all datasets simultaneously to further improve the statistics.To explore pos- Model: const*tbabs*po*highecut. a : the absorbed X-ray flux (in units of 10 −12 erg cm −2 s −1 ) is calculated (using cflux in xspec) in the range 4 − 10 keV.sible better descriptions of the observed spectrum, we tried to fit it with physical models typically used for accreting NSs in Be/XRBs: a power law with Fermi-Dirac cutoff (fdcut, Tanaka 1986), the negative-positive exponential model (npex, Mihara 1995), Comptonization of soft photons in a hot plasma (comptt, Titarchuk 1994), double comptt (see, e.g., Doroshenko et al. 2012).All these models give similar or worse χ 2 /d.o.f.than po*highecut, and no substantial improvement in the residuals.
Therefore, we tried to get a better fit by adding new components to the base model po*highecut.We obtained better fits, both for χ 2 /d.o.f. and for the residuals, by adding a blackbody component or a Gaussian in absorption.The results of these fits are shown in Table 5 and Fig. 7.In Table 5, the cross-normalization constants of the XMM-Newton instruments with respect to NuSTAR reflect the slightly lower X-ray flux of SAX J1324.4−6200 during the XMM-Newton observation (see Table 4).To evaluate the chance probability of improvement of the fit by adding the blackbody or Gaussian in absorption component we simulated, for each case, 5 × 10 4 datasets without the extra components using the XSPEC routine simftest (see, e.g., Ducci & Malacaria 2022, and references therein).We found that the probability that the observed data are consistent with the model without the extra component is < 0.002% for both cases of Gaussian in absorption and blackbody as extra components.The radius of the blackbody component ranges between ∼ 0.07 km and ∼ 0.7 km, depending on the distance of the source (d = 1.5 − 15 kpc).The blackbody emission could be interpreted as the thermal radiation from a hot X-ray spot on the NS surface.In the other case, the Gaussian in absorption could be interpreted as a cyclotron resonance scattering feature (CRSF).In this case, using the law E cyc ≈ 11.6B 12 (keV) (e.g., Staubert et al. 2019), which links the centroid energy of the fundamental CRSF with the magnetic field strength of the pulsar (B 12 is the magnetic field   7. XMM-Newton (black: pn; red: MOS1; green: MOS2), NuSTAR (blue: module A; cyan: module B), and Swift/XRT (orange) spectra of SAX J1324.4−6200(rebinned using setplot rebin in XSPEC).Top panel: spectra are fitted simultaneously with an absorbed power law with high energy cutoff and a Gaussian in absorption at ∼ 6.9 keV.Second panel: Residuals relative to the top panel.Third panel: residuals of the fit of the spectra with an absorbed power law with high energy cutoff and a blackbody.Bottom panel: residuals of the fit of the spectra with an absorbed power law with high energy cutoff.See Table 5 for the best-fit parameters.
strength in units of 10 12 G), we infer a NS surface magnetic field strength of B ≈ 6 × 10 11 G.

Radio pulsations search
No radio pulsations compatible with the rotational ephemeris of SAX J1324.4−6200 were found in the Murriyang data.The nondetection can be used to calculate an upper limit of the pulsed flux density of the pulsar.We do so using the radiometer equation adapted for pulsars (see, e.g., Lorimer & Kramer 2004), adopting a system temperature of 21 K, an antenna gain of 0.735 K/Jy, using a pulse duty cycle of 10%, and a signal-to-noise ratio of 8.With such parameters, we find that SAX J1324.4−6200 must have a radio pulsed mean flux density lower than 9 µJy.

γ−ray binary scenario
We begin by considering the scenario proposed by Harvey et al. (2022): SAX J1324.4−6200 is a γ−ray binary, where the accretion is turned off and the observed emission is produced by the collision of the pulsar wind with the stellar wind of the  4. companion star.In γ−ray binaries, most of the rotational energy of the pulsar is thought to be carried away by the pulsar wind, and a fraction of this energy is released as radiation when this wind interacts with the outflow from the companion star.The loss of rotational energy from the pulsar is therefore the energy reservoir that is used to produce the observed radiation (see, e.g., Dubus 2013).This can be estimated as Ėrot = (2π) 2 I Ṗ/P 3 ≈ 7.2 × 10 31 erg s −1 , where I is the NS moment of inertia calculated assuming M NS = 1.4M ⊙ and R NS = 12 km, and we adopted the spin period P detected with NuSTAR and the spin period derivative Ṗ obtained using all the X-ray observations of SAX J1324.4−6200(Sect.3.1).The distance of SAX J1324.4−6200 is in the ranges 1.5-8 kpc or 5-15 kpc according to Kaur et al. (2009) and Mereghetti et al. (2008), respectively.Therefore, given an observed 0.2-60 keV unabsorbed flux of ∼ 1.8 × 10 −11 erg cm −2 s −1 , the observed Xray luminosity (L x,obs ≈ 4.8 − 480 × 10 33 erg s −1 ) exceeds the spin-down energy loss by orders of magnitude, ruling out the colliding winds mechanism as the origin of the source emission.

Non-accreting magnetar scenario
Retaining the non-accretion scenario considered in Sect.4.1, the magneto dipole formula for pulsars (see, e.g., Manchester & Taylor 1977) 2022;Chandrasekhar & Fermi 1953).In light of this, we decided to also explore the non-accreting magnetar scenario.Recently, it has been proposed that the γ−ray binary LS 5039 contains a nonaccreting magnetar (Yoneda et al. 2020).For this source it was hypothesised that it can tap into its magnetic energy to produce the observed emission.In particular, Yoneda et al. (2020) outlined the hypothesis that the interaction between the stellar wind of the companion star and the magnetar magnetosphere, leads to magnetic reconnections that can release magnetic energy as radiation.Following the same line of reasoning as Yoneda et al. (2020), we find that if the compact object in SAX J1324.4−6200 is also a magnetar, the available energy budget from the magnetic energy would be the dissipation of the magnetic energy Ėmag = R 3 NS B 2 /(3τ) erg s −1 (see, e.g., Dall'Osso et al. 2012).Assuming a magnetic field strength of ∼ 10 15 G (which is 10 times smaller than the value obtained from the magneto dipole formula for pulsars but closer to the typical values encountered in magnetars) and a magnetic field decay time-scale dominated by the so-called Hall term (τ ≈ 10 3 − 10 4 yrs, see, e.g., Dall'Osso et al. 2012;Viganò et al. 2012Viganò et al. , 2013;;Turolla et al. 2015), the magnetic energy loss would be Ėmag ≈ 10 35 − 10 36 erg s −1 .This Ėmag value is promisingly high for the case of SAX J1324.4−6200.However, the relatively low τ adopted here does not guarantee a magnetar-scale field and a magnetic energy actively dissipated in a magnetar with an age of ∼ 10 5 yrs (see, e.g., Dubus 2013).This age corresponds to the most plausible evolutionary stage of a young binary system with a pulsar and a Be star.The presence of magnetars with these properties in γ−ray binaries and high-mass X-ray binaries is currently an intense matter of debate (see, e.g., Bozzo et al. 2008;Xu et al. 2022;Yoneda et al. 2020;Torres et al. 2012;Popov 2023, and references therein).In light of this, if we assume τ ∼ 10 5 yrs, we obtain Ėmag ≈ 10 34 erg s −1 , which is still large enough to power the radiation output of SAX J1324.4−6200.

Accretion torque
We consider the case in which the X-ray luminosity is caused by the accretion of matter and the long-term spin-down is due to a torque generated by the exchange of angular momentum between the pulsar and the accreting matter.To assess quantitatively whether L x , P, and Ṗ are compatible with this hypothe-sis, and which magnetic field strength is required, we compared our measurements with the predictions of three well-established models: Ghosh & Lamb (1979), Wang (1996), andShakura et al. (2012).In the first two models, the accretion is mediated by a disc.In the Wang model, we assumed that the dissipation timescale of the toroidal component of the magnetic field is determined by the reconnection outside of the disc (see Wang 1995).For the inner disc radius, we adopted the prescription described in Bozzo et al. (2009).We set both η (a screening factor due to the currents induced on the surface of the accretion disc) and γ max (the maximum magnetic pitch angle), which appear in equation 19 in Bozzo et al. (2009), equal to 1.Note that these parameters are poorly known by the current theory, as discussed in Bozzo et al. (2009).Moreover, both the Ghosh & Lamb (1979) and Wang (1995) prescriptions for the NS magnetospheric radius and accretion torques tend to display unphysical behaviors at low mass accretion rates, such that results in this regime should be taken with caution (Bozzo et al. 2009(Bozzo et al. , 2018)).
In the third accretion-torque model considered here, the long-term spin-down shown by SAX J1324.4−6200 could be the consequence of the quasi-spherical accretion of the stellar wind from the companion star onto the surface of the NS.Shakura et al. (2012, and references therein) pointed out that, if a pulsar in a binary system is wind-fed and the X-ray luminosity is moderate/low (L x 10 36 erg s −1 ), a hot quasistatic and convective shell forms around the magnetosphere of the NS.This shell can mediate the transfer of angular momentum, and the NS can spin up or down, depending on the difference between the angular velocity of the matter near the magnetospheric boundary and that of the magnetosphere itself.In contrast to the free-fall accretion regime, which produces a more erratic torque reversals (see, e.g., Bildsten et al. 1997;Malacaria et al. 2020, and references therein), the mechanism proposed by Shakura et al. (2012) can explain the predominance of a long-term spin-up or spin-down, although it can be interrupted by episodic short events of opposite spin rate and other fluctuations (see, e.g., González-Galán et al. 2012).As noted in Postnov et al. (2014), for pulsars that exhibit long-term spindown (i.e., in non-equilibrium), it is possible, under some reasonable simplifications, to obtain a lower bound on their magnetic field strength from the spin-down rate, which can be expressed as: ωsd ≈ −0.75 × 10 −12 Π sd µ 13/11 30 Ṁ3/11 In the equation above, Π sd is a combination of dimensionless parameters of the theory: it varies from ∼ 4 to ∼ 10 ( Postnov et al. 2015).Here, we assume Π sd = 5 (due to the lack of information about the binary separation and properties of the stellar wind).Then, the magnetic dipole moment is expressed as µ 30 = µ/(10 30 G cm 3 ), Ṁ16 = Ṁ/(10 16 g s −1 ) is the mass accretion rate, and P 100 is the spin period in units of 100 s.
The results are summarized in Fig. 9.The plot shows the X-ray luminosity on the x-axis and the spin period derivative on the y-axis.The black rectangle corresponds to the measurements of Ṗ and the X-ray luminosity, the latter depending on the distance of SAX J1324.4−6200.The curves intersecting the black rectangle were obtained using the Ghosh & Lamb (1979), Wang (1995), andShakura et al. (2012) models, and represent the limiting solutions that can intersect the black rectangle.The magnetic fields of these solutions are B 10 12 − 10 13 G, which are typical of accreting pulsars in high mass X-ray binaries.The Ṗ − L x solutions obtained with the Ghosh & Lamb (1979) and Wang (1995) models for d = 15 kpc show steep slopes (see Fig. 9).This means that we would expect frequent torque reversals for small variation of the mass accretion rate.On the contrary, in Sect.3.1 we have shown that Ṗ of SAX J1324.4−6200 is remarkably stable since its discovery, i.e. for a time interval of ≥ 29 years.Also L x (a proxy of the mass accretion rate) did not vary significantly: SAX J1324.4−6200never showed an X-ray outburst nor flare (typical of most of Be/XRBs) in any observation, nor it was caught during one of such bright events by X-ray monitors such as INTEGRAL, Swift/BAT, RXTE/ASM.All these observational evidences suggest that the mass accretion rate on the pulsar must be highly stable.We note the small flux variability between the NuSTAR+XRT and XMM-Newton observations (see Table 4).This could be due to a fluctuation of the mass accretion rate and, consequently, of the accretion torque rate, and might therefore also account for the NuSTAR spin period measurement that deviates from the long-term spin-down trend shown in Fig. 3. Another important point is that the overall properties of the X-ray spectrum of SAX J1324.4−6200(a relatively hard X-ray spectrum, with a possible blackbody component with kT bb 1.1 keV, Sect.3.2) are in agreement with the spectra observed in accreting high-mass X-ray binaries (see, e.g., Kretschmar et al. 2019).
The pulsed fraction of SAX J1324.4−6200 reported in Table 3 is relatively high, in agreement with the typical pulsed fractions observed in other accreting X-ray pulsars with high magnetic field, although it is remarkably stable with the energy, while in the other accreting systems it shows a general increase with the energy (see, e.g.Lutovinov & Tsygankov 2009;Ferrigno et al. 2023).
There exists another persistent Be/XRB that exhibits an uninterrupted long-term spin-down similar to SAX J1324.4−6200:CXOU J225355.1+624336.Its spin-down has a rate comparable with that of SAX J1324.4−6200 and it has been ongoing for at least 21 years (La Palombara et al. 2021;Esposito et al. 2013).This similarity, together with the overall X-ray properties of SAX J1324.4−6200emerging from our X-ray analysis, strenghtens the hypothesis that SAX J1324.4−6200belongs to the elusive class of persistent Be/XRBs.

γ−ray emission from an accreting pulsar
Only a few other accreting pulsars in Be/XRBs have a γ−ray counterpart candidate (for a list of XRBs with a γ−ray counterpart candidate see, e.g., Ducci et al. 2023).Different types of models have been proposed to explain the possible γ−ray emission from accreting pulsars.Here, we consider two of them.The first model, developed by Bednarek (2009a,b), introduced concepts that were further developed by Papitto et al. (2014) (see also Torres et al. 2012;Papitto et al. 2012) to explain the γ−ray emission from the low-mass X-ray binary XSS J12270−4859.The model proposed by Bednarek (2009a,b) considers a binary system composed of an accreting and strongly magnetized NS and a companion massive star (OB type).Under certain conditions, the interaction between the NS magnetosphere and the dense wind from the donor star results in the creation of a magnetized and turbulent transition region around the pulsar.Within this region, electrons are accelerated to high energies, and subsequently produce γ−rays (up to few GeV) through synchrotron process and inverse Compton scattering, in response to the Xray radiation emitted from the NS surface.The maximum energy budget available for accelerating the electrons is limited by the energy that can be extracted by the matter interacting with the magnetosphere corotating with the NS.Bednarek (2009a) showed that the upper-limit for the available power for the ac- celeration of electrons is: where B 12 is the magnetic field strength at the polar cap of the NS, in units of 10 12 G, Ṁ16 is the mass accretion rate Ṁacc in units of 10 16 g s −1 , and η is the fraction of the power that can be effectively converted to relativistic electrons and subsequently to gamma-ray radiation.To calculate Ė from Eq. 3, we derived the mass accretion rate from the formula L x ≈ GM ns Ṁacc /R ns , assuming L x within the range ≈ 4.8 − 480 × 10 33 erg s −1 , M ns = 1.4 M ⊙ , and R ns = 12 km.Then, we adopted the corresponding limits for the magnetic field strength of the NS that we derived from the accretion torque models of Ghosh & Lamb (1979), Wang (1995), andShakura et al. (2012): B ≈ 10 12 − 9 × 10 13 G.
From Eq. 3, we obtain that the maximum power available for accelerating the electrons is Ė ≤ 4 × 10 31 ≪ L γ ≈ 10 33 erg s −1 .Ė is further reduced if we consider the efficiency conversion factor η ≈ 0.1 (Bednarek 2009a).Therefore, this mechanism cannot explain the observed γ−ray emission.
The second model we are considering here was originally proposed by Cheng & Ruderman (1989) (see also Bisnovatyi-Kogan et al. 1980) and subsequently improved by others (Cheng et al. 1991;Cheng & Ruderman 1991;Cheng et al. 1992;Romero et al. 2001;Orellana et al. 2007;Ducci et al. 2023).The key concept of this model is that γ−ray photons are the result of cascades initiated by the decay of π 0 , which originate from protons accelerated in the magnetosphere of a pulsar fed by an accretion disc.Here we consider the model version presented in Ducci et al. (2023) where the evolution of cascades inside and outside the accretion disc takes into account pair and photon production processes that involve interactions with nuclei, X-ray photons from the accretion disc, and the magnetic field.This model provides results above 10 GeV, so it cannot be directly compared with the Fermi/LAT detection reported by Harvey et al. (2022), which is in the 0.3 − 10 GeV energy range (see Fig. 10).Nevertheless, the model predictions can be compared with the Fermi/LAT upper limits at higher energies (10 − 300 GeV) reported in Harvey et al. (2022) and with the 5σ upper limit obtained from the H.E.S.S survey of the Galactic plane at energies > 1 TeV (H.E. S. S. Collaboration et al. 2018).The H.E.S.S. upper limit can be obtained from the sensitivity map reported in (H.E. S. S. Collaboration et al. 2018) at the position corresponding to that of the γ-ray source detected by Harvey et al. (2022).Following the model in Ducci et al. (2023), we simulated a grid of γ-ray spectra, expected for a source with the properties of SAX J1324.4−6200,i.e. with an X-ray luminosity in the range L x ≈ 5 − 500 × 10 33 erg s −1 , a magnetic field strength in the range 10 12 − 9 × 10 13 G, and spin period of 175 s.Among the possible solutions, we considered those for the "strong shielding" case.This is an approximation in which the X-ray photons produced by accretion at the stellar surface (including the accretion column) are strongly shielded by the accreting matter.In the case of "weak shielding", the overall reduction of the potential drop over the region where the protons are accelerated would, in the case of SAX J1324.4−6200,lead to solutions with γ-ray emission much lower than in the case of strong shielding and thus of much less interest, as we will see below.In general, all simulated spectra in the case of strong shielding give γ-ray fluxes well below the Fermi/LAT and H.E.S.S. upper limits (< 2 orders of magnitude).An extrapolation of the Harvey et al. (2022) detections to higher energies (i.e.above 10 GeV), or conversely an extrapolation of the simulated data to lower energies, shows that the observed emission is > 100 times brighter than that predicted by the model.To illustrate this result more clearly, we show in Fig. 10, as an example, a comparison between the observed data and the simulated spectrum obtained by assuming L x = 5 × 10 35 erg s −1 (and thus d = 15 kpc) and B = 4 × 10 13 G.The simulated spectrum takes into account the apparent flux increase due to a beaming factor of b f = 0.03, which is also provided by the simulations.The other simulated spectra, not shown here, display an analogous low γ-ray emission.Therefore, although a direct comparison between observed and simulated data below 10 GeV is not possible, based on the constraints given by the upper limits above 10 GeV, and assuming that the observed (simulated) spectrum can be extrapolated smoothly at higher (lower) energies (i.e.without unphysical sharp jumps in fluxes greater than two orders of magnitude), it is reasonable to conclude that also this mechanism cannot explain the intensity of the γ-ray emission reported in Harvey et al. (2022).

Extragalactic AGN scenario
In this section, we consider the hypothesis that the γ-ray source detected by Harvey et al. (2022) and positionally associated with SAX J1324.4−6200 could be an extragalactic AGN, observed through the Galactic plane.To test how likely is this hypothesis, we estimated the expected number of AGNs that Fermi/LAT could detect within a circle with a 95% radius of R 95 = 0.1859 • .We based our calculation on the number of AGNs reported in the Fourth LAT AGN Catalog (4LAC-DR3, Ajello et al. 2023)7 within ±5 • from the Galactic plane.We call the density of the AGNs detected by Fermi/LAT within this area ρ AGN .We considered only AGNs observed through the Galactic plane because we expect a lower density of these background sources than at higher Galactic latitudes, due to effects of absorption.We assumed that AGNs are uniformly distributed and that the circle πR 2 95 is small enough to neglect other effects (such as the curvature of the sky).Using the definition of the Poisson distribution, and given that the expected number of AGNs in 6200 was observed with the Ultra Wideband Low (UWL, Hobbs et al. 2020) receiver of the Murriyang radio telescope (Parkes, NSW, Australia) on 2023 March 30 for three hours.The data were recorded over the entire 3384 MHz bandwidth of the UWL receiver centered at a frequency of 2368 MHz.

Fig. 5 .Fig. 6 .
Fig. 5. Left panel: NuSTAR (black: module A; red: module B) and Swift/XRT (green) spectra of SAX J1324.4−6200,fitted simultaneously with an absorbed power law with high energy cutoff.Right panel: XMM-Newton (black: pn; red: MOS1; green: MOS2) spectra of SAX J1324.4−6200,fitted simultaneously with an absorbed power law with high energy cutoff (E c and E f frozen to the best fit values obtained from the XRT+NuSTAR fit).The lower panels show the residuals of the fit.

Fig.
Fig.7.XMM-Newton (black: pn; red: MOS1; green: MOS2), NuSTAR (blue: module A; cyan: module B), and Swift/XRT (orange) spectra of SAX J1324.4−6200(rebinned using setplot rebin in XSPEC).Top panel: spectra are fitted simultaneously with an absorbed power law with high energy cutoff and a Gaussian in absorption at ∼ 6.9 keV.Second panel: Residuals relative to the top panel.Third panel: residuals of the fit of the spectra with an absorbed power law with high energy cutoff and a blackbody.Bottom panel: residuals of the fit of the spectra with an absorbed power law with high energy cutoff.See Table5for the best-fit parameters.

Table 1 .
Summary of the X-ray observations.

Table 2 .
Measurement of the spin period in the NuSTAR and XMM-Newton observations analysed in this work.

Table 4 .
Best-fit spectral parameters from the simultaneous fit of the Swift/XRT and NuSTAR (XRT+NuSTAR) and XMM-Newton data sets with an absorbed power law with high energy cutoff, and a crosscalibration renormalization constants.

Table 5 .
Best-fit spectral parameters for the simultaneous fit of the Swift/XRT, NuSTAR, and XMM-Newton spectra with three different models.
provides an estimate for the magnetic field strength of SAX J1324.4−6200 of B ≈ 3 × 10 16 G.Such magnetic field would be stronger than the typical values inferred and measured from magnetars (see, e.g, Kaspi & Beloborodov 2017), but still below the maximal virial value for a NS of B max ≈ 10 18 G (Mushtukov & Tsygankov Fig.9.Ṗ − L x functions from theGhosh & Lamb (1979, GL79, dashed  lines),Wang (1995, W95, dot-dash lines), andShakura et al. (2012,  S12, dot lines) models corresponding to the minimum and maximum values of the magnetic dipole moment which provide solutions from these models which are in agreement with the observed spin period derivative of SAX J1324.4−6200 and its X-ray luminosity as a function of the distance (black rectangle).