Observing and Modeling the Gamma-Ray Emission from Pulsar/Pulsar Wind Nebula Complex PSR J0205+6449/3C 58

3C 58 is an extended flat-spectrum radio source that was identified as a supernova remnant (SNR) first in radio (Weiler & Seielstad 1971) and then in optical by Hα observations (van den Bergh 1978). Because of its flat radio spectrum and bright-filled center, 3C 58 was classified as a pulsar wind nebula (PWN, or plerion, Weiler & Panagia 1978) long before the central pulsar, PSR J0205+6449, was discovered (Murray et al. 2002). Subsequent radio imaging observations continued to reveal a center-brightened morphology and a compact size of 6' × 9' (Green 1986; Reynolds & Aller 1988; Bietenholz et al. 2001; Bietenholz 2006), which is consistent with the morphology observed in infrared and X-ray bands (Bocchino et al. 2001; Slane et al. 2004, 2008). Due to the jet-torus structure (Slane et al. 2004), filaments and knots (Fesen et al. 2008) observed in 3C 58 resemble those seen in the Crab Nebula, 3C 58 was proposed to be a "Crab-like" PWN. The distance of 3C 58 was estimated to be 3.2 kpc (Roberts et al. 1993), though, there is an on-going discussion regarding this value, given that a recent H i measurement suggests a nearer distance of just 2 kpc (Kothes 2013). 3C 58 was proposed to be the remnant of SN 1181, observed by Chinese and Japanese astronomers (Stephenson 1971; Stephenson & Green 2002). However, there is controversy regarding the viability of this connection (see Fesen et al. 2008, Table 3 and Kothes 2013, Table 1 for a discussion).

PSR J0205+6449 is a pulsar with a rotational period of 65 ms located in 3C 58. It was discovered by Chandra in X-rays many years after 3C 58 was classified as a PWN. Timing parameters give a surface magnetic field estimation of 3.6 × 10¹² G, a characteristic age of 5400 years, and a very high spindown luminosity, 2.7 × 10³⁷ erg s⁻¹, making it the third most energetic of the known Galactic pulsars (Murray et al. 2002). Camilo et al. (2002) reported the detection of its radio pulsation, which leads the X-ray pulse by an∼0.10 spin phase. Because of the high spindown power of PSR J0205+6449, it was expected to shine in gamma-rays. PSR J0205+6449 was among the first gamma-ray pulsars detected by Fermi-LAT (Abdo et al. 2009). Gamma-ray emission from the PWN 3C 58 was reported in the second Fermi-LAT Pulsar Catalog (Abdo et al. 2011, which we shall refer to as 2PC hereafter). 3C 58 was also detected above 10 GeV in the First and Third Catalog of Hard Fermi-LAT Sources, suggesting its potential nature as a TeV gamma-ray source (Ackermann et al. 2013, 1FHL; Ajello et al. 2017, 3FHL). Several imaging atmospheric Cherenkov telescopes have observed 3C 58 (Whipple, Hall et al. 2001; VERITAS, Aliu 2008; MAGIC I, Anderhub et al. 2010) but it was only recently detected by MAGIC II (Aleksić et al. 2014a). PWN models for 3C 58 have been presented by Bednarek & Bartosik (2003, 2005), Bucciantini et al. (2011), and Torres et al. (2013); in the latter paper, a comparison among these models is provided. Here, in an effort to better understand the radiation from the complex, we report on further analysis/modeling of PSR J0205+6449/3C 58 using more than eight years of Fermi-LAT data and the newest response functions.

The Fermi-LAT (Atwood et al. 2009) data used for this paper range from 2008 August 4 (MJD 54682) to 2016 September 20 (MJD 57651), covering 97 months and extending the three years' coverage of the 2PC. The Fermi Science Tools⁶ 10-00-05 release was used to analyze the data. The data selection and analysis method adopted in this paper are similar to those in Li et al. (2017). We selected events from the "Pass 8" event class, and used "P8R2 V6 Source" instrument response functions (IRFs). All gamma-ray photons within an energy range of 0.1–300 GeV and a circular region of interest (ROI) of 10° radius centered on PSR J0205+6449 were considered. Gamma-ray photons were selected only with a zenith angle <90° to reject contaminating gamma-rays from the Earth's limb.

