A DETAILED X-RAY INVESTIGATION OF PSR J2021+4026 AND THE γ-CYGNI SUPERNOVA REMNANT

The Cygnus region, which is located ∼80° from the Galactic center, is one of the most complex γ-ray structures on the Galactic plane. Many interesting sources are found to reside in this region, including an OB association, a microquasar, supernova remnants, and rotation-powered pulsars. Since the launch of the Fermi Gamma-Ray Space Telescope, many new γ-ray pulsars have been uncovered in the Cygnus region either in blind searches or by folding the γ-ray data with the timing ephemeris determined through radio observations (Abdo et al. 2013). Among them, the most intriguing one is PSR J2021+4026, which belongs to the growing class of radio-quiet γ-ray pulsars.

PSR J2021+4026 is a bright γ-ray source that was detected at a significance of >10σ with only the first three months Large Area Telescope (LAT) data (Abdo et al. 2009a; 0FGL). Its timing ephemeris was subsequently reported (Abdo et al. 2009b, 2010). It has a rotation period and first period derivative of P = 265 ms and $\dot{P}=5.48\times 10^{-14}$ s s⁻¹, respectively, which imply a spin-down age of ∼77 kyr, surface dipolar magnetic field strength of ∼4 × 10¹² G, and a spin-down power of ∼10³⁵ erg s⁻¹.

Recently, the long-term γ-ray monitoring of PSR J2021+4026 has discovered its γ-ray flux at energies >100 MeV has suddenly decreased by ∼18% near MJD 55850 (Allafort et al. 2013). This flux jump was accompanied by the change in the γ-ray pulse profile and the spin-down rate. These attributes make PSR J2021+4026 the first variable γ-ray pulsar ever observed.

Despite the efforts devoted to searching for its radio counterpart, no radio pulsar associated with PSR J2021+4026 has been detected (Becker et al. 2004; Trepl et al. 2010; Ray et al. 2011). Without radio detection, observations in the X-ray regime are particularly important for constraining the nature of the emission of radio-quiet γ-ray pulsars. Soon after the detection of PSR J2021+4026 was reported by Abdo et al. (2009a), we searched for its X-ray counterpart by using all the available archival data (Trepl et al. 2010). Within the 95% confidence circle of PSR J2021+4026 at that time (0FGL J2021.5+4026; see Abdo et al. 2009a), we reported the identification of the X-ray source 2XMM J202131.0+402645/CXOU 202130.55+402646.9 from XMM-Newton (ObsID: 150960801) and Chandra (ObsID: 5533) data as the possible counterpart of PSR J2021+4026. This source was found to be the only non-variable X-ray object without any optical/IR counterpart within the γ-ray error circle. The association between the pulsar and this source is reinforced by the fact that the X-ray position is consistent with the optimal γ-ray timing solution (Trepl et al. 2010). These results have been confirmed by a follow-up 56 ks Chandra observation during a dedicated investigation (ObsID 11235; Weisskopf et al. 2011). This observation also enables us to examine the X-ray spectrum of this promising counterpart of PSR J2021+4026 and we found that it possibly contains both thermal and non-thermal contributions.

Although all the aforementioned investigations strongly suggest the association between 2XMM J202131.0+402645 and PSR J2021+4026, the physical connection between this X-ray source and the γ-ray pulsar could not be confirmed unambiguously until the X-ray pulsation was recently discovered by us with a deep XMM-Newton observation (Lin et al. 2013). The consistency between the detected X-ray periodicity and the γ-ray pulsation in the same epoch has eventually pinpointed the long-sought connection between 2XMM J202131.0+402645 and PSR J2021+4026. In this paper, we reported a further X-ray analysis of PSR J2021+4026 by disentangling its pulsed and unpulsed components. Also, by re-examining archival on-axis Chandra observations, we have searched for the possible pulsar wind nebula (PWN) around PSR J2021+4026.

Apart from PSR J2021+4026, our XMM-Newton observation also covers the central region and part of the southeastern rim of the supernova remnant (SNR) γ-Cygni (G78.2+2.1), which has been suggested to be associated with PSR J2021+4026 (Trepl et al. 2010). G78.2+2.1 is rather extended. Its radio and X-ray shells have a size of ∼1° (Leahy et al. 2013). The X-ray emission from G78.2+2.1 has been investigated several times with ROSAT, ASCA, and Chandra (Lozinskaya et al. 2002; Uchiyama et al. 2002; Aliu et al. 2013; Leahy et al. 2013). The most recent X-ray imaging spectroscopic analysis of G78.2+2.1 was reported by Leahy et al. (2013). By using the archival Chandra data, the authors have examined the diffuse X-ray emission from the northern part and the central region of G78.2+2.1. While the column absorption for these two spatial regions is found to be comparable, the plasma temperature of the central region, ∼10⁷ K, is suggested to be higher. However, the uncertainties of the spectral parameters are too large to draw a firm conclusion. This is can be ascribed to the low surface brightness of the central region (see Figure 4 in Leahy et al. 2013) and the relatively inferior collecting power of Chandra.

The southeastern rim contributes ∼60% of the radio flux of this SNR, which is known as DR4 (Downes & Rinehart 1966). While this is the brightest part in the radio, it is very dim in the X-ray regime (see Figure 1 in Uchiyama et al. 2002). This part has been excluded in the updated investigation by Leahy et al. (2013), as there is no existing Chandra data that cover this region. With our deep XMM-Newton observation, we are able to provide tighter constraints on the nature of the emission of the central and southeastern parts of G78.2+2.1 and its possible association with PSR J2021+4026.

2.1. Observations and Data Reduction

2.2. Spatial Analysis

The X-ray color images as obtained by MOS1/2 and PN are shown in Figures 1 and 2 respectively (red: 0.5–1 keV; green: 1–2 keV; blue: 2–10 keV). In order to correct for the non-uniformity across the detector and the mirror vignetting, each image has been normalized by the exposure map generated by using the XMMSAS task eexpmap for the corresponding detector. The exposure-corrected images have been adaptively smoothed so as to attain a minimum signal-to-noise ratio of 3. This deep observation allows us to search for new X-ray point sources in this field. In Figure 1, point sources with various X-ray hardness can be seen. For determining their positions and count rates, we performed a source detection by using maximum likelihood fitting on MOS1, MOS2, and PN data individually with the aid of the XMMSAS task edetect_chain. We set the detection threshold to be 4σ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By visual inspection, we removed the weak sources that are potentially false detections from the source lists resulting from individual cameras. These include several sources close to the edge of the FoV (one from MOS1, one from MOS2, and two from PN) as well as a few slightly extended sources coinciding with the diffuse emission that are likely to be clumps of the supernova remnant (three from MOS1 and two from MOS2). The screened lists are subsequently correlated and merged by using the XMMSAS task srcmatch. In the case that the position of a source obtained from two detections is consistent within its 3σ uncertainties, the two detections are merged as a single entry. Their source properties are summarized in Table 1 ⁶. Including PSR J2021+4026 (i.e., source 18), 42 distinct point sources are detected within the 30' FoV around the aim-point in this observation. Among all these X-ray point-like objects in this field, source 8 (XMM J202154.7+402855) is the brightest. A further analysis found that its X-ray flux is significantly variable, which is worth a deeper investigation. The spectral and temporal analyses of XMM J202154.7+402855 are reported in the Appendix.

