XMM-Newton and Chandra Observations of the Unidentified Fermi-LAT Source 3FGL J1016.5–6034: A Young Pulsar with a Nebula?

Figure 1 shows the XMM-Newton EPIC mosaic image of the J1016 field. The bright X-ray source is located toward the edge of the 3FGL 95% positional uncertainty ellipse of the Fermi-LAT source (shown in the left panel of Figure 1). In addition to a number of faint unresolved sources, there is evidence of faint large-scale emission located ∼5' south of the bright X-ray source (shown by the solid white circle in the right panel of Figure 1). Figure 2 shows the vicinity of CXOU J101546 as seen by Chandra. The subarcsecond angular resolution of the Chandra resolves the bright X-ray source into a point source surrounded by compact diffuse emission seen out to a distance of ∼30''. The radial profile of the extended emission is shown in Figure 3. The large-scale extended emission is too faint to be detected in the short Chandra exposure. A deeper Chandra observation is needed to study the properties of both the compact and large-scale extended emission.

Zoom In Zoom Out Reset image size

We have jointly fit the EPIC pn, MOS1, and MOS2 spectra of CXOU J101546, having a combined total of 5650 ± 140 net counts (the background contribution is ∼8%–11%). We first fit an absorbed power-law (PL) model to the data. The best-fit model has an absorbing hydrogen column density N_H = (2.2 ± 0.2) × 10²¹ cm⁻² and photon index ${\rm{\Gamma }}={2.29}_{-0.05}^{+0.06}$, with a reduced chi-squared ${\chi }_{\mathrm{red}}^{2}=1.17$ for ν = 229 degrees of freedom. In order to further improve the fit, we added an absorbed blackbody (BB) component to the PL model. This model (shown in Figure 4) has a best-fit absorbing hydrogen column density N_H = (2.4 ± 0.4) × 10²¹ cm⁻², temperature kT = 0.20 ± 0.02 keV, BB radius $R={430}_{-90}^{+170}{d}_{1\mathrm{kpc}}\ {\rm{m}}$, and photon index Γ = 1.8 ± 0.1 with a reduced chi-squared ${\chi }_{\mathrm{red}}^{2}=0.98$ for ν = 227 degrees of freedom.

Zoom In Zoom Out Reset image size

To evaluate the statistical significance of adding the additional BB component, we have used a likelihood ratio test,¹² with 30,000 simulated spectra. We find that the BB+PL model is preferred over the PL-only model at a >4σ level, implying that the BB component is statistically necessary to adequately fit the spectra. XMM-Newton lacks the angular resolution required to resolve the point source from the compact extended emission, and, therefore, this spectrum is a combination of both components. The best-fit model gives an absorbed source flux of F_X = (8.2 ± 0.2) × 10⁻¹³ erg cm⁻² s⁻¹ and unabsorbed thermal and PL component fluxes of ${2.4}_{-0.3}^{+0.4}\times {10}^{-13}$ erg cm⁻² s⁻¹ and (8.1 ± 0.3) × 10⁻¹³ erg cm⁻² s⁻¹, respectively, in the 0.5–10 keV band.

Chandra is able to resolve the point source (having 33 ± 6 net counts) from the surrounding extended emission (having 57 ± 8 net counts), but the short duration of the Chandra observation limits our ability to explore the spectra. Therefore, we chose to freeze N_H at the best-fit value obtained from the XMM-Newton spectral fits to reduce the number of free parameters. For the point source, the best-fit absorbed BB model¹³ has ${kT}={0.35}_{-0.05}^{+0.06}\ \mathrm{keV}$, $R={120}_{-30}^{+60}{d}_{1\mathrm{kpc}}\ {\rm{m}}$, and an absorbed flux of (1.6 ± 0.3) × 10⁻¹³ erg s⁻¹ cm⁻² (in 0.5–10 keV). The spectrum of the extended emission is well fit by an absorbed PL with Γ = 1.7 ± 0.3 and absorbed flux, ${F}_{X}=({5.1}_{-0.9}^{+1.0})\times {10}^{-13}$ erg s⁻¹ cm⁻² (in 0.5–10 keV). The unabsorbed flux of the PL component is $({6.0}_{-0.9}^{+1.0})\times {10}^{-13}$ erg s⁻¹ cm⁻² (in 0.5–10 keV). After fitting using W-stat, we calculate the Churazov weighted χ² values to assess the goodness of fit (Churazov et al. 1996). The weighted ${\chi }_{\nu }^{2}$ values are ${\chi }_{\nu }^{2}=1.01$ and ${\chi }_{\nu }^{2}=0.916$ for the point source and extended emission spectral fits, respectively, suggesting that the model fits the data reasonably well.

