PSR J2021+4026 IN THE GAMMA CYGNI REGION: THE FIRST VARIABLE γ-RAY PULSAR SEEN BY THE Fermi LAT

Pulsars are the largest high-energy γ-ray source class in the Milky Way, with 117 characterized in the second Fermi Large Area Telescope (LAT) pulsar catalog (2PC; Abdo et al. 2013). Early γ-ray observations suggested glitch-associated pulse and flux changes in the Crab (Greisen et al. 1975) and Vela (Grenier et al. 1988) pulsars, unconfirmed by additional observations (Nolan et al. 2003). Steady γ-ray fluxes and pulse profiles on timescales longer than those needed for pulsar detections has been an axiom (e.g., Nolan et al. 2012, hereafter 2FGL).

PSR J2021+4026 (hereafter J2021+4026) was discovered in a blind frequency search using LAT data (Abdo et al. 2009). Its spin frequency f ∼ 3.8 Hz and frequency derivative $\dot{f}\sim -8 \times 10^{-13}$ Hz s⁻¹ point to a young, energetic pulsar (characteristic age τ_c = 77 kyr, spindown power $\dot{E}_\mathrm{SD}\sim 10^{35}$ erg s⁻¹). Radio and optical searches did not yield any plausible counterparts, while deep observations with Chandra and XMM-Newton led to an association with the X-ray source S20 (Weisskopf et al. 2011; Trepl et al. 2010), from which X-ray pulsations were recently detected (Lin et al. 2013). J2021+4026 is seen within the radio shell of the supernova remnant (SNR) G78.2+2.1 (e.g., Ladouceur & Pineault 2008), also an extended γ-ray source with ∼06 radius (Lande et al. 2012). A tentative association of these two sources implies a pulsar distance of ∼1.5 kpc.

Using AGILE observations from 2007 November to 2009 August, Chen et al. (2011) reported variability in 1AGL J2022+4032, positionally coincident with J2021+4026, but concluded that it was more likely due to another source along the line of sight. This Letter reports the Fermi LAT observation of a discrete change in the γ-ray flux and frequency derivative of J2021+4026 and further results from its long-term monitoring.

We analyzed ∼52 months of LAT data from 2008 August 4 to 2012 December 11. We selected P7REP_SOURCE "photon" class events (Bregeon et al. 2013; Ackermann et al. 2013b) in a 15°radius Region Of Interest (ROI) with energies from 100 MeV to 300 GeV and zenith angles ⩽100°. We excluded time intervals when the LAT rocking angle was >52° or the zone defined by the zenith cut intersected the ROI.

We characterized the spectrum of J2021+4026 through a binned likelihood fit over a 14° × 14° region centered on the pulsar, with 01 angular grid. We used logarithmically spaced energy bins from 100 MeV to 300 GeV (16 below 10 GeV, 8 above). The combined likelihood technique (e.g., Ackermann et al. 2012b) treats photons converting in the front and back tracker sections separately to exploit the former's higher angular resolution (Atwood et al. 2009). We used the P7REP_SOURCE_V15 LAT instrument response functions (IRFs) and associated diffuse emission models⁵⁴ and LAT Science Tools v09r32p01. The isotropic background and the residual Earth limb emission were kept fixed. The Galactic diffuse emission was multiplied by a power law with free normalization and spectral index. Our model includes the known extended sources in the region, i.e., the Cygnus-X cocoon (Ackermann et al. 2011) and SNR G78.2+2.1 (Lande et al. 2012), not included in previous analyses (2FGL, 2PC). We also included 2PC pulsars, and additional 2FGL sources with high significance (average Test Statistic⁵⁵ TS > 100), or within 4° of J2021+4026. Known flaring sources in the region (Ackermann et al. 2013a) are included in the model. Source spectra were modeled using the functional forms described in the catalogs with all the spectral parameters free. We then searched a TS map of the region for excesses, used as seeds to determine positions and spectra of new sources (see 2FGL). We found five, modeled with log-parabola spectra, at the epoch J2000 positions $(\mathrm{R.A}.,\mathrm{{\rm d}ec{\rm l.}}) = (310{.\!\!^\circ}52, 42{.\!\!^\circ}06)$, (30666, 4005), (30965, 4222), (30867, 4305), (31244, 4438).