Table 1. X-Ray Sources Detected in the Southeastern Field of G78.2+2.1 as Labeled in Figure 1

Sourcea	R.A. (J2000)	Decl. (J2000)	r1σb	Net Count Rate			Remarkc
No.	(h:m:s)	(d:m:s)	(arcsec)	MOS1	MOS2	PN
				(counts ks−1)	(counts ks−1)	(counts ks−1)
1	20:22:21.77	40:31:17.93	0.43	7.11 ± 0.57	7.27 ± 0.60	⋅⋅⋅	⋅⋅⋅
2	20:22:15.89	40:28:30.17	0.44	3.01 ± 0.46	3.45 ± 0.39	⋅⋅⋅	⋅⋅⋅
3	20:22:09.89	40:29:58.75	0.67	1.06 ± 0.28	1.64 ± 0.28	⋅⋅⋅	⋅⋅⋅
4	20:22:02.69	40:26:08.79	0.47	1.17 ± 0.21	1.25 ± 0.30	⋅⋅⋅	⋅⋅⋅
5	20:21:58.23	40:30:51.61	0.48	2.93 ± 0.31	2.51 ± 0.29	⋅⋅⋅	⋅⋅⋅
6	20:21:57.57	40:26:48.07	0.52	1.25 ± 0.20	1.48 ± 0.23	⋅⋅⋅	W44
7	20:21:54.78	40:24:34.84	0.69	1.34 ± 0.22	1.38 ± 0.24	⋅⋅⋅	⋅⋅⋅
8	20:21:54.65	40:28:55.28	0.21	11.69 ± 0.49	9.83 ± 0.49	⋅⋅⋅	⋅⋅⋅
9	20:21:52.95	40:24:19.73	0.78	0.69 ± 0.19	0.77 ± 0.18	⋅⋅⋅	⋅⋅⋅
10	20:21:50.50	40:18:32.05	0.55	3.62 ± 0.38	5.37 ± 0.50	⋅⋅⋅	⋅⋅⋅
11	20:21:48.07	40:23:40.65	0.77	0.94 ± 0.21	1.10 ± 0.22	⋅⋅⋅	⋅⋅⋅
12	20:21:42.89	40:23:54.06	0.72	0.75 ± 0.17	0.66 ± 0.18	⋅⋅⋅	W38
13	20:21:38.20	40:24:43.18	0.60	1.40 ± 0.19	1.45 ± 0.22	⋅⋅⋅	W34
14	20:21:38.05	40:29:36.34	0.30	0.88 ± 0.18	1.30 ± 0.22	3.30 ± 0.42	W33
15	20:21:37.18	40:29:58.72	0.32	2.53 ± 0.26	2.50 ± 0.27	7.65 ± 0.52	W32
16	20:21:34.81	40:28:35.74	1.23	0.00 ± 0.00	0.00 ± 0.00	1.34 ± 0.32	W30
17	20:21:33.30	40:29:09.82	0.47	0.86 ± 0.18	1.10 ± 0.25	3.21 ± 0.39	W27
18	20:21:30.48	40:26:46.30	0.22	4.79 ± 0.32	4.45 ± 0.33	12.90 ± 0.64	W20
19	20:21:29.95	40:29:48.68	0.38	1.70 ± 0.23	1.52 ± 0.23	4.89 ± 0.46	W19
20	20:21:20.07	40:17:27.78	0.46	5.52 ± 0.50	5.68 ± 0.55	⋅⋅⋅	⋅⋅⋅
21	20:21:11.10	40:28:04.43	0.21	9.13 ± 0.44	7.82 ± 0.43	⋅⋅⋅	W3
22	20:21:10.21	40:30:53.78	0.55	1.64 ± 0.25	1.97 ± 0.26	⋅⋅⋅	⋅⋅⋅
23	20:21:06.54	40:18:44.75	0.49	2.55 ± 0.33	3.23 ± 0.42	⋅⋅⋅	⋅⋅⋅
24	20:21:01.74	40:34:42.45	0.38	5.74 ± 0.52	5.06 ± 0.51	⋅⋅⋅	⋅⋅⋅
25	20:20:59.24	40:31:42.71	0.63	1.37 ± 0.26	1.24 ± 0.26	⋅⋅⋅	⋅⋅⋅
26	20:20:57.37	40:28:27.62	0.34	4.49 ± 0.37	4.42 ± 0.40	⋅⋅⋅	⋅⋅⋅
27	20:20:57.01	40:33:26.54	1.10	1.40 ± 0.31	1.11 ± 0.28	⋅⋅⋅	⋅⋅⋅
28	20:20:54.62	40:28:21.02	1.10	0.71 ± 0.19	1.12 ± 0.26	⋅⋅⋅	⋅⋅⋅
29	20:20:53.15	40:30:24.83	0.51	2.02 ± 0.31	2.23 ± 0.35	⋅⋅⋅	⋅⋅⋅
30	20:20:52.33	40:24:29.23	0.29	3.04 ± 0.35	3.54 ± 0.49	⋅⋅⋅	⋅⋅⋅
31	20:20:52.09	40:28:26.67	0.46	4.53 ± 0.42	3.61 ± 0.40	⋅⋅⋅	⋅⋅⋅
32	20:20:51.11	40:30:31.22	0.84	1.30 ± 0.27	1.32 ± 0.31	⋅⋅⋅	⋅⋅⋅
33	20:20:40.65	40:27:03.54	0.40	5.89 ± 0.53	4.51 ± 0.49	⋅⋅⋅	⋅⋅⋅
34	20:20:33.41	40:18:25.60	0.69	5.31 ± 0.82	6.74 ± 0.92	⋅⋅⋅	⋅⋅⋅
35	20:22:00.96	40:33:41.04	0.89	⋅⋅⋅	1.33 0.28	⋅⋅⋅	⋅⋅⋅
36	20:21:55.92	40:38:31.92	0.85	⋅⋅⋅	0.84 0.30	⋅⋅⋅	⋅⋅⋅
37	20:21:51.36	40:32:30.48	1.31	⋅⋅⋅	1.48 0.25	⋅⋅⋅	⋅⋅⋅
38	20:21:37.44	40:36:40.68	2.40	⋅⋅⋅	1.10 0.31	⋅⋅⋅	⋅⋅⋅
39	20:21:25.44	40:36:14.40	1.53	⋅⋅⋅	1.33 0.30	⋅⋅⋅	⋅⋅⋅
40	20:21:25.20	40:38:15.00	0.78	⋅⋅⋅	1.37 0.30	⋅⋅⋅	⋅⋅⋅
41	20:21:20.16	40:34:55.92	1.28	⋅⋅⋅	0.68 0.20	⋅⋅⋅	⋅⋅⋅
42	20:21:07.20	40:38:53.52	0.74	⋅⋅⋅	1.29 0.35	⋅⋅⋅	⋅⋅⋅

Notes. ^aSee Figure 1. ^b1σ positional uncertainty. ^cWx indicates the source detected independently by Weisskopf et al. (2011) with x corresponding to the label in their Figure 1.

Download table as:  ASCIITypeset image