The spectrum of the large-scale extended emission seen in the XMM observation was extracted from the white solid circle, south of the point source, shown in Figure 1. We only use the pn data because the extended emission partially overlaps with one of MOS 1's damaged CCDs and MOS2 has fewer net counts than the pn. Following the approach taken by, for example, Younes et al. (2016), we first extracted the background spectrum from the region shown by the dashed white circle in Figure 1. Then we fit the background spectrum with a model containing a bremsstrahlung component and a PL component.¹⁴ We exclude energies between 1.4 and 1.6 keV to avoid complications related to fitting the narrow Al Kα instrumental line at ∼1.49 keV. We have checked that the exclusion of this line does not dramatically impact the resulting fit. We have also excluded energies above 7 keV because the extended emission becomes strongly background dominated. A total of 2056 ± 45 counts remain in the background spectrum after these energy cuts are applied. For the purposes of our extended emission analysis, we are not concerned with the physical interpretation of the background model as long as it fits the background spectrum reasonably well. The best-fit background model has a chi-squared χ² = 81 for ν = 73 degrees of freedom and adequately describes the background spectrum.

After fitting the background spectrum, the extended source spectrum (356 ± 46 net counts with the same energy cuts as above) is fit with a combination of an absorbed PL model and the best-fit background spectral model. In this fit, the best-fit background model normalizations are scaled to account for the differences in the region area and the exposure map between the background and source regions and frozen to these scaled values. The N_H for the source PL model was frozen at the best-fit value for the point source (i.e., N_H = 2.4 × 10²¹ cm⁻²). The best-fit PL spectrum of the extended emission has Γ = 2.1 ± 0.2 and an unabsorbed flux F = (1.6 ± 0.2) × 10⁻¹³ erg cm⁻² s⁻¹ in the 0.5–10 keV energy range with a chi-squared χ² = 52 for ν = 55 degrees of freedom. Following the same procedure, we have also fit the source spectrum with a bremsstrahlung model instead of a PL model and find a best-fit temperature ${kT}={2.4}_{-0.5}^{+0.7}\ \mathrm{keV}$ and unabsorbed flux F = (1.4 ± 0.2) × 10⁻¹³ erg cm⁻² s⁻¹ in the 0.5–10 keV energy range with a chi-squared χ² = 48 for ν = 55 degrees of freedom.

We have searched for periodicity using the ${Z}_{n}^{2}$ (for n = 1 and 2) test (Buccheri et al. 1983). In the 10⁻⁴–0.19 Hz range we used the combined XMM-Newton EPIC pn, MOS1, and MOS2 barycenter corrected data extracted from a circle (r = 665) centered on CXOU J101546, which contains 1258, 1329, and 4413 photons in MOS1, MOS2, and pn, respectively, in the 0.5–10 keV energy range. In addition, we used the pn data to search in the 0.19–6.8 Hz frequency range afforded by the higher time resolution of the pn detector. We have also performed the search in three energy bands—0.5–10 keV, 0.5–2 keV, and 2–10 keV. The most significant signal is ${Z}_{1}^{2}=32.57$ in the energy range of 2–10 keV found at a frequency of 5.871150(6) Hz using the pn-only data.¹⁵ This signal is significant only at the 2.6σ level,¹⁶ which is not enough to be considered a detection. Unfortunately, Chandra does not provide a high enough time resolution or number of counts to further test this periodicity. Observations with higher time resolution with EPIC pn operated in Small Window mode are needed to look for pulsations at higher frequencies typical for young pulsars.

We have also searched for X-ray variability in the source during our XMM observation by making a 500 s binned light curve. We find no evidence of variability in this light curve.

The 3FGL catalog,¹⁷ which is used to construct our model, was built from only 4 yr of Fermi data and was assembled using the Pass 7 IRF. Fermi has now been observing the sky for >10 yr , and the Pass 8 IRFs provide a number of dramatic improvements in the data (Atwood et al. 2013).

We use a standard binned likelihood analysis¹⁸ with a spatial binning of 005 and a spectral binning of four bins per decade of energy. Our analysis begins by including a search for new sources that may have appeared due to the increased sensitivity of the current LAT data. We find three new sources within 3° of the 3FGL position of J1016 having a test statistic (TS) larger than 40,¹⁹ the nearest of which is at a distance of ∼12. The new sources are all located within the Galactic plane (b > −2°) and could be caused by imperfections in the Galactic diffuse emission model (Acero et al. 2016). Therefore, we leave the detailed parameters of these new sources for the next FGL catalog. For the remainder of our analysis, we include these sources in our model.