J2021+4026's position was set to Chandra S20's. Its γ-ray spectrum was modeled using a power law with an exponential cutoff (PLEC; Equation 1), with b free for the phase-averaged spectral analysis and b = 1 (PLEC1) for the phase-resolved analysis in Section 5:

$$\begin{equation} \frac{{d}N}{{d}E}= N_0 \left(\frac{E}{E_0}\right)^{-\Gamma } \exp \left[-\left(\frac{E}{E_\mathrm{c}}\right)^b\right]. \end{equation} \tag{ 1 }$$

PLEC1 represents expectations for high-altitude magnetospheric emission. Phase-averaged spectra are usually better fit with b < 1 due to the superposition of several PLEC1 components with different photon indices and cutoff energies (Abdo et al. 2010).

Following Chen et al. (2011), we searched for flux variability around J2021+4026 at energies E > 100 MeV and E > 1 GeV, applying the method in 2FGL. We first fit the data over the entire time range. Then, we divided the range into 7 and 30 day time bins and refit, allowing free normalizations for all sources and fixing the other spectral parameters to their long-term average. The fit was then repeated in each time bin by also fixing the source of interest's normalization to its long-term average. The Galactic diffuse normalization was fit in each time bin, verifying a posteriori its compatibility with a constant. Following 2FGL Equation (4), the fit maximum likelihood values established the probability P that the observed fluctuations are stochastic only. A 2% flux systematic error accounts for exposure uncertainties between the different epochs.

We applied this procedure to J2021+4026, to SNR G78.2+2.1 (seen in the same direction, but separable from the pulsar at high energies due to its extension⁵⁶) and to J2021+3651, located 35 away, with spectrum and flux similar to J2021+4026's. Both SNR G78.2+2.1 and J2021+3651 show constant fluxes (P > 0.4 in all cases), while J2021+4026 shows significant variability: P_30 days = 7 × 10⁻⁸ (1 × 10⁻¹⁰), P_7 days = 6 × 10⁻³ (2 × 10⁻⁴) at energies >100 MeV (>1 GeV).

Figure 1 shows the >1 GeV energy flux for the three sources in 30 day bins. J2021+4026 shows an abrupt ∼20% flux decrease near MJD 55850 (Table 1), confirmed at energies >100 MeV and for 7 day bins. We exclude that this drop is due to systematic effects since there is no analogous drop for the two other sources, observed simultaneously. No significant changes in the Fermi observing strategy occurred near MJD 55850. We verified that fixing J2021+4026's normalization to its average yields negative residuals consistent with a point-like source at the pulsar position in all energy bands after MJD 55850.

Zoom In Zoom Out Reset image size

Table 1. J2021+4026's Properties Before and After the Jump

Time Range	54682–55850	55850–56273
(MJD)
Number of days	1167	423
Fγa >0.1 GeV	8.33 ± 0.08	6.86 ± 0.13
Fγa >1 GeV	3.57 ± 0.05	2.74 ± 0.06
$\dot{f}$b	−7.6978 ± 0.0007	−8.166 ± 0.002
δP1c	0.19 ± 0.02	0.13 ± 0.02
Δ12d	0.505 ± 0.005	0.565 ± 0.006
δP2c	0.176 ± 0.007	0.174 ± 0.006
Δ1BRd	0.229 ± 0.008	...
δBRc	0.11 ± 0.02	...
P1/P2e	0.54 ± 0.06	0.24 ± 0.03
BR/P2e	0.16 ± 0.03	...
Constant/P2e	1.83 ± 0.14	1.09 ± 0.06