Besides the point-source population, faint diffuse X-ray structures have also been seen in this observation. At the southeastern edge of the MOS1/2 image (Figure 1), extended soft emission is highlighted by the elliptical region. This region partially covers the structure R1 examined by Uchiyama et al. (2002). In the central region, an extended structure around PSR J2021+4026 can be seen in both Figures 1 and 2. A close-up of this feature as seen by the PN camera is shown in Figure 2. The solid line ellipses in Figures 1 and 2 illustrate the source regions for extracting the spectra from the southeastern and the central parts, respectively, of G78.2+2.1 (see Section 2.4).

Apart from reporting the discovery of X-ray pulsation, Lin et al. (2013) also found that the phase-averaged spectrum can be described by a blackbody-plus-power-law (BB+PL) model. One possible origin for the PL component is the PWN. Although this XMM-Newton observation provides superior photon statistics for spectral and timing analysis, its relatively wide point-spread function (PSF) does not allow us to search for the possible compact PWN around the pulsar.

In investigating the Chandra ACIS image, Weisskopf et al. (2011) briefly mentioned a possible feature associated with PSR J2021+4026  which may be indicative of a PWN. By fitting the X-ray image of PSR J2021+4026 with a circular Gaussian plus a constant background, the authors concluded that it is consistent with being point-like and placed the upper limit of its extent at ≲ 6''. Although this may be useful in quantifying the extent of bright and symmetric nebula (e.g., plerionic emission associated with young pulsars like Crab), such a method can overlook a faint asymmetric extended feature. This motivates us to reexamine this archival data to characterize the properties of this feature in detail. The adaptively smoothed ACIS image of a 05 × 05 field around PSR J2021+4026 is shown in Figure 3. It clearly shows a nebula-like structure that extends to west from PSR J2021+4026. For the PWN associated with a fast-moving pulsar, the extended X-ray emission is typically aligned with the direction of proper motion (see Gaensler 2005). Although the proper motion of PSR J2021+4026 is unknown, the orientation of its associated extended X-ray PWN indicates that it might be moving eastward. Assuming that the birth place of PSR J2021+4026 is not far away from the nominal center of G78.2+2.1 given by Green (2009; i.e., the white cross in Figure 1), we speculate the direction of the proper motion by extrapolating from the nominal remnant center to the current pulsar position. The speculated direction is illustrated by the white arrow in Figure 3. It is interesting to notice that it deviates from the symmetric axis of the nebula by only ∼5°.

Zoom In Zoom Out Reset image size

For quantifying the extent of such elongated compact nebula, we computed its brightness profile along its orientation (Hui & Becker 2007, 2008; Hui et al. 2012). We estimated the counts in 18 consecutive boxes with a size of 1'' × 10'' from the raw image with a pixel size of 05 × 05 along the extended feature (see the inset of Figure 4). The observed brightness profile is shown in Figure 4. To estimate the background level, we sampled the source-free regions around PSR J2021+4026 within a 2' × 2' FoV. The average background level and its 1σ uncertainties are indicated by the horizontal solid line and dotted lines, respectively. The nebular feature apparently extends for ∼10'' to the west before it falls to the estimated background level, which clearly exceeds the upper limit placed by Weisskopf et al. (2011) through a symmetric Gaussian fitting.

Zoom In Zoom Out Reset image size

To further examine the nature of its emission, we have extracted the photons from this feature within a box of 9' × 6' centered at R.A. = $20^{\rm h}21^{\rm m}30\buildrel{\mathrm{s}}\over{.}19$ decl. = $+40^{\circ }26^{^{\prime }}45{\buildrel{\prime\prime}\over{.}} 2$ (J2000). Only 25 counts are collected from this observation, which forbids any meaningful spectral analysis. However, we can still estimate its X-ray color and compare it with PSR J2021+4026. Following Trepl et al. (2010) and Weisskopf et al. (2011), we divide the energy range into three bands: soft (S: 0.5–1 keV), medium (M: 1–2 keV), and hard (H: 2–8 keV). The X-ray colors of this extended feature are estimated to be (H − S)/T = 0.41 ± 0.21 and M/T = 0.59 ± 0.16. In comparison to the X-ray colors of PSR J2021+4026, (H − S)/T = 0.02 ± 0.03, and M/T = 0.75 ± 0.02 (Weisskopf et al. 2011), the X-ray emission from the extended feature is apparently harder, which possibly indicates its non-thermal nature.

2.3. Analysis of the Pulsed X-Ray Emission from PSR J2021+4026

2.3.1. Pulsed Spectrum of PSR J2021+4026

From the phase-averaged spectral analysis, it has been shown that a single component model is not able to describe the observed data beyond ∼3 keV (Weisskopf et al. 2011; Lin et al. 2013). Statistically, an additional hard component is required at a confidence level of >99.995%. A BB+PL model fits the data reasonably well (see Figure 3 in Lin et al. 2013). Although the PWN as seen by Chandra is likely to contribute at least a part of this hard component, one cannot rule out the possibility that the observed non-thermal X-rays originate from the pulsar magnetosphere and thus have contributions to the observed pulsation. To distinguish these two scenarios, one has to determine in which rotational phases this PL component dominates. This motivates us to examine the pulsed spectrum of PSR J2021+4026 with the PN data.

We extracted the source spectrum in a circular region of a 20'' radius centered at the X-ray position of PSR J2021+4026 as determined by the source detection algorithm. The adopted extraction region corresponds to an encircled energy function of ∼76%. For determining the pulse phase of each photon, their arrival times were first corrected to the solar system barycenter with the aid of the XMMSAS task barycen by using the JPL DE405 earth ephemeris. For assigning the pulse phase to each event, we adopted the temporal parameters as determined by Lin et al. (2013).

We divided the rotational phases into two regimes, which are illustrated by the shaded regions in Figure 5. The phase intervals 0.2–0.7 and 0.85–1.2 are defined as the "on-pulse" and "off-pulse" components, respectively. Assuming the off-pulse component contributes a steady DC level across the entire rotational phase, the pulsed spectrum was then obtained by subtracting the off-pulse component from the on-pulse component. The response files were produced by the XMMSAS task rmfgen and arfgen. The spectrum is binned so as to have >50 counts per spectral bin. We used XSPEC 12.7.0 for all the spectral analysis reported in this work. The spectral fits were performed in 0.5–10 keV.

Zoom In Zoom Out Reset image size

The pulsed spectrum is found to be softer than the phase-averaged spectrum. Majority of the pulsed X-rays have energies <3 keV (see Figure 6). At energies >3 keV, ≳ 72% of the collected photons are contributed by the off-pulse component. We further found that the pulsed spectrum can be well described by a simple absorbed BB model (χ² = 12.86 for 15 degrees of freedom) without requiring additional components. The best-fit model yields a column density of $N_{H}=(9.1^{+5.2}_{-3.5})\times 10^{21}$ cm⁻², a temperature of $kT=0.23^{+0.06}_{-0.05}$ keV, a BB emitting region with a radius of $R=500^{+1053}_{-282}d_{2}$ m, where d₂ represents the distance to PSR J2021+4026 in units of 2 kpc. As the uncertainty of N_H is large, we fixed it at the value derived from the phase-averaged analysis (Lin et al. 2013), N_H = 7 × 10²¹ cm⁻², and constrained the BB temperature and radius to be $kT=0.26^{+0.03}_{-0.02}$ keV and $R=318^{+101}_{-77}d_{2}$ m, respectively (χ² = 13.51 for 16 degrees of freedom). For a conservative estimate, all the quoted errors of the spectral parameters are 1σ for two parameters of interest (i.e., Δχ² = 2.30 above the minimum).

