3.1 Radio spectral properties

Fig. 2 shows the mean integrated flux density of PSR J0837 $-$ 2454 over the Parkes observing campaign. This pulsar shows a significant amount of variation in the flux density over both short and long timescales, which is likely indicative of both refractive and diffractive scintillation from the intervening interstellar medium (ISM). On the other hand, neither the DM nor RM show any evidence for significant variations over the same timespan. Using a median $\textrm{RM} = 360.72$ rad m $^{-2}$ and a DM of 143.1 pc cm $^{-3}$ , we can derive an average magnetic field strength along the line of sight,

\begin{equation} \displaystyle B_{\parallel} = 1.232 \frac{\textrm{RM}}{\textrm{DM}} = 3.1 \, \mu \textrm{G}.\end{equation}

To measure the spectrum of PSR J0837 $-$ 2454, we do not require the full native resolution of 1 MHz offered by the UWL receiver. Thus we downsampled our data to 8 frequency bins and calculated the mean integrated flux density in each bin. The spectrum calculated with this method for every epoch is shown in Fig. 3 using transparent square markers. Given the large flux density variation as a function of time, we also averaged the per-epoch spectrum to get the mean spectrum, which is shown using solid purple markers in Fig. 3, where the uncertainty is given by the standard deviation of the measured flux density in each bin. A power-law fit ( $S_{\nu} \propto \nu^{\alpha}$ ) to this mean spectrum yields a spectral index of $\alpha = 0.13 \pm 0.6$ , though we note that this flat spectrum is likely due to effects of refractive scintillation, similar to that seen in, e.g. PSR J1740+1000 (McLaughlin et al. Reference McLaughlin, Arzoumanian, Cordes, Backer, Lommen, Lorimer and Zepka2002).

Figure 2. Mean integrated flux density over time from Parkes Murriyang UWL observations of PSR J0837 $-$ 2454 starting on 2021 October 3 (MJD 59490). Errors shown are 1 $\sigma$ .

Figure 3. Radio spectrum from Parkes Murriyang UWL observations of PSR J0837 $-$ 2454. Light colored points indicate mean integrated flux densities for each epoch of observation, while dark purple points indicate the average at each frequency. Errors shown are 1 $\sigma$ .

3.2 Updated radio timing model

Using our radio data, we also generated a new timing solution for PSR J0837 $-$ 2454. Since there was no significant variation in the DM relative to the fiducial DM used in the fold-mode observations, we held the DM fixed to this value. Given that there are no significant DM variations, we time and frequency scrunched each observation and, along with the standard profile from Pol et al. (Reference Pol2021), derived a corresponding time-of-arrival (TOA) for each epoch using pat.

The timing model and the corresponding residuals are shown in Table 1 and Fig. 4, respectively. We used PINT (Luo et al. Reference Luo2019, Reference Luo2021) to time the pulsar, including modeling the white (EFAC only) and red noise (Susobhanan et al. Reference Susobhanan2024). We used the timing model from Pol et al. (Reference Pol2021) to initialize our timing model after which the parameters were allowed to vary. We tried including more frequency derivatives, parallax, and proper motion to the timing model, but none of these parameters were found significant via an f-test.

Similar to Pol et al. (Reference Pol2021), we notice significant red noise in the residuals for this pulsar. Without an explicit red noise process in the timing model, the red noise was absorbed in the fits to the sky-position of the pulsar, resulting in underestimating the uncertainty on the measured right ascension and declination for the pulsar by $\sim$ 4 and $\sim$ 6, respectively. To accurately reflect the uncertainties on these parameters, we include a red noise model with only a single harmonic, as determined using the WaveX model in PINT (Susobhanan et al. Reference Susobhanan2024). We also attempted to phase connect the TOAs from these observations with those collected in Pol et al. (Reference Pol2021), but we were unsuccessful due to the large red noise present in this pulsar.

3.3 VLBA astrometric results

The 7–13 $\sigma$ detections described in Section 2.1.2 enable us to measure the position of PSR J0837 $-$ 2454 with a nominal accuracy of 0.5–1 milliarcsec in right ascension and 1–2 milliarcsec in declination at each epoch. However, these statistical uncertainties do not encompass any residual calibration errors, which can be substantial, particularly at low observing elevations as is necessary for the VLBA at these declinations (e.g. Deller et al. Reference Deller2019). In this instance, the significant angular broadening of the calibrator sources also means that longer integration times are required to achieve sufficient signal-to-noise ratio on the calibrators, which further degrades the astrometric quality as the interpolation time increases.