Notes. Statistical uncertainties only. For details on parameters, see Sections 4 and 5. ^a10⁻¹⁰ erg cm⁻² s⁻¹. ^bAt the reference epoch for the two timing solutions, 10⁻¹³ Hz s⁻¹. ^cPeak FWHM (E > 0.1 GeV). ^dPhase lag between peaks (E > 0.1 GeV). ^eRatios of the peak amplitudes or constant-level-to-P2 amplitude (E > 0.1 GeV).

Download table as:  ASCIITypeset image

Figure 1 also suggests a steady flux increase for J2021+4026 before the drop near MJD 55850. A χ² fit of a linear function of time versus flux in 30 day bins⁵⁷ before MJD 55850 gives a 4% ± 2% yr⁻¹ flux increase, preferred over the constant flux hypothesis at the ∼3.3σ level for both >100 MeV and >1 GeV. The χ² test applied to J2021+3651 favors a constant flux. We further assessed this trend independently of any functional dependency using the Kendall rank correlation test. We obtain a Kendall coefficient τ = 0.78 (0.71) for J2021+4026 for >100 MeV (>1 GeV): the probability⁵⁸ of this coming from stochastic fluctuations of a steady flux is 0.0005 (0.01), indicating a monotonic increase with time at the ∼3.5σ (∼3.2σ) level. For J2021+3651, τ = 0.11 (0.09): the same probability is 0.38 (0.41).

To investigate the origin of the flux drop, we monitored the evolution of the pulsar timing parameters. We divided the entire time range into 60 day bins, where pulsations are clearly detectable, yet we can neglect the timing noise and approximate the frequency evolution as a linear function,

$$\begin{equation} f(t)=f_0 + f_1\times (t-t_0)\,. \end{equation} \tag{ 2 }$$

For each bin, we used the $Z^{2}_{n}$ (n = 4) test (Buccheri et al. 1983) to search the f₀–f₁ space for periodicity (f₀ and f₁ represent the frequency f and frequency derivative $\dot{f}$, respectively, in each bin). Figure 1 shows that near MJD 55850, $\dot{f}$ suddenly decreases by ∼5 × 10⁻¹⁴ Hz s⁻¹, i.e., ∼4% of the initial value. This is reflected as a change of the frequency evolution slope, while f does not change appreciably. The $\dot{f}$ change is simultaneous with the flux decrease, strongly suggesting that the flux change is from the pulsar itself rather than another source along the line of sight. This is strengthened by the results for higher energy and narrower time bins (Figure 2), suggesting that the flux variation occurred within a week or less. We also explored three-day and one-day binning, but count rates are too low to measure when and how quickly the flux change occurred. The data hint that it happened within a few days after MJD 55850.

Zoom In Zoom Out Reset image size

The frequency derivative discontinuity resembles a glitch in $\dot{f}$ (e.g., Cordes & Downs 1985). However, the phenomenology differs from radio and γ-ray glitches (e.g., Espinoza et al. 2011; Pletsch et al. 2012): glitches are not usually associated with a flux change and are followed by a recovery, not detected for J2021+4026 prior to MJD 56200.

Doppler shift due to pulsar motion in a binary system cannot explain the change in $\dot{f}$. If we assume the pulsar moved in the same direction for the ∼3 yr before the jump, i.e., half of a circular orbit with radius 6 (1) a.u., that would yield a fractional frequency change due to the Doppler shift of 10⁻⁵ (10⁻⁴), compared to the observed ∼4% variation. Reproducing the observed change for a six-year orbit requires a highly eccentric orbit with an unrealistically small minor axis of 0.01 a.u. Therefore, the $\dot{f}$ change is likely related to some phenomenon in the pulsar magnetosphere.