Zoom In Zoom Out Reset image size

We have also attempted to perform a phase-resolved spectroscopy for each phase bin in Figure 5 in order to investigate how the spectral properties vary across the rotational phase. However, the photon statistic for individual phase bins is generally too small for a constraining analysis. The high instrumental background has further exacerbated the situation. Therefore, such analysis will not be considered further for this observation.

2.3.2. Multi-epoch X-Ray Spectral Analysis of PSR J2021+4026

As mentioned in the Introduction, PSR J2021+4026 is the first variable radio-quiet γ-ray pulsar where its γ-ray flux as spin-down properties suddenly change around MJD 55850 (2011 October 16). Allafort et al. (2013) argued that the variability of the γ-ray pulsed emission is due to certain global change in the magnetosphere. Since the X-ray emission from the hot polar cap results from the bombardment of the backflow current from the outergap (Cheng & Zhang 1999), it is not unreasonable to speculate that the change in the γ-ray properties might be accompanied with the change in the X-ray properties.

In the context of the outergap model (Cheng & Zhang 1999), the γ-ray luminosity is given by $L_{\gamma }\sim f^{3}\dot{E}$, where f is the fractional size of the outergap. Therefore, the observed change in the γ-ray luminosity, δL_γ/L_γ ∼ 18% (Allafort et al. 2013), implies that the gap size changed by δf/f ∼ δL_γ/3L_γ ∼ 6%. Since our analysis suggests that the pulsed X-ray emission is thermal, this can be produced through the polar cap heating by the return particle flux $\dot{N_{p}}=f\dot{N_{\rm GJ}}$, where $\dot{N_{\rm GJ}}$ is the Goldreich–Julian particle flux (Goldreich & Julian 1969). The X-ray luminosity can thus be estimated by $L_{X}\sim \dot{N_{p}}E_{p}$ where E_p is the typical particle energy deposited on the stellar surface. Therefore, the expected change in L_X is found to be $\delta L_{X}/L_{X}\sim \delta \dot{N_{p}}/\dot{N_{p}}\sim \delta f/f\sim 6\%$.

Since the Chandra (MJD 55435) and the XMM-Newton (MJD 56028) observations used in this study were performed before and after the γ-ray flux jump, it is instructive to examine if PSR J2021+4026 exhibited any X-ray variability. In order to investigate whether the X-ray spectral properties of PSR J2021+4026 vary, we examined the phase-averaged spectra of PSR J2021+4026 as obtained by Chandra ACIS-S3 and the EPIC camera on XMM-Newton. For Chandra, we extracted the source spectrum from a circular region with a radius of 2'' centered at the pulsar position. The background spectrum was sampled from an annular region with inner/outer radii of 25/4'' around the pulsar. For XMM-Newton, we followed the same procedure adopted by Lin et al. (2013) in preparing the phase-averaged spectrum.

We jointly fitted an absorbed BB+PL model to the Chandra and XMM-Newton spectra with the column absorption tied together throughout the analysis. In order to minimize the number of free parameters, we examined the variability of each spectral component one at a time. First, with the PL component in different epochs tied together, we allowed the BB component of the spectrum in different epochs to vary independently during the fitting. This yielded a column absorption of $N_{H}=6.4^{+0.8}_{-1.8}\times 10^{21}$ cm⁻², a photon index of Γ = 1.5 ± 0.8, and a PL normalization of $3.6^{+6.8}_{-2.6}\times 10^{-6}$ photons keV⁻¹ cm⁻² s⁻¹ at 1 keV. For the BB component, the best-fit temperature and emission radius for the epoch MJD 55435 are kT = 0.22 ± 0.04 keV and $R=350^{+359}_{-94}d_{2}$ m, respectively. After the γ-ray flux change, the BB parameters are found to be $kT=0.24^{+0.04}_{-0.02}$ keV and $R=288^{+193}_{-27}d_{2}$ m in the epoch MJD 56028.

For inspecting the possible variability of the PL component, the BB components of both spectra were tied. It yielded a column absorption of $N_{H}=6.4^{+0.8}_{-1.8}\times 10^{21}$ cm⁻², a BB temperature of kT = 0.24 ± 0.04 keV, and a emission radius of $R=298^{+22}_{-96}d_{2}$ m. The photon index and PL normalization before the γ-ray flux change are found to be $\Gamma =1.0^{+2.0}_{-1.0}$ and $1.4^{+11.1}_{-1.4} \times 10^{-6}$ photons keV⁻¹ cm⁻² s⁻¹ at 1 keV. The corresponding parameters after the γ-ray flux change are Γ = 1.8 ± 0.8 and $5.4^{+10.1}_{-3.7}\times 10^{-6}$ photons keV⁻¹ cm⁻² s⁻¹ at 1 keV.

Within the tolerance of the quoted statistical uncertainties, neither the pulsed thermal X-ray component nor the unpulsed non-thermal X-ray component is found to be variable in these two epochs.

2.3.3. Analysis of the X-Ray Pulse Profile of PSR J2021+4026

The X-ray pulse profile can also be used to investigate the global properties of the neutron star. As the polar cap sweeps across our line-of-sight, modulation in the soft X-ray regime can be seen (see Hui & Cheng 2004; Pechenick et al. 1983). Since the gravity of a neutron star is tremendous, the shape of the thermal X-ray pulse profile is determined by the near-field space-time curvature. The effect of gravity on the trajectory of emitted photons, which depends on the mass-to-radius ratio of the neutron star, must be considered in modeling the light curve.

Following Pechenick et al. (1983), we simulated the X-ray pulse profile resulting from the general relativistic calculation and compared it with the observational result. We chose our coordinates so that the observer is on the positive z-axis at r = r₀ where r₀ → ∞ (see Figure 6 in Hui & Cheng 2004). We described the stellar surface by angular spherical coordinates θ and ϕ where θ is measured from the z-axis defined above. For the photon emitted at an angle δ from the stellar surface, it will seem to the observer that they are emitted at an angle $\theta ^{^{\prime }}$ from the z-axis as a result of gravitational light bending. The relationship between θ and $\theta ^{^{\prime }}$ is given by

$$\begin{equation} \theta =\int _{0}^{\frac{GM}{Rc^{2}}} \left[\left(\frac{GM}{bc^{2}}\right)^{2}- (1-2u)u^{2}\right]^{-1/2}du, \end{equation} \tag{ 1 }$$

where $b=r_{0}\theta ^{^{\prime }}$ is the impact parameter of the photon and u = GM/c²r.

For a neutron star, (GM/c²R) must be less than 1/3. Therefore, a photon emitted from the surface that reaches the observer must have an impact parameter b ⩽ b_max where b_max = R(1–2GM/c²R)^−1/2 (Pechenick et al. 1983). The condition b = b_max sets the maximum value of θ, namely θ_max.

Considering a polar cap of an angular radius of α centered at θ = θ₀, a function h(θ; α, θ₀) is then defined as the range of ϕ included in the "one-dimensional slice" at θ of the polar cap (see Figure 6 in Hui & Cheng 2004). If θ₀ + α ⩽ θ_max ⩽ 180° and θ₀ − α ⩾ 0, then h(θ; α, θ₀) is defined as