The spectral results presented in this work were calculated by performing a binned maximum likelihood fit (30 logarithmically spaced bins in the 0.1–300 GeV range) using the Science Tool gtlike. A spectral-spatial model was constructed including Galactic and isotropic diffuse emission components ("gll_iem_v06.fits," Acero et al. (2016), and "iso_P8R2_CLEAN_V6_v06.txt," respectively⁷ ) as well as known gamma-ray sources within 15° of the 3C 58, based on the Fermi-LAT Third Source Catalog (3FGL, Acero et al. 2015). The spectral parameters of the sources within 3° of our target were left free. The spectral parameters of other sources included were fixed at the 3FGL values. In the pulsar phase-resolved analysis, photons within a specific spin phase interval are selected. To account for that, the prefactor parameters of all sources were scaled to the width of the spin phase interval. The significance of the sources is evaluated by the test statistic (TS). This statistic is defined as TS = −2 $\mathrm{ln}$ (L_max,0/L_max,1), where L_max,0 is the maximum likelihood value for a model in which the source studied is removed (the "null hypothesis"), and L_max,1 is the corresponding maximum likelihood value with this source incorporated. The larger the value of TS, the less likely the null hypothesis is correct (i.e., a significant gamma-ray excess lies on the tested position). The square root of the TS is approximately equal to the detection significance of a given source. The pointlike analysis package (Kerr 2011a) was used to produce the TS maps in this paper. In the analysis, the systematic errors have been estimated by repeating the analysis using modified IRFs that bracket the effective area⁸ (Ackermann et al. 2012), and artificially changing the normalization of the Galactic diffuse model by ±6% (2PC). The latter dominates the systematic errors. Energy dispersion correction has been adopted in spectral analysis.⁹ In this paper, the first (second) uncertainty shown corresponds to the statistical (systematic) error.

To search for the possible spatial extension of 3C 58 in the off-peak gamma-ray emission, we followed the method described in Lande et al. (2012). The source is modeled to be spatially extended with a symmetric disk model. We fitted its position and extension with the pointlike analysis package. The extension significance was defined as TS_ext = −2($\mathrm{ln}$L_point/L_disk), in which L_disk and L_point were the gtlike global likelihood of the extended source hypotheses and the point source, respectively. A threshold for claiming the source to be spatially extended is set as TS_ext > 16, which corresponds to a significance of ∼4σ.

Photons from PSR J0205+6449 within a radius of 065 and a minimum energy of 200 MeV were selected, which maximized the H-test statistic (de Jager et al. 1989; de Jager & Büsching 2010). Adopting the most updated ephemeris for PSR J0205+6449 (M. Kerr & D. Smith 2017, private communication), we assigned pulsar rotational phases to each gamma-ray photon that passed the selection criteria, using Tempo2 (Hobbs et al. 2006) with the Fermi plug-in (Ray et al. 2011).

The light curve of PSR J0205+6449 was divided into two parts, an off-peak and an on-peak interval. We began by deconstructing the pulsed light curve into simple Bayesian Blocks using the same algorithm described in the 2PC, by Jackson et al. (2005) and Scargle et al. (2013). To produce Bayesian Blocks on the light curve, we have extended the data over three rotations, by copying and shifting the observed phases to cover the phase range from −1 to 2. We have defined the final blocks to be between phases 0 and 1. To increase the statistics and in accordance with the current results, we have adopted a wider interval for the off-peak phases than that used in the 2PC (0.35 of the total). The off-peak interval in our analysis is then defined to be at ϕ = 0.0–0.184, 0.291–0.574, and 0.786–1.0, yielding 0.681 of the the total revolution. We also tested a conservative selection for the off-peak phases, selecting them as ϕ = 0.0–0.144 and 0.825–1.0, which is defined as the lowest Bayesian block with 10% reduction on either side (2PC), yielding 0.319 of the total revolution. It leads to consistent results. The on-peak phases are thus located at ϕ = 0.184–0.291 and ϕ = 0.574–0.786. Figure 1 shows the pulsar spin light curves, using a photon weighting technique based on the method of Kerr (2011b). An additional discussion of Figure 1 is presented in Section 5.

Zoom In Zoom Out Reset image size