The jump causes phase coherence loss. We therefore built two timing solutions using LAT γ rays (Ray et al. 2011). We used 32 day intervals to determine pulse times-of-arrival (TOAs).⁵⁹ We obtained 36 (13) before (after) the jump. We used the TEMPO2 package (Hobbs et al. 2006) to fit these TOAs using a model with absolute phase, frequency, and its first three derivatives at the reference epoch. The rms of the timing residuals of the post-jump timing solution is 2.1 ms. The pre-jump solution needed whitening with sinusoidal waves to achieve a 3.0 ms residual rms. We verified that this is due only to the different lengths of the time ranges. The timing solutions⁶⁰ confirm the sudden $\dot{f}$ change near MJD 55850 (Figures 1 and 2). Owing to similarities with Geminga, we shifted photon phases to center the second highest peak at 0.1, resulting in a half-period shift relative to 2PC.

We studied the pulse profile and the spectrum before and after MJD 55850, repeating the likelihood analysis of Section 2 for the two time intervals independently. Then we selected γ rays within 2°of the pulsar and used the best-fit spectrum to assign each photon a weight, the probability of being associated with the pulsar. We assigned each photon a phase using the timing solutions described in Section 4, and thus built weighted pulse profiles for different energy bands (Figure 3). The profiles show two main peaks, P1 and P2, interconnected by a bridge (BR), where a third peak appears before the jump, especially at E > 1 GeV. An off-peak region, OP, follows P2.

Zoom In Zoom Out Reset image size

We fit the pulse profile peaks with three Gaussians, adding a constant to account for steady emission (Figure 3). The third peak in BR is included only when detected with >3σ significance. Table 1 summarizes the fit results and other pulsar characteristics. Across the jump the constant component decreases compared to P2's amplitude, and the P1–P2 lag (Δ₁₂) increases. The peak height ratio P1/P2 shows hints of a decrease, while the third peak significance decreases⁶¹ from 8σ to 2.8σ (5.5σ to 2.2σ) for >100 MeV (>1 GeV).

To determine the pulsar spectral energy distribution, we subdivided the data before and after the jump into four phase intervals: 0–0.25 (P1), 0.25–0.4 (BR), 0.4–0.8 (P2), and 0.8–1 (OP). For each, we determined the pulsar spectrum over 10 logarithmically distributed energy bins from 100 MeV to 10 GeV, approximating the spectrum within each bin with a power law with spectral index 2. We also determined the spectral energy distribution over the entire energy band using a PLEC1 model.

As shown in Figure 4, across the jump the flux varies at all phases but P2, strengthening the association of the flux drop with the pulsar as opposed to another source. The OP spectrum is always well-described by a power law with an exponential cutoff at ∼2 GeV, indicating a magnetospheric origin over all phases. There is an indication of a decrease in E_c for P1 (∼2σ).

Zoom In Zoom Out Reset image size

We detected a "jump," a sudden decrease of J2021+4026's flux above 100 MeV of ∼20% associated with a ∼4% increase in spindown rate on a timescale shorter than one week. The jump is also accompanied by changes in the pulse profile. Furthermore, we found evidence for a small, steady flux increase preceding the jump. The temporal correlation between spindown and flux changes strongly indicates that these phenomena are related to the pulsar. While mode changes and other intermittent behavior are well known for some radio pulsars (e.g., Lyne et al. 2010), this is the first time such behavior has been seen at γ-ray energies.

J2021+4026 belongs to a small set of unusual LAT pulsars—PSR J0633+1746 (Geminga), J1836+5925 and J2021+4026—the sources in 2PC with the brightest magnetospheric emission at all spin phases. They are all radio-quiet, with phase lags between the main peaks Δ > 0.5, higher than typical (Abdo et al. 2013). Finally, although only Geminga has a parallax distance (we rely on the SNR association of J2021+4026 and X-ray spectral arguments for J1836+5925), if we adopt the common assumption that the γ-ray pulse is effectively uniform on the sky, beaming factor $f_\Omega =1$, then all three have large efficiencies $\eta = 4\pi f_\Omega F_\gamma d^2/ \dot{E}_\mathrm{SD} \ge 1$ (d is the distance and F_γ is the energy flux; see Table 1). J2021+4026 is the most extreme of the three, with η = 2.3.