$$\begin{equation} h(\theta; \alpha,\theta _{0})=\left\lbrace \begin{array}{@{}l} 2\cos ^{-1}\left(\frac{\cos \alpha -\cos \theta _{0}\cos \theta }{\sin \theta _{0}\sin \theta }\right) \\ [1pt] \qquad\qquad {\rm for}\quad \theta _{0}-\alpha \le \theta \le \theta _{0}+\alpha; \\ 0 \\ \qquad\qquad {\rm for}\quad \theta $ outside the range $\theta _{0}\pm \alpha. \end{array}\right. \end{equation} \tag{ 2 }$$

For generating the light curves, θ₀ is expressed as a function of time/rotational phase. Let β be the angle between the axis of the rotation of the star and the line joining the center of the polar cap and the center of the star, and γ is the angle between the axis of rotation and the z-axis, then

$$\begin{equation} \cos \theta _{0}=\sin \left(\beta \right)\sin \left(\gamma \right)\cos \left(\Omega t\right) + \cos \left(\beta \right)\cos \left(\gamma \right), \end{equation} \tag{ 3 }$$

where Ω is the rotational frequency of the star. The relative brightness can be expressed as a function of θ₀, M/R, and α

$$\begin{eqnarray} A\left(\theta _{0};M/R,\alpha \right)&=&\left(1-\frac{2GM}{c^{2}R}\right)^{2}\left(\frac{GM}{c^{2}R}\right)^{2}\nonumber\\ &&\times\int _0^{x_{{\rm max}}} h\left(x; \alpha,\theta _{0}\right)x\,dx, \end{eqnarray} \tag{ 4 }$$

where x = (c²b/GM) and x_max = (c²b_max/GM).

With a view toward minimizing the number of free parameters in the modeling, we utilized the results from another analysis. From the best-fit BB model of the pulsed spectrum, the polar cap size is found to be ∼320 m if PSR J2021+4026 locates at a distance of 2 kpc. Assuming a neutron star radius of R ∼ 10 km, we fixed the angular radius of the polar cap at α  ∼  2°. From modeling the γ-ray light curve, Trepl et al. (2010) suggest that the viewing angle can possibly be in a range of 83°–87°. For a given mass-to-radius ratio, (GM/c²R), the effect of varying the viewing angle in such a small range in the pulse profile is negligible. Also, for a polar cap with a small angular radius of ∼2°, it is likely that only one pole will cross the line-of-sight. This scenario is supported by the observed single broad peak. With these constraints, we minimized the number of free parameters by assuming a simple orthogonal rotator (γ = β = 90°) with a single pole contribution. This leaves the (GM/c²R) to be the only parameter for modeling the X-ray light curve. The best-fit model yields (GM/c²R) = 0.21 and a goodness-of-fit of χ² = 27.1 for 31 degrees of freedom. For R ∼ 10 km, it implies a neutron star mass of M ∼ 1.4 M_☉. The comparison of the best-fit model and the observed light curve is shown in Figure 7. For a conservative estimate, the 90% confidence interval for one parameter of interest (i.e., Δχ² = 2.71 above the minimum) is found to be 0.17  <  (GM/c²R)  <  0.25, which corresponds to M ∼ 1.2–1.7 M_☉ for R ∼ 10 km. Analysis with deeper follow-up observations can provide a tighter constraint on the mass-to-radius ratio of this neutron star.

Zoom In Zoom Out Reset image size

2.4. Imaging Spectroscopy of the Central and Southeastern Regions of G78.2+2.1

For investigating the X-ray emission from G78.2+2.1, we only focused on the extended features with relatively high surface brightness, which are highlighted by the solid-line ellipses in Figure 1 (referred as southeastern region hereafter) and Figure 2 (referred as central region hereafter). Before the spectra of these extended structures were extracted, contributions from all the resolved point sources were first subtracted from the data. The background spectra for the southeastern region and the central region were sampled from the nearby low count regions as illustrated by the dashed ellipse in Figure1 and dashed circle in Figure 2, respectively. The response files for this extended source analysis are generated by rmfgen and arfgen with uniform spatial averaging. The spectra obtained from different cameras were binned dynamically so as to achieve a comparable signal-to-noise ratio.

By inspecting the X-ray spectrum of the central region (see Figure 8), some emission line features such as Mg at ∼1.4 keV and Si at ∼1.9 keV can be clearly seen. This prompts us to examine the spectrum with an absorbed collision ionization equilibrium (CIE) plasma model (XSPEC model: VEQUIL). To examine whether the metal abundance of G78.2+2.1 deviates from the solar values, we thawed the corresponding parameters individually to see if the goodness of fit can be improved. The best-fit model yields a plasma temperature of kT ∼ 0.6 keV and a column absorption of N_H ∼ 10²² cm⁻². There is also an indication that Mg is overabundant in comparison to the solar values. However, such a model cannot provide an adequate description of the observed data even with metal abundance Mg open as a free parameter (χ² = 391.53 for 159 degrees of freedom). By examining the fitting residuals, systematic deviations at energies beyond ∼2 keV are noted. We suspected that an extra PL component might be required. With the additional PL component, the goodness-of-fit has been significantly improved (χ² = 209.36 for 157 degrees of freedom). It yields a column absorption of $N_{H}=8.2^{+0.8}_{-1.0}\times 10^{21}$ cm⁻², a plasma temperature of $kT=0.59^{+0.02}_{-0.03}$ keV, a Mg abundance of $1.9^{+0.4}_{-0.3}$ with respect to the solar value, an emission measure of ∫_Vn_en_HdV = (2.8 ± 0.9) × 10¹¹D² cm⁻³, and a PL index of $\Gamma =1.6^{+0.5}_{-0.2}$, where n_e, n_H, V, and D are the electron density (cm⁻³), hydrogen density (cm⁻³), volume of interest (cm³), and the source distance (cm).

Zoom In Zoom Out Reset image size

For further improving the spectral modeling, we examined the fitting residuals of the CIE+PL fit. We noticed that there is scattering of the residuals at energies greater than ∼2 keV, which probably stemmed from the additional PL. As demonstrated by Huang et al. (2014), the residuals in the hard band could possibly result from the residual soft proton contamination in individual cameras after the data screening. Therefore, instead of originating from the particle acceleration, the additional PL component merely provides a phenomenological description for such residual soft proton background with the PL index and normalization varying among different EPIC cameras. By disentangling the PL component in MOS1, MOS2, and PN, we found that the goodness-of-fit can be further improved (χ² = 175.04 for 153 degrees of freedom). In view of the different best-fit PL index inferred from different cameras, we conclude that the residuals in the hard band are contributed by the residual background. Under this consideration, the best-fit absorbed CIE component yields $N_{H}=7.9^{+0.8}_{-0.5}\times 10^{21}$ cm⁻², $kT=0.60^{+0.02}_{-0.03}$ keV, a Mg abundance of 2.0 ± 0.4 with respect to the solar value, and $\int _{V}n_{e}n_{H}dV=2.5^{+0.8}_{-0.5}\times 10^{11}D^{2}$ cm⁻³.