The off-peak emission of PSR J0205+64493C 58 was previously reported in the 2PC using three years of data with P7V6 IRFs (P7V6 is one of the previous versions of LAT IRFs¹⁰ ), yielding a TS value of 33.7 and a flux level of (1.75 ± 0.68)×10⁻¹¹ erg cm⁻² s⁻¹ in the 0.1–316 GeV range. This emission was thought to be the pulsar wind nebula radiation. 3FGL was based on four years of Pass 7 LAT data while in this analysis the data set is doubled (eight years of Pass 8 data). Thus the region around PSR J0205+6449/3C 58 may not be well modeled solely by 3FGL sources. For a better spectral-spatial modeling of this region, we added three additional point sources (assumed to be described by a simple power law) to the spatial model (Figure 2), using a similar method as that of Caliandro et al. (2013). The best positions of the additional sources were determined with pointlike as R.A. = 3119 ± 0.05, decl. = 6691 ± 0.04 (NEW 1); R.A. = 3359 ± 0.05, decl. = 6252 ± 0.04 (NEW 2); R.A. = 3530 ± 0.04, decl. = 6292 ± 0.04 (NEW 3). Assuming a power-law spectral shape (dN/dE =N₀(E/E₀)^−Γ cm⁻² s⁻¹ GeV⁻¹), the off-peak analysis of PSR J0205+6449/3C 58 results in a TS value of 202 and a flux level of (1.24 ± 0.13 ± 0.17)×10⁻¹¹ erg cm⁻² s⁻¹ in the 0.1–300 GeV energy range, which is consistent with 2PC results. The spectral index is estimated as 2.04 ± 0.07 ± 0.15 (Table 1), which is softer than the value reported in the 2PC (1.61 ± 0.21) but still within 3σ errors. Figure 2 (left panel) shows the TS map of the PSR J0205+6449/3C 58 region during off-peak phases. In off-peak phases, we also modeled PSR J0205+6449/3C 58 with a power law having an exponential cutoff (dN/dE = N₀(E/E₀)^−Γ exp(−E/E₀) cm⁻² s⁻¹ GeV⁻¹). The two models are compared using the likelihood ratio test (Mattox et al. 1996). The ΔTS¹¹ between the two models is 0.03, which indicates that a cutoff is not significantly preferred. This result is also consistent with the 2PC. The best-fit spectral parameters and corresponding TS values are listed in Table 1, while the spectral energy distributions (SEDs)¹² along with the best-fit power-law model are shown in Figure 3 (left panel).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The extension of PSR J0205+6449/3C 58 during off-peak phases was explored as well in the 2PC analysis, but has not been favored over a point-like morphology. Here, using pointlike, we have fitted an extended disk to the off-peak gamma-ray emission of PSR J0205+6449/3C 58, yielding a ΔTS_ext = 0.1; the disk is not favored either. The localization of the off-peak emission determined with pointlike is R.A. = 3140, decl. = 6483, with a 95% confidence error circle of radius 0025, which is consistent with PSR J0205+6449/3C 58. Considering the flat spectrum and the absence of a spectral cutoff, it is natural to propose that the off-peak gamma-ray emission of PSR J0205+6449/3C 58 is dominated by the PWN 3C 58, though a weak magnetospheric component cannot be completely ruled out at low energies. 3C 58 was also detected in the 1FHL and the 3FHL catalogs, and the reported spectra are consistent with our results. The detected morphology of 3C 58 being point-like in 0.1–300 GeV is not unexpected. The arcmin-sized extension in radio and X-rays (Bietenholz et al. 2001; Bocchino et al. 2001) is smaller than the Fermi-LAT PSF (e.g., 01 at 10 GeV).

We have considered a power law with an exponential cutoff for modeling the on-peak emission of PSR J0205+6449. The best-fit parameters are shown in Table 1. The right panel of Figure 2 shows the on-peak TS map of the PSR J0205+6449 region. Its SED is shown in Figure 3. At high energies, the SED is consistent with that derived for the off-peak, which indicates that the PWN 3C 58 dominates the flux. Alternative spectral shapes, like a power law with a subexponential cutoff (dN/dE = N₀(E/E₀)^−Γexp(−E/E₀)^b cm⁻² s⁻¹ GeV⁻¹, leaving the exponential index b free) yield a ΔTS = 5 and are thus not preferred.