All of these attributes point to peculiarity in the γ-ray beaming. They are most easily understood in the context of the classical outer gap (OG) model. Romani & Watters (2010) show that such large peak lag implies small magnetic inclinations α < 30° and near-equatorial viewing angles 80° < ζ < 100°. For this geometry the pulsars should be radio-quiet, the OG geometry predicts $f_\Omega \approx 0.1\hbox{--}0.2$ (η < 1) and the Earth line-of-sight skims nearly tangentially to the peak caustics, producing complex peak structure and strong off-peak emission (Romani & Watters 2010). Also, two-pole caustic (TPC) models (Dyks & Rudak 2003) can produce strong off-peak components for a wider range of geometries (most with α < 30°). These models tend to have single broad pulses at small ζ, but large ζ models can be double pulsed. Thus, the preferred geometry is similar to that of the OG case, and should also be radio-quiet. These models have $f_\Omega \approx 0.5\hbox{--}0.75$, making it harder to accommodate the observed γ-ray flux. If classical TPC solutions are extended to higher altitude, then one may recover the broad equatorial pulses and small $f_\Omega$ (M. Pierbattista et al. 2013, in preparation). The nearly aligned rotator viewed at high inclination scenario is independently confirmed for Geminga due to X-ray observations of its rotating hot spot (Caraveo et al. 2004).

When emission from near the light cylinder dominates the pulse, the concentration of the γ-ray beam to a narrow equatorial strip gives high apparent η and allows small changes in magnetic field morphology or even in α to move a substantial fraction of the γ-ray beam over the line of sight, giving large fractional changes to the pulse profile and $f_\Omega$. For young pulsars, we expect the γ-ray luminosity to scale with $\sqrt{\vphantom{A^A}\smash{{{\dot{E}_\mathrm{SD}}}}}$ (e.g., Harding 1981; Abdo et al. 2013). The decrease in flux rate associated to an increase in spindown rate after the jump strengthens the case that beaming must play a key role.

Therefore, we can speculate that the jump of J2021+4026 represents a shift in the magnetic field structure, leading to either an effective α change or an effective current change. These may be precipitated by a reconfiguration of field line footpoints at the surface, i.e., in the crustal layers, that modifies the overall magnetic dipole torque on the star. There is no reason to expect that the resulting spindown increase should enhance the solid-angle integrated luminosity of the pulsar γ-ray emission, since the principal effect is that of a modified beaming. If the slow variation in the pulsar flux before the jump is substantiated by additional study, this might plausibly be associated with more gradual changes in geometry, for example, from force-free precession (e.g., Jones 2012).

The very sensitivity of the beaming to currents and geometry for the equatorial, small α, OG, or TPC models complicates the interpretation of the observations in terms of magnetosphere configurations. Alternative tests of these scenarios may rely on non-γ-ray constraints on spin geometry, e.g., from X-ray imaging of the synchrotron termination shock (Ng & Romani 2004), or, if radio-loud examples can be detected, from polarization studies.

Radio and X-ray observations have shown that mode changes, and variability in general, are key to understanding pulsars (Lyne et al. 2010; Hermsen et al. 2013) and therefore to investigating their fundamental physics (Alpar et al. 1984; Cordes et al. 2004). The "jump" in J2021+4026 breaks the axiom of pulsars as steady γ-ray emitters, opening new avenues for pulsar variability studies at γ-ray energies, where the bulk of their spindown energy is emitted.

The Fermi LAT Collaboration acknowledges support from a number of agencies and institutes for both the development and the operation of the LAT as well as scientific data analysis. These include NASA and DOE in the United States, CEA/Irfu and IN2P3/CNRS in France, ASI and INFN in Italy, MEXT, KEK, and JAXA in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council, and the National Space Board in Sweden. Additional support from INAF in Italy and CNES in France for science analysis during the operations phase is also gratefully acknowledged.