The LAT spectrum of J1016 was best fit by an LP spectral model in the 3FGL catalog (Acero et al. 2015) due to J1016's significant spectral curvature (i.e., curvature_significance = 5.8). However, we have carried out the spectral analysis steps described below using both an LP and PL model and find TS values of ∼200 and ∼193, respectively. Following the definition of the curvature significance used for the 3FGL catalog (see Section 3.3 in Acero et al. 2015), we find that the source no longer has significant spectral curvature, possibly due to the additional exposure time, improvements in event reconstruction, and/or changes to the Galactic diffuse model. Therefore, we use the PL model for our analysis. To determine the position of J1016 more accurately we restrict the photon energies to 0.3–300 GeV (for the source localization only) because the LAT point-spread function becomes very broad at low energies and the source is relatively faint. The updated position (R.A. = 153999 and decl. = −60495) and corresponding 2σ error ellipse (${r}_{\mathrm{semimajor}}=0\buildrel{\circ}\over{.} 090$, ${r}_{\mathrm{semiminor}}=0\buildrel{\circ}\over{.} 069$, θ_ellipse = 17258) are shown in the left panel of Figure 1. Following the 3FGL catalog approach, we account for systematic uncertainties by multiplying both ellipse axes by a factor of 1.05 and adding an additional 0005 in quadrature to each 95% ellipse axis, leading to 95% statistical plus systematic uncertainties of ${r}_{\mathrm{semimajor}}=0\buildrel{\circ}\over{.} 095$, ${r}_{\mathrm{semiminor}}=0\buildrel{\circ}\over{.} 072$ (Acero et al. 2015). The error ellipse size has shrunk substantially compared with that from 3FGL catalog, and the position of CXOU J101546 is consistent with the 2σ positional uncertainty. For completeness, we note that the LP fit gives a different position (R.A. = 154.099, decl. = −60.487) that is offset by ∼3' from the PL fit position, but the two error ellipses still overlap considerably and both of them include CXOU J101546 (see Figure 1).

Next we fit the PL spectral model, fixing the source positions, while freeing all components of the Galactic and isotropic diffuse emission, as well as the normalizations and spectral parameters of any point sources located within 3° of the refined position of J1016. We then produced and inspected the TS and residual maps resulting from the best-fit PL model to ensure that there was no significant excess emission left in the vicinity of J1016. The best-fit PL model for J1016 has Γ = 2.85 ± 0.07 and GeV flux F = (1.8 ± 0.1) × 10⁻¹¹ erg s⁻¹ cm⁻² (in 0.1–300 GeV) with a TS = 192.73. Lastly, we produced the spectral energy distribution (SED) of J1016, which is calculated by fitting the flux normalization in a number of energy bins using a fixed-index PL spectral model (the photon index is fixed at the best-fit value Γ = 2.85). The SED has a flux F = (3.2 ± 0.7) × 10⁻¹¹ erg s⁻¹ cm⁻² in 0.1–100 GeV (see Figure 5).

Zoom In Zoom Out Reset image size

For the timing analysis we used the same filtering as described in Section 2.3 but restricted the events to be within an r < 07 radius from the position of CXOU J101546. We have also applied a barycentric time correction to all LAT events using the spacecraft orbit file with the gtbary tool. We then employed the approach described in Pletsch & Clark (2014) combined with the weighting scheme suggested by Bruel (2019) and tested by Smith et al. (2019) to perform a blind search for pulsations in the frequency and frequency derivative domain (f = 0–500 Hz and $\dot{f}=-{10}^{-9}\mbox{--}0\,{\rm{s}}$ s⁻¹). In the Pletsch & Clark (2014) algorithm we used a window size of T_coh = 6 days. Our search did not reveal any significant pulsations. We also checked the X-ray periodicity candidate frequency of 5.871150(6) Hz using a modified ${Z}_{n}^{2}$ (for n = 1 and 2) including a modification to account for the frequency derivative but did not see any substantial periodic signal in the LAT data at this frequency.

To search for variability, we have also constructed the Fermi-LAT light curves with weekly and monthly time bins. By fitting a straight line to the data, we find that the source is variable at the 1.5σ level, which is not significant. Therefore, we conclude that the source is not variable at GeV energies.