Hubble Space Telescope Detection of the Millisecond Pulsar J2124−3358 and its Far-ultraviolet Bow Shock Nebula

To identify the pulsar in the HST images (see Figure 1), we had to check the astrometry of our data. We used stars from the Guide Star Catalog (GSC v.2.3) within the field of view of the WFC3/UVIS to improve the astrometry of the F475X image. Other catalogs with more precise coordinate measurements (e.g., 2MASS, UCAC4) had too few objects in the UVIS field of view. Moreover, most of these objects are bright stars, saturated in the UVIS image, which prevents accurate centroiding and renders them unsuitable for the astrometric correction. Before the correction, we found the systematic offset between the GSC v.2.3 and the UVIS stars to be ≈04 (consistent with the typical HST pointing error⁹ ). We then selected 20 stars from the GSC with obvious UVIS counterparts suitable for accurate centroid determination. As these stars did not have proper motion measurements, we could not adjust their coordinates for the epoch of the HST observations. This can explain four outliers (i.e., UVIS stars having ≳3 times the mean offset from the catalog position), which we interpret as a proper motion effect. These stars were removed from our final sample of reference stars. We then used the IRAF task DAOFIND to find the detector coordinates of the centroids for the remaining 16 reference stars and supplied them to the IRAF task CCMAP ¹⁰ to calculate new astrometric solution by matching these stars with their GSC counterparts. The refined astrometric solution has standard deviations σ_α = 017 and σ_δ = 016 for the sample of 16 stars (at 68% confidence).

In the SBC image only five sources are seen. Since none of them have counterparts in catalogs, we had to align the F125LP and F140LP images with the UVIS image. Three of the sources, seen in both the UVIS and SBC images, were used in the alignment process. However, all three sources appear to be extended in both images (they are most likely background galaxies). This hinders accurate centroid determination and leads to an uncertainty¹¹ of ≈01 in aligning the SBC images with the UVIS image.

In the astrometry-corrected UVIS image we found a faint source located at α = 21:24:43.841(16) and δ = −33:58:45.01(18). These coordinates are offset by ${\rm{\Delta }}\alpha \cos \delta =-0\buildrel{\prime\prime}\over{.} 01\pm 0\buildrel{\prime\prime}\over{.} 19$ and Δδ = −018 ± 018 from the pulsar's radio coordinates (α = 21:24:43.841662(24) and δ = −33:58:45.1897(5)) expected at the epoch of the UVIS observation (MJD 56976). Since the offset is within the alignment uncertainty, we conclude that this is a very viable candidate for the pulsar counterpart.¹²

To find the optimal aperture for the photometry in SBC/F125LP image, we calculated the signal-to-noise ratio (S/N) for a set of varying circular apertures centered on the source. The background was estimated from an annular region with r_in = 20 and r_out = 50 pixels (064 and 16, respectively) centered on the source. We found that the r = 5 pixels (016) aperture provides a maximum S/N ≈ 9. This aperture corresponds to about 64% of the encircled energy, according to the ACS Instrument Handbook.¹³ Repeating the process for the F140LP image, we found an optimal aperture (S/N ≈ 3) of r = 4 pixels (013) corresponding to the encircled energy fraction of 58%. The pulsar detection is less significant because the exposure time of the F140LP observations is a factor of 6 shorter and because the throughput is reduced compared to F125LP. Note that the pulsar photometry was obtained on individual images to avoid additional broadening of the source caused by the imperfect alignment.

We repeated this procedure for the UVIS/F475X image, where the background was taken from an annulus with r_in = 8 and r_out = 12 pixels (032 and 048, respectively). We found the optimal extraction aperture of r = 3 pixels (012), which corresponds to the encircled energy fraction of ≈78% and provides an S/N ≈ 19.