Table 1. Parameters for PSR J0837 $-$ 2454.

Figure 4. Post-fit residuals, along with 1 $\sigma$ uncertainties, from the Parkes timing campaign starting on 2021 October 3 (MJD 59490). Red noise is clearly visible, with the overall RMS of 100 $\mu$ s in the residuals being comparable to the timing solution obtained in Pol et al. (Reference Pol2021). We are unable to phase connect the TOAs from this observing campaign to those from Pol et al. (Reference Pol2021).

In Table 2, we show the results of a straightforward least squares fit to the reference position, proper motion, and parallax for PSR J0837 $-$ 2454 based on the VLBA results, and we plot the model against the observed positions as a function of time in Fig. 5. However, significant unmodelled systematic errors are present in the data, as can be seen from the $\chi^2$ of the fit. We discuss the degree to which these results (and a future extension to the VLBA astrometry) can constrain the distance to the pulsar in Section 4.1. Finally, we note that the uncertainty in the reference position of the pulsar, which is not quoted in Table 2, is dominated by the uncertainty in the in-beam calibrator positions, which themselves were referenced to the primary phase calibrator located at a larger angular separation – this uncertainty is likely to be on the order of 5 milliarcsec in the current dataset.

Table 2. VLBA astrometric results for PSR J0837 $-$ 2454.

Figure 5. The offset of PSR J0837 $-$ 2454 from the fitted reference position in right ascension (top) and declination (bottom) as a function of time. As discussed in the text, the error bars are from the image plane fit to the pulsar position only and are underestimated (particularly in declination).

3.4 X-ray pulsation search and pulse profile

Fig. 6 shows power spectra of the pn data generated using powspec in FTOOLS (Blackburn Reference Blackburn, Shaw, Payne and Hayes1995). A peak at 1.58845 Hz (629.545 ms) is evident and corresponds to the spin frequency of PSR J0837 $-$ 2454. Note that Pol et al. (Reference Pol2021) measured a spin frequency of 1.58878890 Hz ( $P=629.41024\mbox{ ms}$ ) and spin frequency time derivative of $({-}8.808\pm0.004)\times10^{-13}\mbox{ Hz s}^{-1}$ [ $\dot{P}=(3.490\pm0.002)\times10^{-13}\mbox{ s s}^{-1}$ ] at MJD 55588.628331, such that the expected spin frequency at the time of the XMM-Newton observation is $1.5884\pm0.0007$ Hz ( $629.5\pm0.4$ ms), which matches that from the power spectrum here.

Figure 6. Power spectra from pn data for 0.5–1 keV (top) and 0.2–10 keV (bottom). The peak at 1.58845 Hz is the spin frequency of PSR J0837 $-$ 2454.

Next, we perform acceleration searches using PRESTO (Ransom, Eikenberry, & Middleditch Reference Ransom, Eikenberry and Middleditch2002). Data are folded at the candidate pulse frequency using prepfold, and a refined frequency is determined. We find the strongest pulsations at a frequency of 1.58845213(9) Hz or period of 629.54368(4) ms (at MJD 60059.0334841) in the energy range 0.5–1 keV. We then use photonphase in PINT (Luo et al. Reference Luo2021) to calculate the rotation phase for each photon. The resulting pulse profiles in various energy bands are shown in Fig. 7. Pulsations are only clearly evident at $E<1\mbox{ keV}$ . For the sinusoidal pulse profiles at 0.2–10 keV and 0.5–1 keV, the pulsed fractions [ $=\mbox{(max-min)/(max+min)}$ ] are 17% and 39%, respectively.

Figure 7. Pulse profile of PSR J0837 $-$ 2454, obtained by folding the pn data in the energy ranges 0.2–0.5, 0.5–1, 1–2, 2–10, and 0.2–10 keV at the 629.544 ms spin period. Errors shown are 1 $\sigma$ .

3.5 X-ray spectral analysis