We have also examined the central part of G78.2+2.1 with a non-equilibrium ionization (NEI) model (XSPEC model: VNEI). With an additional PL component applied to account for the residual soft proton contamination in the individual camera, we found that the NEI model results in a further improved goodness-of-fit (χ² = 167.95 for 152 degrees of freedom). It yields a column absorption of $N_{H}=7.4^{+1.0}_{-1.4}\times 10^{21}$ cm⁻², a plasma temperature of $kT=1.6^{+0.6}_{-0.3}$ keV, a Mg abundance of 1.6 ± 0.2 with respect to the solar value, an emission measure of $\int _{V}n_{e}n_{H}dV=7.4^{+1.0}_{-0.8}\times 10^{10}D^{2}$ cm⁻³, and an ionization timescale of $n_{e}t=1.8^{+0.6}_{-0.4}\times 10^{10}$ s cm⁻³, where n_e and t are the electron density and time elapsed, respectively, since the gas has been shock-heated. We note that the inferred plasma temperature is significantly higher than the CIE fit. This can be due to the effect that the ionization states for a plasma in NEI at a given temperature are lower than those in the CIE situation (e.g., see Figure 11 in Vink 2012). Assuming a CIE condition can thus result in an underestimation of the plasma temperature. The NEI condition is further indicated by the best-fit ionization timescale, which is significantly less than that (n_et ∼ 10¹² s cm⁻³) required to reach CIE (see Vink 2012). In view of these, the NEI scenario is preferred.

For fitting the spectrum in the southeastern region (see Figure 9), we have also considered both CIE and NEI models. Same as the situation in the analysis of the diffuse X-rays from the central region, additional PL component were included for modeling the residual soft proton contamination in MOS1 and MOS2 individually. The CIE yields $N_{H}=4.4^{+1.0}_{-0.7}\times 10^{21}$ cm⁻², $kT=0.63^{+0.02}_{-0.03}$ keV, a Mg abundance of 3.3 ± 0.8 with respect to the solar value, and $\int _{V}n_{e}n_{H}dV=5.1^{+2.5}_{-1.1}\times 10^{11}D^{2}$ cm⁻³ (χ² = 305.11 for 242 degrees of freedom). In comparison, the NEI model results in an improved goodness-of-fit (χ² = 261.06 for 241 degrees of freedom). It yields N_H = (7.4 ± 0.6) × 10²¹ cm⁻², $kT=1.1^{+0.5}_{-0.3}$ keV, a Mg abundance of 1.2 ± 0.1 with respect to the solar value, $\tau =n_{e}t=2.4^{+1.9}_{-0.7}\times 10^{10}$ s cm⁻³, and $\int _{V}n_{e}n_{H}dV=8.2^{+3.4}_{-1.9}\times 10^{11}D^{2}$ cm⁻³. Both the goodness-of-fit and the small value of n_et resulting from this fit suggest that the remnant emission in this region is also in an NEI condition.

Zoom In Zoom Out Reset image size

We have reported the results from a detailed X-ray analysis of PSR J2021+4026 and G78.2+2.1. The column absorption deduced from the X-ray spectra of PSR J2021+4026 (see Section 2.3) is consistent with that deduced from various parts of the diffuse emission (Section 2.4). We also note that it is consistent with the neutral hydrogen density inferred from the H i absorption spectrum (Leahy et al. 2013). These results indicate that the pulsar emission, diffuse X-ray emission, and the radio shell are essentially at the same distance. Hence, the association between PSR J2021+4026 and G78.2+2.1 is supported by our investigation.

Given the pulsar–SNR association and assuming that the birth place of PSR J2021+4026 is not far away from the geometrical center of G78.2+2.1, we estimated the projected velocity of the pulsar. The angular separation between PSR J2021+4026 and the geometrical center is ∼01 (see Figure 1). At a distance of 2 kpc, this corresponds to a physical separation of 1.4 × 10¹⁴ km. Together with a Sedov age of ∼8000 yr deduced for G78.2+2.1 (Leahy et al. 2013), the magnitude of the projected velocity of PSR J2021+4026 is expected to be v_p ∼ 550 km s⁻¹, which is not unreasonable for the known pulsar population (Hobbs et al. 2005). The projected direction of the pulsar motion is indicated by the arrow in Figure 3. The speculated pulsar motion can possibly be checked by a dedicated γ-ray timing analysis of the full time-span Fermi-LAT data in further studies.

Such speculated pulsar velocity should be far exceeding the local speed of sound. For a pulsar moving supersonically, it is expected to drive a bow shock through the ambient medium. The pulsar wind particles will be accelerated and produce synchrotron X-ray emission. With the motion of the pulsar, this will result in a cometary-like nebula as the extended structure found in the high-resolution Chandra image (see Figures 3 and 4). In this case, the termination shock radius R_s is determined by the ram pressure balance between the relativistic pulsar wind particles and the circumstellar medium at the head of the shock:

$$\begin{equation} R_{s}\simeq \left(\frac{\dot{E}}{2\pi \rho v_{p}^{2}c}\right)^{1/2}\,{\sim}\,3\times 10^{16}\dot{E}_{34}^{1/2}n^{-1/2}v_{p,100}^{-1}\, {\rm cm}, \end{equation} \tag{ 5 }$$

where v_{p, 100} is the velocity of the pulsar in units of 100 km s⁻¹, $\dot{E}_{34}$ is the spin-down luminosity of the pulsar in units of 10³⁴ erg s⁻¹, and n is the number density of the circumstellar medium in units of cm⁻³.

For constraining n, we utilized the results inferred by the NEI plasma model fit of the central extended X-ray emission of G78.2+2.1  which surrounds PSR J2021+4026 (see Figure 2). The best-fit emission measure of this feature allows us to estimate the hydrogen density n_H and the electron density n_e in this circumstellar region. Assuming that n_e and n_H are uniform in the extraction region and a distance of 2 kpc, the emission measure can be approximated by n_en_HV  ∼  2.8 × 10⁵⁴ cm⁻³. We further assumed the geometry of an oblated spheroid for the spectral extraction region; the volume of interest is V  ∼  4 × 10⁵⁵ cm³. For a fully ionized plasma with ∼10% He (n_e  ∼  1.2n_H), n_H is estimated as ∼0.24 cm⁻³. Together with the spin-down power of $\dot{E}_{34}=10$ and our speculated v_{p, 100}  ∼  5.5, Equation (5) implies a termination radius of R_s  ∼  3.5 × 10¹⁶ cm. It corresponds to a stand-off angle of ∼1'' ahead of the pulsar at a distance of 2 kpc. Comparing this estimate to the angular size of the cometary-like feature behind the pulsar (see Figure 3), the ratio of termination shock radii between the directions immediately behind and directly ahead of the pulsar is estimated to be ∼10, which is comparable to the ratios observed in other fast-moving pulsars such as PSR J1747-2958 (Gaensler et al. 2004).

At 2 kpc, the physical size of the synchrotron X-ray nebula associated with PSR J2021+4026 is l_pwn  ∼  3 × 10¹⁷ cm. This implies that the timescale for the pulsar to traverse its nebula is τ_pwn = l_pwn/v_p  ∼  170 yr. The magnetic field strength of the nebula can be estimated by assuming that τ_pwn is comparable to the synchrotron cooling timescale of the electrons:

$$\begin{equation} \tau _{\rm syn}=\frac{6\pi m_{e}c}{\gamma \sigma _{T}B^{2}}\simeq 10^{5}\left(\frac{h\nu _{X}}{\rm keV}\right)^{-\frac{1}{2}}B_{\mu G}^{-\frac{3}{2}}\, {\rm yr}, \end{equation} \tag{ 6 }$$