The net source count rate was calculated as ${C}_{s}={C}_{t}\ -{C}_{b}({A}_{s}/{A}_{b})$, where C_t is the total count rate in the source region of area A_s, and C_b is the background count rate in the region of area A_b. Correspondingly, the source count rate uncertainty was calculated as $\delta {C}_{s}={[{C}_{t}+{C}_{b}{({A}_{s}/{A}_{b})}^{2}]}^{1/2}\,{t}_{\exp }^{-1/2}$, where t_exp is the exposure time. The net count rates for the optimal apertures are given in Table 2.

We then corrected the rates for the finite aperture sizes and followed standard photometry procedures to convert count rates to flux densities.¹⁴ The observed mean flux densities for F125LP, F140LP, and F475X are 27 ± 4 nJy, 25 ± 11 nJy, and 22 ± 2 nJy, respectively. For a plausible reddening E(B − V) = 0.03 (see Section 2.3) these values correspond to dereddened flux densities of 35 ± 6 nJy, 33 ± 15 nJy, and 24 ± 2 nJy, for F125LP, F140LP, and F475X, respectively.

To interpret the photometry results, we need to know the interstellar extinction toward J2124. The extinction coefficient A_ν is proportional to the color excess E(B − V), which, on average, is proportional to the hydrogen column density N_H. Zavlin (2006) reported N_H = (1–3) × 10²⁰ cm⁻² toward J2124 based on fits to the X-ray data. The dispersion measure, DM = 4.6 pc cm⁻³, leads to a similar N_H = 1.4 × 10²⁰ cm⁻² (for an assumed 10% degree of ISM ionization). The relation ${N}_{{\rm{H}}}=2.2\times {10}^{21}{R}_{V}E(B-V)$ cm⁻² (Gorenstein 1975; Güver & Özel 2009), with the commonly assumed R_V = 3.1, corresponds to E(B − V) = 0.015–0.045. In our analysis below, we explore an even broader range of color excess, E(B − V) = 0.01–0.08, to account for the additional uncertainty introduced by the N_H-to-E(B − V) conversion. The upper bound on color excess, E(B − V) = 0.08, in the direction to J2124 (l = 109, b = −454) is based on the Galactic extinction maps by Schlegel et al. (1998). To perform dereddening, we use the extinction curves adopted from Clayton et al. (2003).

The photometric data in only three spectral bands (one of which, F140LP, is within the other, F125LP) can be fit with many spectral models. We first fit the data with an absorbed PL model, ${f}_{\nu }={f}_{{\nu }_{0}}{(\nu /{\nu }_{0})}^{\alpha }\times {10}^{-0.4{A}_{\nu }}$, which is often used to describe magnetospheric emission (Pavlov et al. 1997). We fit the PL slope α and normalization ${f}_{{\nu }_{0}}$ at ${\nu }_{0}=6.07\times {10}^{14}$ Hz (corresponding to the pivot wavelength of the F475X filter) for four fixed values of color excess: E(B − V) = [0.01, 0.03, 0.05, 0.08]. An example of such a fit is shown in Figure 4, and the PL parameters for the four E(B − V) values are given in Figure 5. The slope and normalization values on the E(B − V) grid are α = [0.20 ± 0.02, 0.28 ± 0.02, 0.35 ± 0.02, 0.46 ± 0.02] and ${f}_{{\nu }_{0}}$ = [22.6 ± 0.2, 24.1 ± 0.3, 25.7 ± 0.3, 28.4 ± 0.3] nJy (at 6.07 × 10¹⁴ Hz), respectively. Their dependences on color excess can be approximated by linear functions in this E(B − V) range: α ≈ 0.164 + 3.71 E(B − V) and ${f}_{{\nu }_{0}}\approx 21.6+82.9\,E(B-V)$ nJy. The corresponding dereddened fluxes in the FUV and optical bands are F(1250–2000 Å) = [2.6 ± 0.3, 2.9 ± 0.3, 3.4 ± 0.4, 4.4 ± 0.5] × 10⁻¹⁶ and F(3800–6000 Å) = [6.6 ± 0.4, 7.1 ±0.4, 7.6 ± 0.4, 8.4 ± 0.5] × 10⁻¹⁷ erg s⁻¹ cm⁻².

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We then considered a single thermal (BB) model with a fixed radius R = 12 km, a typical radius of an NS as seen by a distant observer. We fit both the BB temperature T and color excess E(B − V) for distances d = 340–500 pc. Because the FUV band falls in the Rayleigh–Jeans regime for the obtained temperatures, we can express the final results as $T=(3.6\pm 0.2)\times {10}^{6}{d}_{410}^{2}$ K and E(B − V) = 0.51 ± 0.01. Since this color excess is much higher than the Galactic value along the line of sight and the temperature is implausibly high for the entire NS surface (contradicts the X-ray data), this model can be ruled out.