We have also modeled the phase-averaged emission of PSR J0205+6449 with a power law having an exponential cutoff. The best-fit parameters are shown in Table 1 and the SED is shown in Figure 3 (right panel). Adopting the best-fit spectral model derived, we calculated the probabilities for photons to come from PSR J0205+6449 within a radius of 3°, using the tool gtsrcprob, and produced a weighted pulsed light curve based on these photons (Kerr 2011b). The bottom panel of Figure 1 shows the folded, pulsed light curve above 100 MeV. The remaining panels of the same figure show the light curve in narrower energy bands. The light curve shows two distinct peaks, which is consistent with the profile reported by Abdo et al. (2009) and the 2PC. To locate the two peaks, we fitted the light curve with two asymmetric Lorentzian functions plus a constant (Figure 1). The fitted constant accounts for the background, which, as we have just shown, is dominated by the PWN 3C 58. The first (P1) and second (P2) peaks are at 0.234 ± 0.003 and 0.719 ± 0.001, respectively. The separation between the two peaks is 0.485 ± 0.003, which is consistent with Abdo et al. (2009). The phase reference used in this paper is different from Abdo et al. (2009). By shifting ∼0.149 spin phase, P2 would be aligned with the profile in Abdo et al. (2009). In that case, the gamma-ray profile we observed is in good alignment with the X-ray profile but is offset from the radio pulse by ∼0.085.

The strength of P1 and P2 is calculated as the sum of the weighted counts during the corresponding on-peak phase (ϕ₁ = 0.184–0.291 and ϕ₂ = 0.574–0.786, respectively) minus the background. The relative strength of P1 and P2 decreases significantly from low to high energies (Figure 4, left panel), as first reported by Abdo et al. (2009). A similar trend was observed in Vela, Crab, Geminga, B1951+32, and J0007+7303 pulsars (Kanbach 1999; Thompson 2001; Aleksić et al. 2014b; Li et al. 2016), which shows a spectral energy dependence of the gamma-ray light curve. We carried out a spectral analysis for the two peaks. In the corresponding on-peak phase for the two peaks (ϕ₁ = 0.184–0.291 and ϕ₂ = 0.574–0.786), the pulsar spectrum is modeled as a power law with an exponential cutoff. The right panel of Figure 4 shows the spectral parameters of PSR J0205+6449 in the two on-peak phases. The lower cutoff energy of P1 when compared with that obtained for P2 (Figure 4, right panel) explains the energy evolution of the P1/P2 ratio.

Zoom In Zoom Out Reset image size

In order to search for Crab-like flares (Abdo et al. 2011; Buehler et al. 2012) from 3C 58, we produced light curves with a 30-day time bin in the 0.1–1 GeV and 1–10 GeV bands during the off-peak phase and in 10–300 GeV during all spin phases. The spectral index is fixed at the best-fit value along the off-peak phases, as listed in Table 1. All the light curve data points are below the detection threshold of TS = 25; therefore, we find no evidence of flaring on this timescale.

Using 8.5 years of Fermi-LAT data and a contemporaneous ephemeris, we carried out a detailed analysis of PSR J0205+6449 both during its off-peak and on-peak phase intervals.

During the off-peak phases, PSR J0205+6449/3C 58 is significantly detected, having a TS value of 202. Its spectrum can be modeled by a simple power law. No extension is detected. The flat spectrum and the nondetection of a spectral cutoff argue for a PWN origin of the off-peak gamma-ray emission of PSR J0205+6449/3C 58. The top panel of Figure 5 shows a theoretical model of the nebula, based on a time-dependent integration of the dynamical evolution of both the nebula and the SNR, the radiation of particles, and the particle population. For details on the model, see Martin et al. (2016) and the appendix in Torres (2017). The spin-down power, the particle injection, the energy losses, and the magnetic field all depend on time, and their dependence is accounted for in the model. The particle content of the nebula is obtained from the balance of energy losses, injection, and escape. We include losses by synchrotron, inverse-Compton (Klein Nishina inverse Compton with the cosmic-microwave background as well as with IR/optical photon fields), self-synchrotron Compton, and bremsstrahlung, devoid of any radiative approximations, and compute the radiation produced by each process. The model also considers the dynamical influence of the reverse shock traveling backward toward the pulsar, and compressing the nebula. However, given that the pulsar is young, this effect is not found to be relevant and results would be very similar if neglected altogether: the nebula is freely expanding. We find a good agreement with data considering a nebula at 3.2 kpc (the same distance we use below for the computation of the pulsar magnetospheric power) and an age of 2500 years. The fitting model features a broken power law for injected electrons at the termination shock, with a low (high) energy index of 1.1 and 2.94, and an energy break at Lorentz factor 9 × 10⁴. These values of the parameters (as well as the magnetic energy fraction) are in agreement with our earlier analysis (Torres et al. 2013). The magnetic fraction (the fraction of spindown energy that goes into the magnetic field) is 0.2; thus 3C 58 is a particle-dominated nebula; nevertheless, the magnetic field has one of the highest energy reservoirs when compared with all other nebulae of similar age, perhaps with the exception of CTA 1 (see Torres et al. 2014). Using the current data, we have seen that there is a degeneracy regarding which inverse Compton contribution dominates at high energies. We have explored about 1000 models varying the energy densities and temperatures for the NIR and FIR photon backgrounds and the best-fitting one has similar contributions of both. Thus, we can find models where one or the other dominates without changing the overall fit significantly (e.g., within a factor of 1.3 of the miminum χ²). Further data would be needed to distinguish among these possibilities; in particular, 3C 58 will be a bright source for the Cerenkov Telescope Array. Observations with this facility will help determine the peak and the fall-off of the gamma-ray emission, distinguishing between NIR- or FIR-dominated scenarios.