Fig. 8 shows the combined MOS (MOS1+MOS2) and pn spectra. We perform model fits to these spectra using Xspec 12.14.0 (Arnaud Reference Arnaud, Jacoby and Barnes1996). For the spectral model, we include several components. First, we use constant to account for a possible instrumental difference between MOS and pn spectral normalisations, and we fix its value to 1 for the pn spectrum and allow it to vary for the combined MOS spectrum. Next, we include a component to account for photoelectric absorption by the interstellar medium, i.e. tbabs with abundances from Wilms, Allen, & McCray (Reference Wilms, Allen and McCray2000) and cross-sections from Verner et al. (Reference Verner, Ferland, Korista and Yakovlev1996). To model the intrinsic spectrum of PSR J0837 $-$ 2454, we first consider a single component comprising a power law (PL; powerlaw), blackbody (BB; bbodyrad), or neutron star atmosphere. For the last, we use a magnetic partially ionized hydrogen atmosphere model (nsmaxg; Ho, Potekhin, & Chabrier Reference Ho, Potekhin and Chabrier2008; Ho Reference Ho, Petit, Jardine and Spruit2014) and fix the model parameters of neutron star mass and radius to $M=1.4\,{\rm M}_\odot$ and $R=10\mbox{ km}$ , respectively. We also tried a non-magnetic fully ionized hydrogen atmosphere model (nsatmos; Heinke et al. Reference Heinke, Rybicki, Narayan and Grindlay2006), a non-magnetic partially ionized hydrogen, helium, or carbon atmosphere model (nsx; Ho & Heinke Reference Ho and Heinke2009), and a magnetic partially ionized carbon, oxygen, or neon atmosphere model (Mori & Ho Reference Mori and Ho2007), but none of these provide good fits and thus will not be discussed further.

Figure 8. Spectra of PSR J0837 $-$ 2454 from MOS1+MOS2 and pn data. Top panel shows data with 1 $\sigma$ errors (crosses) and spectral model (dashed lines for each component). Bottom panel shows $\chi^2=\mbox{(data-model)/error}$ . The spectral model shown here is a $B=10^{13}\mbox{ G}$ thin partially-ionized hydrogen atmosphere with two temperature components; a two temperature blackbody model yields similar results (see text and Table 3 for details).

The results of fits with a single component model for the intrinsic spectrum of PSR J0837 $-$ 2454, i.e. either power law, blackbody, or atmosphere model, are given in Table 3. We see that the power law model not only has the worst fit of the three, but the best-fit power law index $\Gamma\approx 8$ is unrealistic. Meanwhile, thermal models provide better fits if the distance to the pulsar is taken to be less than a few kpc. For example, in order for the best-fit blackbody emission radius $R^\infty$ [ $=R(1-2GM/c^2R)^{-1/2}$ ] to be less than $\sim$ $19\mbox{ km}$ (corresponding to a maximum neutron star radius $R\approx 13-15\mbox{ km}$ ), then the distance must be $d<4\mbox{ kpc}$ (after accounting for gravitational redshift). Similarly assuming that the entire stellar surface emits at a single uniform temperature, the best-fit atmosphere model gives a distance of only 90 pc. We note that nsmaxg models with $B=2$ , 4, or $7\times10^{12}\mbox{ G}$ yield fits that are about as good as that of $10^{13}\mbox{ G}$ and have best-fit parameters $N_\textrm{H}$ increasing, $T_\textrm{eff}^\infty$ decreasing, and d decreasing with lower values of B, while fits are poor for models with lower and higher magnetic fields, i.e. $10^{10}$ – $10^{12}\mbox{ G}$ and $(2-3)\times10^{13}\mbox{ G}$ . The absorbed 0.3–10 keV flux is $7.2\times10^{-14}\mbox{ erg s}^{-1}$ cm $^{-2}$ , so that the unabsorbed 0.3–10 keV X-ray luminosity $L^\infty_\textrm{X}\sim3.3\times10^{31}\mbox{ erg s}^{-1}(d/0.9\mbox{ kpc})^2$ and bolometric luminosity $L^\infty\sim7.6\times10^{31}\mbox{ erg s}^{-1}(d/0.9\mbox{kpc})^2$ .

Table 3. Results of spectral modeling.

Errors are $1\sigma$ .

^aXspec models: PL=powerlaw, BB=bbodyrad.

^bNeutron star mass $M=1.4\,{\rm M}_\odot$ and radius $R=10\mbox{ km}$ are fixed, such that a redshift of 1.3 is used to convert between unredshifted and redshifted values. Distance is allowed to vary with fit, while fraction of surface emitting is fixed to 1. Magnetic field $B=10^{13}\mbox{ G}$ is assumed (see text for other values).