Last, we fit a two-component BB+PL model. Given the small number of data points, we fit the BB temperature T and PL normalization ${f}_{{\nu }_{0}}$ on a 3D grid α = [−1, 0], d = [340, 410, 500] pc, and E(B − V) = [0.01, 0.03, 0.05, 0.08]. The resulting BB temperatures and PL normalizations are plotted in Figure 6. Examples of best-fit model spectra at most plausible d = 410 pc, E(B − V) = 0.03 are shown in Figure 7 for two values of PL slope, α = 0 and −1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A bow-shaped structure spatially coincident with the known Hα bow shock (see Gaensler et al. 2002; Brownsberger & Romani 2014) is clearly visible in the F125LP image shown in Figure 1 (top left). The FUV bow shock appears to be asymmetric, being thicker and brighter east of the pulsar. We measure the apex distance (from the pulsar to the leading edge of the bow shock) to be ≈25, consistent with the measurements in Hα by Brownsberger & Romani (2014), Gaensler et al. (2002).

To extract the bow-shock flux, we defined the source and two background regions, "inside" and "outside" of the bow shock, shown in Figure 1 (top right). The background regions were selected based on the F125LP image, in which the bow shock was most clearly detected. The regions are close to the bow-shock extraction region to mitigate effects from the inhomogeneous detector background while excluding the pulsar. Because of the larger apparent thickness of the east side of the bow shock, more space was left between the bow-shock extraction region and the background regions. We took the background variations into account by selecting background regions on both sides of the bow shock: one "inside" and one "outside" the bow shock. These regions were used for flux measurements in both F125LP and F140LP. A combined background from the two regions was extracted. The net source count rate was calculated similarly to the pulsar, with the area A_b being the combined area of the "outside" and "inside" background regions. The F125LP and F140LP count rates and their uncertainties are given in Table 2. From the measured count rates we calculated the absorbed fluxes in the F125LP and F140LP wavelength ranges: F(1250–2000 Å) = (8.1 ± 0.4) × 10⁻¹⁵ erg s⁻¹ cm⁻² and F(1350–2000 Å) = (7.0 ± 0.7) ×  10⁻¹⁵ erg s⁻¹ cm⁻², for a flat f_λ spectrum. The unabsorbed luminosity can be estimated as L(1250–2000 Å) ≈ 1.9 × 10²⁹ erg s⁻¹, for E(B − V) = 0.03, and d = 410 pc.

The bow shock is not detected in the F475X images (see Figure 8). Because of the apparent non-uniformities in the background, which significantly exceed the statistical fluctuations, we chose to measure the upper limits by sampling the background counts from ten different regions in the vicinity of the bow shock and calculating the standard deviation from the mean value. The areas of the background regions were equal to that of the bow-shock region used for FUV flux measurements. The standard deviation can be considered as a conservative 1σ upper limit. The corresponding 3σ upper limit on the bow-shock flux in the area of 52 arcsec² is F(3750–6000 Å) ≈ 2.9 × 10⁻¹⁵ erg s⁻¹ cm⁻², for a flat f_λ spectrum and E(B − V) = 0.03.

Zoom In Zoom Out Reset image size