where γ is the Lorentz factor of the wind, σ_T is the Thompson cross section, ν_X = 3γ²eB/2m_ec is the characteristic synchrotron frequency, and B_μG is the magnetic field in the shocked region in units of microgauss. This suggests that the nebular magnetic strength is at the order of ∼15 μG for a characteristic energy of hν_X ∼ 1 keV.

The magnetic field estimate can further enable us to compute the electron synchrotron cooling frequency ν_c:

$$\begin{equation} \nu _{\rm c}=\frac{18\pi em_{e}c}{\sigma _{T}^{2}\tau ^{2}_{\rm syn}B^{3}}, \end{equation} \tag{ 7 }$$

which is estimated to be ∼1.7 × 10¹⁹ Hz (i.e., hν_c  ∼  70 keV). Since this is far exceeding the observed frequencies, it suggests that the X-ray emission of the nebula is in a slow cooling regime (Chevalier 2000; Cheng et al. 2004). In this regime, electrons with the energy distribution, N(γ)∝γ^−p, are able to radiate their energy in the trail with the photon index Γ = (p + 1)/2. The index p due to shock acceleration typically lies between 2 and 3 (see Cheng et al. 2004 and references therein). This would result in a photon index ∼1.5–2.0. This is consistent with the observed photon index of Γ  ∼  1.5 for the unpulsed non-thermal spectral component of PSR J2021+4026.

For the X-ray pulsation of PSR J2021+4026, our analysis of the pulsed spectrum confirmed its thermal origin. We noted that the best-fit BB radius is comparable with the polar cap size of $\sqrt{({2\pi R}/{cP})}\,{\sim}\,300$ m by adopting a dipolar field geometry, a neutron star radius of R =10 km, and a rotational period of P = 265 ms. This suggests that the thermal emission originates from the hot polar cap. Such an assertion is further demonstrated by the agreement between the X-ray pulse profile and the simulated modulation by a hot spot on the surface of a canonical rotating neutron star (Figure 7).

The most remarkable properties of PSR J2021+4026 are the sudden changes in its spin-down rate, pulse profile, and flux in γ-ray on a timescale shorter than a week (Allafort et al. 2013). The authors speculated that such abrupt changes resulted from a shift in the magnetic field structure that in turn leads to the change of either magnetic inclination and/or the effective current (see discussion in Allafort et al. 2013). These can be precipitated by a reconfiguration of the magnetic field line footprints on the stellar surface. The polar cap size is defined by the footprint of the last open-field lines and its temperature is determined by the backflow current from the accelerating region. Therefore, according to the scenario proposed by Allafort et al. (2013), one should expect a correlated change in the thermal X-ray flux and/or the X-ray pulse profile.

We have attempted to look for such an expected X-ray change across the γ-ray jump by a joint analysis of Chandra and XMM-Newton data. Although we did not find any conclusive variability, we would like to point out that the significance of the analysis is limited by the small photon statistic of Chandra data. Also, the poor temporal resolution of this Chandra data does not allow any investigation of the X-ray pulsation. For a follow-up investigation of this unique pulsar, we encourage a long-term coordinated X-ray and γ-ray monitoring with XMM-Newton and Fermi, which could provide a better understanding the nature of its variability.

Apart from PSR J2021+4026, we have examined the diffuse X-ray emission of G78.2+2.1 in the FoV of our XMM-Newton observation (see Section 2.4). Leahy et al. (2013) have also analyzed the central region with Chandra data and obtained N_H = (7.5–11.1) × 10²¹ cm⁻², kT = 0.6–2.7 keV and n_et = (1.7–12) × 10¹⁰ s cm⁻³. Within the tolerance of 1σ uncertainties, our results are consistent with theirs. With the much improved photon statistic of our XMM-Newton data, we constrained the spectral parameters N_H, kT, and n_et to an accuracy of ∼32%, ∼56%, and ∼56%, respectively. For comparison, the uncertainties of the corresponding parameters reported by Leahy et al. (2013) are ∼39%, ∼203%, and ∼312%.

For the southeastern rim, Uchiyama et al. (2002) have investigated its X-ray properties with ASCA (R2 region in their work). They have modeled the spectrum with a CIE model and obtained a temperature of kT = 0.53 ± 0.07 keV, which is consistent with our CIE estimate. However, the quality of ASCA data did not allow the authors to discern whether the X-ray emission is in a CIE or an NEI state.

In our study, we confirmed that the remnant emission from our investigated regions are in an NEI state. This is probably due to the low electron density and the time elapsed since the gas has been shock-heated is not long enough for the plasma to reach the equilibrium. The best-fit emission measures and the ionization timescales allow us to estimate these quantities. From the above discussion, the electron density of the central region is found to be n_e  ∼  0.3 cm⁻³. With the best-fit ionization timescale of n_et  ∼  1.8 × 10¹⁰ s cm⁻³, the elapsed time since the arrival of the shock is estimated as t  ∼  1900 yr. For the diffuse emission in the southeastern region, the emission measure and the volume of interest are n_en_HV  ∼  3.1 × 10⁵⁵ cm⁻³ and V  ∼  10⁵⁷ cm⁻³, respectively. This implies an electron density of n_e  ∼  0.2 cm⁻³ in this region. Together with the ionization timescale of n_et  ∼  2.4 × 10¹⁰ s cm⁻³ inferred for this region, this suggests that the gas has been shock-heated ∼3800 yr ago. Such age estimates are significantly smaller than the Sedov age of G78.2+2.1 (Leahy et al. 2013). This might indicate that these plasma have been heated by the reverse shock(s) that returned to the remnant center not long ago. We would like to point out that this interpretation stems from the spectral fitting with a relatively simple NEI model. For example, the temperatures of different plasma constituents do not need to be equilibrated for such a small ionization timescale. A more sophisticated modeling of the observed remnant spectrum around the center can help to confirm whether or not the reverse shocks have arrived yet.

C.Y.H. is supported by the National Research Foundation of Korea through grant 2014R1A1A2058590. L.C.L. is supported by the Ministry of Science and Technology of Taiwan through grant NSC 101-2112-M-039-001-MY3. C.P.H. and Y.C. are supported by the Ministry of Science and Technology through the grant NSC 102-2112-M-008-020-MY3. J.T. and K.S.C. are supported by a 2014 GRF grant of the Hong Kong government under HKU 17300814P. A.K.H.K. is supported by the Ministry of Science and Technology of Taiwan through grants 100-2628-M-007-002-MY3, 100-2923-M-007-001-MY3, and 103-2628-M-007-003-MY3.

In this XMM-Newton observation, XMM J202154.7+402855 is detected by MOS1/2 (i.e., source 8 in Figure 1 and Table 1). In examining its light curve, we found that this source is significantly variable and resembles a flash-like event (see Figure 10). A rapid rise in intensity of XMM J202154.7+402855 occurred at ∼40 ks after the start of the investigation at ∼MJD 56028.31. Its count rate increased by a factor of ∼40 above the quiescent level with a timescale of ∼1 ks and became the brightest one among all the point sources detected in this observation. After reaching the peak, its count rate returned to the quiescent level in about an hour. In view of its interesting temporal behavior, we carefully examined timing and spectral properties of this newly detected flash-like event.

Zoom In Zoom Out Reset image size