Zoom In Zoom Out Reset image size

For the on-peak interval, PSR J0205+6449 can be modeled by a power law with an exponential cutoff. We explored the existence of a subexponential cutoff, but no improvement was found. PSR J0205+6449 shows a two-peak pulse profile. The ratio of P1 and P2 decreases significantly with energy (Figure 4, left panel). This is consistent with the cutoff energy of P1 being lower than that obtained for P2 (Figure 4, left panel).

The most common interpretation of the magnetospheric radiation for this and all other gamma-ray pulsars is that it originates in synchro-curvature radiation in a high-altitude gap (alternatively, see, e.g., Cerutti et al. 2015). To consider this statement further, we have applied the model discussed in detail in Vigano & Torres (2015) and Vigano et al. (2015). In this model, the detectable radiation coming from the magnetosphere can be estimated by integrating the single-particle synchro-curvature spectrum along the traveled distance in the gap, convolved with an effective particle distribution, dN/dx, where x is the distance along the magnetic field line (see Vigano et al. 2015 for more details on the computation of the synchro-curvature power). In this model, two parameters, the electric field ${E}_{| | }$ accelerating particles, and the degree of uniformity of the particle's distribution along the trajectory, x₀/R_lc, define the spectrum completely. The normalization N₀ moves the spectrum up and down in luminosity, without modifying its shape. We have applied this model to the phase-averaged SED (i.e., Fermi-LAT data only and a fixed magnetic gradient b = 2.5, as in Vigano et al. 2015) and find a good match for the values shown in Figure 5 (bottom panel). The fitting results are consistent with those in Vigano et al. (2015). Our current analysis gives $\mathrm{log}{N}_{0}={30.90}_{-0.30}^{+0.21}$ compared to the value of $\mathrm{log}{N}_{0}={29.93}_{-0.33}^{+0.67}$ in Vigano et al. (2015). Whereas applying the model to the on-peak phase spectrum would only lead to comparable results (given that the spectral index of both the phase averaged and on-peak spectra are the same within errors, and the cutoff energy is only slightly different), the phase-averaged fit allows for a comparison with other gamma-ray pulsars for which only averaged spectral results are currently available. PSR J0205+6449's fit parameters (that are well in agreement with those found using earlier data—see Table 2 of Vigano et al. 2015) confirm the reported trend relating ${E}_{| | }$ to x₀/R_lc (see the first panels in Figure 2 of Vigano et al. 2015), emphasizing that synchro-curvature dominated radiation is likely behind the pulsations. We note that PSR J0205+6449 is part of the soft gamma-ray pulsar catalog and it has also been detected in X-rays (Kuiper & Hermsen 2015). For a complete model of the pulsed spectrum considering also the X-ray data, and further discussion, see Torres (2018).

The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515.

We acknowledge the assistance from Dr. M. Kerr and Dr. D. Smith with the gamma-ray ephemeris for PSR J0205+6449 and Dr. P. Saz Parkinson for discussions. We acknowledge the support from The National Key Research and Development Program of China (2016YFA0400800), the grants AYA2015-71042-P, SGR 2014-1073 and the National Natural Science Foundation of China via NSFC-11473027, NSFC-11503078, NSFC-11673013, NSFC-11733009, XTP project XDA 04060604, and the Strategic Priority Research Program "The Emergence of Cosmological Structures" of the Chinese Academy of Sciences, grant No. XDB09000000, as well as the CERCA Programme of the Generalitat de Catalunya. Work at NRL is supported by NASA.