Similar to our analysis of the FUV bow shock around J0437 (Rangelov et al. 2016), we measured the ratio of the count rates in F125LP and F140LP filters for the J2124 bow shock. We obtained a ratio of C₁₂₅/C₁₄₀ = 1.3 ± 0.2. This is lower than C₁₂₅/C₁₄₀ = 1.63 ± 0.08 for the J0437 bow shock, but the associated uncertainty is large.

The obtained photometric measurements can be used to constrain the optical-FUV spectral energy distribution (SED) for J2124 and to look for thermal emission, similar to that of J0437.

The slope α ≳ 0.2 inferred from the nonthermal (single PL) fit is outside the typical range of −1 ≲ α ≲ 0 found for the optical-UV magnetospheric emission seen in other pulsars (Kargaltsev & Pavlov 2007; Mignani 2011). However, we cannot rule out the possibility that the FUV-optical emission of J2124 is entirely magnetospheric, with α > 0, because virtually nothing is known about optical spectra of non-accreting recycled pulsars.¹⁵ The PL fit implies the V-band luminosity ${L}_{{\rm{V}}}\approx 5\times {10}^{26}{d}_{410}^{2}$ erg s⁻¹, which corresponds to an optical efficiency¹⁶ ${\eta }_{V}\equiv {L}_{V}/{\dot{E}}_{\mathrm{int}}\sim 2\times {10}^{-7}$, similar to the efficiencies of middle-aged pulsars (see Figure 4 in Danilenko et al. 2012).

A thermal interpretation for combined FUV and optical emission (single BB fit) can be ruled out because it requires a color excess much higher than the Galactic value along this line of sight, as well as an implausibly high temperature, ∼3.6 × 10⁶ K, of the entire NS surface, which would result in an X-ray luminosity much higher than observed.

Finally, there remains a good possibility that the optical emission is predominantly nonthermal while the FUV emission is largely thermal, similar to what is observed for nearby middle-aged pulsars (Kargaltsev & Pavlov 2007). The optical efficiency in the PL+BB fit is similar to (slightly lower than) that of the single PL fit, because the thermal component contribution is very small in the optical (see Figure 7). It follows from Figure 6 that the NS surface temperature of J2124 is most likely in the range T ≈ (0.5–2.1) × 10⁵ K for the considered (conservative) range of color excess values, E(B − V) = 0.01–0.08, and distances d = 340–500 pc. Thus this interpetation of the J2124 spectrum yields the NS surface temperature similar to that of J0437 (Durant et al. 2012).

Zavlin (2006) has shown that two components are required to fit the X-ray spectrum of J2124. However, various combinations of thermal and nonthermal components were found to be statistically acceptable (e.g., two-component H-atmosphere model, H-atmosphere + PL, BB+PL; see Section 1). Extrapolation of the thermal components, originating from small hot areas on the NS surface (presumably polar caps), underpredicts the observed FUV fluxes by several orders of magnitude. The X-ray-emitting heated regions therefore cannot be responsible for the observed FUV emission from the J2124. The extrapolated best-fit X-ray PL components of the BB+PL and H-atmosphere + PL models overpredict the FUV-optical fluxes by a large margin, but this is not unusual for middle-aged pulsars, which often show a spectral break between the X-ray and optical bands. The extrapolated γ-ray spectrum is a few orders of magnitude below the optical-FUV fluxes, which means that the optical-FUV and γ-ray emission are produced by different mechanisms.

If the two-component (thermal plus magnetospheric) nature of the optical-FUV spectrum of J2124 is confirmed in future observations, then J2124 would be the second recycled pulsar, after J0437, with such a high temperature, T ∼ 10⁵ K of the NS surface. It would firmly prove that a heating mechanism operates at least in millisecond pulsars and could distinguish between heating models.