In order to probe the spectral behavior of XMM J202154.7+402855, we divided its spectrum into two components: the quiescent spectrum (events in 0–35,000 s and 80,000–110,000 s) and the X-ray flash spectrum (events in 40,000–70,000 s).⁷ For the quiescent spectrum, we have examined it with various single component model: PL, BB, and Comptonization of soft photons in a hot plasma (Titarchuk 1994). None of this single component model can provide an acceptable description of the data. We proceeded to fit the quiescent spectrum with composite models. We found that the BB+PL model can fit the data reasonably well (χ² = 16.3 with 12 degrees of freedom), which yields a column density of $N_{H}=6.6^{+1.0}_{-1.1}\times 10^{21}$ cm⁻², a photon index of $\Gamma =1.8^{+1.3}_{-1.4}$, a PL normalization at 1 keV of $5.3^{+1.4}_{-1.5}\times 10^{-6}$ photons keV⁻¹ cm⁻² s⁻¹, a BB temperature of $kT=97.2^{+4.3}_{-5.5}$ eV and an emitting radius of $R=7.8^{+2.4}_{-2.9}d_{2}$ km. The unabsorbed flux of XMM J202154.7+402855 in the quiescent state is ∼9 × 10⁻¹³ erg cm⁻² s⁻¹ in 0.5–10 keV.

For investigating the X-ray flash spectrum, we further divided it into three stages: a rising stage (events in 40,000–41,400 s), a rapid declining stage (events in 41,400–45,000 s), and a slow declining stage (events in 45,000–70,000 s). For accounting the quiescent contribution in all these segments, we included the BB and PL components in the spectral fits with the parameters fixed at the best-fit values of the quiescent spectrum. In view of the narrow time-windows in dividing these stages, their photon statistics are lower than the quiescent spectrum. In order to minimize the number of free parameters, we also fixed the column absorption at N_H = 6.6 × 10²¹ cm⁻² as inferred from the quiescent spectrum. On top of the quiescent level, we have added an extra component for modeling the contribution from the flash-like event. We found that an additional Comptonized BB model (Nishimura et al. 1986; XSPEC model: COMPBB) is capable to yield reasonable spectral fits for all three stages. With the electron temperature of the plasma fixed at 50 keV, the best-fit temperature of the Comptonized BB for the rising stage, the rapid declining stage, and the slow declining stage are found to be $86.9^{+21.0}_{-24.1}$ eV, $103.5^{+13.4}_{-12.8}$ eV, and $95.9^{+15.3}_{-14.1}$ eV, respectively. The sum of the unabsorbed fluxes in all three stages is ∼4 × 10⁻¹¹ erg cm⁻² s⁻¹ in 0.5–10 keV.

Since both burst duration and spectral properties of XMM J202154.7+402855 are inconsistent with a Type-I X-ray burst, we proceeded to consider other possible emission scenarios. We have also explored its temporal behavior by fitting a PL model to the fading tail of the event. The results are summarized in Table 2. With all the parameters to be free, the light curve fitting yielded a PL index of −1.27 ± 0.16 and the best-fit function is shown as the red curve in Figure 10. Within its 3σ uncertainty, this value is consistent with that of a tidal disruption event (TDE) that has a time dependence of t^−5/3 (Rees 1988; Phinney 1989; Lodato 2012). With the PL index fixed at −5/3, the best-fit curve (i.e., blue curve in Figure 10) also results in a comparable goodness-of-fit.

Table 2. Best-fit Parameters with Both Free/Fixed PL Indices in $y=a_1+a_2 \times (t-a_3)^{a_4}$ to the Tail of X-Ray Flash

Parameters	a1	a2	a3	a4	$\chi ^2_{\nu }$
	(counts s−1)		(s)		(d.o.f)
a4 free	0.01	1511 ± 2074	40950 ± 215	−1.27 ± 0.16	1.099(62)
a4 fixed	0.01	49112 ± 5128	40388 ± 121	−5/3	1.117(63)

Download table as:  ASCIITypeset image

Considering the possibility of a TDE, we further attempted to search for the quasi-periodic oscillation (QPO) from the data, which was detected from Swift J1644+57 (Reis et al. 2012). A ∼200 s QPO was detected from Swift J1644+57, which is interpreted as the Keplerian frequency of the innermost stable circular orbit of a supermassive black hole (Reis et al. 2012). Motivated by this discovery, we searched for the periodic signal from XMM J202154.7+402855 in a range of 0–400 s centered at 200 s with a resolution of 0.1 s by χ² test. The highest peak obtained from the periodogram is at 115.1 s with $\chi ^2_{15}$ less than 3. Even we only considered those events obtained between the main outburst of 40,000–55,000 s shown in Figure 10, the similar result of the $\chi ^2_{15}=3.2$ was detected at a trial period of 115 s. We also considered detecting the periodic signal using the Lomb–Scargle method (Scargle 1982) on light curves rebinned with 20 s and 40 s. With this method, we concluded that there is no periodicity that can be detected with a power more significant than 90% confidence level. Our analysis indicates that there is no stable periodic signal can be detected from the current observation.

Since QPO might appear intermittently or varies with time, these make the aforementioned periodicity search for the whole light curve difficult. In view of this, we have also searched the possible periodic signal by computing the dynamic power spectrum (Clarkson et al. 2003a, 2003b). We adopted a window size of 1000 s, which is approximately the duration of the flash-like event. In order to depress the effect of the trend, we used the empirical mode decomposition (Huang et al. 1998) to filter the trend and only the de-trended light curve was examined by the dynamical power spectrum. However, there is no significant signal of QPO can be detected from the dynamic Lomb–Scargle periodogram. Together with the non-detection of any periodic signal, the fact that the burst duration of XMM J202154.7+402855 is far shorter than a typical TDE, which lasts for a timescale of months (e.g., Reis et al. 2012) does not favor this scenario.

Another possible source nature of XMM J202154.7+402855 is a flaring early-type star. This requires a search for the optical counterpart for constraining its properties. Utilizing the USNO-B1.0 catalog (Monet et al. 2003), we have identified a bright source, USNO-B1.0 1304-0388936, locates at ∼3'' away from the nominal X-ray position of XMM J202154.7+402855 with magnitudes of B = 14.48, R = 13.65, and I = 12.85. Assuming that it is the optical counterpart of XMM J202154.7+402855, its X-ray-to-optical flux ratio is f_x/f_opt  ∼  10⁻³ during quiescence, which is quite typical for a field star (Maccacaro et al. 1988). To investigate if the positional offset between the X-ray source and the optical counterpart is the result of proper motion, we have also checked the UCAC3 catalog (Zacharias et al. 2010). The source 3UC 261-199420 in UCAC3 has its nominal position differed from USNO-B1.0 1304-0388936 only by ∼05 and its proper motion is very small: μ_RA = 1.9 mas yr⁻¹ and μ_Dec = −0.9 mas yr⁻¹. Calculating the angular shift from the central epoch for the position given by UCAC3 to the epoch of our XMM-Newton observation, we found that it only moves by ∼31 mas. Hence, the offset between XMM J202154.7+402855 and the optical source cannot be reconciled by the proper motion. Although the positional offset is comparable to the absolute point accuracy of XMM-Newton,⁸ further investigation is required to secure the association between XMM J202154.7+402855 and the optical source.