Next, we consider two component models for fitting the measured spectra, on top of constant and tbabs, in particular, two power laws (PL+PL), a blackbody and a power law (BB+PL), or two blackbodies (BB+BB). For two power laws, the best-fit does not yield an improvement in $\chi^2$ , and the second power law only has an upper limit on its normalisation. Meanwhile, the best-fit results for BB+PL and BB+BB are given in Table 3. While we see an improvement in $\chi^2$ for BB+PL compared to only a power law or only a blackbody, with the power law providing the primary flux at $>$ $1\mbox{ keV}$ , the power law index is still a high $\Gamma>5$ . Meanwhile, two blackbodies yield a much improved fit ( $\Delta\chi^2=25.4$ for two extra degrees of freedom, so that a f-test gives a probability of $3\times10^{-6}$ ). The second blackbody with a higher temperature and small emitting region could be interpreted as emission from a hot spot on the neutron star surface that rotates to generate the pulsations detected and described in Section 3.4. The size of the colder blackbody limits the distance to PSR J0837 $-$ 2454 to be $d<1.7\mbox{ kpc}$ . We note for completeness that we also fit the spectra with nsmaxg+PL, with best-fit results shown in Table 3. The small improvement of the best-fit over only nsmaxg shows that the addition of the power law (to better fit the spectra at high energies) is not really needed. This difference with the BB+PL results is because neutron star atmosphere model spectra are harder at high energies than blackbody spectra at the same temperature due to the atmospheric opacity decreasing with energy, such that higher energy photons come from deeper, hotter layers of the atmosphere (see, e.g. Pavlov et al. Reference Pavlov, Shibanov, Zavlin, Meyer, Alpar, Kiziloglu and van Paradijs1995). This also implies that fitting spectra with a single blackbody gives a higher temperature and smaller emitting region compared to fitting with an atmosphere model, as evident from the results given in Table 3.

Finally, we fit the X-ray spectra using a two component neutron star atmosphere model, like the two component blackbody model described above. Such a model with emission from a small hot magnetic polar cap and emission from the remaining cooler surface can produce rotation modulated pulsations. A more robust distance estimate can also be determined. Since a single atmosphere model provides a good fit to the observed spectra without a need for a second component and guided by the good fit of the two blackbody model, the model we use here is a thin atmosphere model. A thin atmosphere is one whose physical thickness is such that the atmosphere is optically thick ( $\tau>1$ ) at low energies and optically thin ( $\tau<1$ ) at high energies. Thus the observed high energy photons arise from a condensed surface beneath the gaseous atmosphere, and this condensed surface emits blackbody-like emission (away from features such as electron and ion cyclotron resonances; see, e.g. van Adelsberg et al. Reference van Adelsberg, Lai, Potekhin and Arras2005; Potekhin et al. Reference Potekhin, Suleimanov, van Adelsberg and Werner2012; Potekhin Reference Potekhin2014). Such thin atmospheres were first used to model the spectrum of the strongly magnetic ( $10^{12}$ – $10^{13}\mbox{ G}$ ) neutron stars RX J0720.4 $-$ 3125 (Motch, Zavlin, & Haberl Reference Motch, Zavlin and Haberl2003) and RX J1856.5 $-$ 3754 (Ho Reference Ho2007; Ho et al. Reference Ho, Kaplan, Chang, van Adelsberg and Potekhin2007). Here we use nsmaxg models with condensed surface spectra from van Adelsberg et al. (Reference van Adelsberg, Lai, Potekhin and Arras2005). The best-fit results are shown in Table 3 and Fig. 8 (for a Thomson depth of 0.7, corresponding to an atmosphere thickness of 3.8 cm). We see that the best-fit is as good as the two blackbody fit and a better fit than that of a single atmosphere. Taking the cool component as being due to the entire neutron star surface, we estimate the distance to PSR J0837 $-$ 2454 to be $d\sim1\mbox{ kpc}$ . Results for somewhat larger depths than 0.7 are roughly within the uncertainties given in Table 3 and tend toward lower temperatures and a smaller distance, while smaller depths produce essentially blackbodies in X-rays.

We note here that the best-fit interstellar absorption $N_\textrm{H}$ is in the range $\sim$ $(1.8-2)\times10^{21}\mbox{ cm}^{-2}$ . These are a factor of 2–3 larger than the Galactic column density of $N_\textrm{H}=7.3\times10^{20}\mbox{ cm}^{-2}$ in the direction of PSR J0837 $-$ 2454 (HI4PI Collaboration et al. 2016). On the other hand, the empirical relation between $N_\textrm{H}$ and dispersion measure from He, Ng, & Kaspi (Reference He, Ng and Kaspi2013) yields $N_\textrm{H}=4\times10^{21}\mbox{ cm}^{-2}$ for the pulsar’s DM of $143.1\mbox{pc cm}^{-3}$ . The relatively large $N_\textrm{H}$ from our spectral fits, as well as the high DM, could be indicating significant material in our line-of-sight from, for example, a supernova remnant. This material may also contribute to the strong scattering seen in radio observations.