Heating mechanisms that may operate in the interiors of ancient NSs include the dissipation of rotational energy due to interactions between superfluid and normal components of the NS (vortex creep; Alpar et al. 1984; Shibazaki & Lamb 1989), release of strain energy stored by the solid crust due to spin-down deformation (Cheng et al. 1992), rotochemical heating (Reisenegger 1995; Fernández & Reisenegger 2005; Petrovich & Reisenegger 2010), rotation-induced deep crustal heating (a variant of rotochemical heating that operates in the crusts of pulsars; Gusakov et al. 2015), and decay or annihilation of dark matter particles, trapped in the NS interiors (de Lavallaz & Fairbairn 2010; Kouvaris & Tinyakov 2010). Of the internal heating mechanisms, the vortex creep, rotochemical heating, and deep crustal heating, all of which are driven by the spin-down of pulsars, are able to account for the high surface temperatures, ∼10⁵ K, of millisecond pulsars (Gonzalez & Reisenegger 2010; Gusakov et al. 2015). Both J2124 and J0437 are much older than the standard NS cooling time, therefore they are expected to be in a quasi-stationary state in which any internal heating and photon cooling balance each other (e.g., Reisenegger 1995; Fernández & Reisenegger 2005). Thus the effective temperature should be a function of the spin-down parameters, $\propto {\dot{{\rm{\Omega }}}}^{1/4}$ for vortex creep, $\propto {({\rm{\Omega }}\dot{{\rm{\Omega }}})}^{2/7}$ for rotochemical heating (assuming modified Urca reactions and neglecting superfluid corrections), and $\propto {({\rm{\Omega }}\dot{{\rm{\Omega }}})}^{1/4}$ for rotation-induced deep crustal heating, where Ω is the NS rotation rate and $\dot{{\rm{\Omega }}}$ is its time derivative. The spin-down parameters (corrected for the Shklovskii effect) of J2124 are quite similar to those of J0437, therefore the ratio of their effective temperatures is expected to be T_J2124/T_J0437 = 0.89, 0.95, 0.96, respectively, for the three heating mechanisms. Since for a radius of 12 km (as assumed here), T_J0437 ≈ 1.8 × 10⁵ K (inferred from the results of Durant et al. 2012), the temperature expected for J2124 is only slightly lower, T_J2124 ∼ (1.6–1.7) × 10⁵ K, which is consistent with temperatures estimated in the PL+BB scenario. To confirm this scenario and estimate the temperature more accurately, the J2124 pulsar must be observed with several filters in the UVOIR range, and the distance to the pulsar should be measured more precisely.

The pulsar J2124 is the second pulsar with a bow shock detected in FUV. The first FUV bow shock was discovered by Rangelov et al. (2016) around J0437. These authors tested two spectral models for FUV emission. A PL continuum model with −3 ≲ α ≲ −1 in the FUV range (α ≳ −1 if the PL spectrum extends to the optical range) could correspond to synchrotron radiation of relativistic electrons leaked from the X-ray PWN region and trapped at the forward bow shock. The other model, favored by Rangelov et al. (2016), was the spectrum emitted by the ISM matter heated and compressed in a collisionless shock (SHELLS model in Bykov et al. 2013); it was consistent with the measured bow-shock fluxes at reasonable values of the shock (pulsar) velocity and upstream ISM density.

The bow shock of J2124 is dimmer and has a lower S/N, which prevents us from performing a similarly detailed analysis. Nevertheless, we can compare general properties of the J2124 shock with those of the J0437 shock. The FUV bow shocks of both pulsars are spatially coincident with the Hα bow shocks (see Figure 1). In contrast to J0437, the J2124 bow shock shows a significant asymmetry of its shape with respect to the proper motion direction, which has been attributed to a local ISM density gradient perpendicular to the pulsar's direction of motion (Gaensler et al. 2002). Interestingly, there is also an asymmetry in shock brightness, opposite in the the Hα and FUV images (compare the top left and bottom right panels in Figure 1); the reason for this asymmetry remains unclear. In addition to the asymmetry, the head (apex region) of the FUV-Hα J2124 bow shock is flatter than the head of the J0437 bow shock. This could be caused by anisotropy of the pulsar wind, such that the wind is concentrated in the pulsar equatorial plane that is nearly perpendicular to the pulsar velocity (Gaensler et al. 2002; Brownsberger & Romani 2014).

The FUV luminosity, L(1250–2000 Å) ≈ 1.9 ×10²⁹${d}_{410}^{2}$ erg s⁻¹, extracted from the A_s = 52 arcsec² area of the J2124 bow shock, is a factor of $\approx 3.8{d}_{410}^{2}$ higher than the FUV luminosity from the A_s = 32 arcsec² area of the J0437 shock. However, the (distance-independent) mean specific intensity, $I=L/(4\pi {d}^{2}{A}_{s})\approx 1.9\times {10}^{-16}$ erg cm⁻² s⁻¹ arcsec⁻² for the J2124 bow shock is a factor of 2.8 lower than that of the J0437 shock. This difference could be partly due to the fact that the width of the extraction region of the J2124 shock was somewhat broader than the actual width of the bright shock area. In addition, the more distant J2124 shock was imaged at larger physical separations in units of length from the apex, where the shock becomes intrinsically dimmer. It is therefore most likely that the specific intensities of the J2124 and J0437 shocks are close to each other at the same (physical) separations from the apices. Since in the ISM shock interpretation the emitted shock flux per unit emitting area is crudely proportional to the upstream density and only weakly depends on the pulsar velocity in the relevant range of parameters (see Figure 6 of Rangelov et al. 2016), which means that the ISM densities ahead of the J0437 and J2124 shocks are close to each other.

For pulsar J0437, Rangelov et al. (2016) used the ratio of the F125LP and F140LP count rates, C_F125/C_F140 = 1.63 ± 0.08, to estimate the slope of the J0437 spectrum in the FUV range, −3 ≲ α ≲ −1 (see Figure 4 in Rangelov et al. 2016). Similarly, for J2124, we have C_F125/C_F140 = 1.3 ± 0.2, from which we infer −10 ≲ α ≲ −4, i.e., a very steeply decreasing f_ν spectrum. This result suffers from the large uncertainty of the F140LP count rate, however. Such a steep spectrum cannot be extrapolated to the optical (F475X) band because the extrapolated spectral flux would greatly exceed the upper limit in that filter, which requires a FUV-optical slope α ≳ 0 (see Figure 9), corresponding to C_F125/C_F140 ≳ 1.85. Unless the count rate uncertainties are greatly underestimated as a result of unaccounted systematic errors, this discrepancy can be considered as an argument for a shock SED more complex than a simple PL.

Zoom In Zoom Out Reset image size

In the shocked ISM emission model the ratio of FUV and Hα fluxes depends on the model parameters. Brownsberger & Romani (2014) report the Hα photon flux ${{ \mathcal F }}_{{\rm{H}}\alpha }\ =5.3\times {10}^{-3}$ photons s⁻¹ cm⁻², or energy flux F_Hα = 1.6 ×10⁻¹⁵ erg s⁻¹ cm⁻². This corresponds to the ratio of F125LP to Hα energy fluxes of ≈6, although this value can be higher (by a factor of up to 2) because we do not know the exact apex areas used by Brownsberger & Romani (2014). For comparison, the flux ratio in the same filters for J0437 is F_FUV/F_Hα = 12 ± 4 (Figure 6 in Rangelov et al. 2016), i.e., generally consistent with J2124. The non-detection of the J2124 bow shock in the F475X image provides additional constraints on the spectral shape and possibly on the emission mechanism. Figure 9 shows a single PL model and a shocked ISM model (normalized to FUV flux measurement) compatible with the 3σ F475X upper limit. The shocked ISM model appears to be able to describe the measurements; however, given the complexity of the model and the scarcity of the data, we do not attempt to fit it.

The detection of the FUV bow shock around two pulsars allows us to assume that such shocks can be detected from many other pulsars, including those from which no Hα shock is seen because of a lack of neutral hydrogen atoms upstream of the shock. Such observations, currently possible only with the HST, would be very useful for studying the ISM properties and the physics of relativistic shocks.