X-ray Studies of the Extended TeV Gamma-Ray Source VER J2019+368

TL;DR

This study uses X-ray observations to analyze the extended PWN associated with VER J2019+368, revealing its properties, constraining electron transport, and modeling its contribution to TeV gamma-ray emission.

Contribution

First detailed X-ray analysis of VER J2019+368's PWN, constraining its distance, magnetic field, and electron transport, and modeling its gamma-ray emission contribution.

Findings

X-ray PWN extent is up to 15' x 10' with consistent spectral index.

X-ray absorption suggests the source is closer than 10 kpc.

Synchrotron and inverse Compton models explain about 80% of TeV flux.

Abstract

This article reports the results of X-ray studies of the extended TeV $\gamma$-ray source VER J2019+368. Suzaku observations conducted to examine properties of the X-ray pulsar wind nebula (PWN) around PSR J2021+3651 revealed that the western region of the X-ray PWN has a source extent of $15' \times 10'$ with the major axis oriented to that of the TeV emission. The PWN-west spectrum was closely fitted by a power-law for absorption at $N({\rm H}) = (8.2^{+1.3}_{-1.1}) \times 10^{21}~{\rm cm^{-2}}$ and a photon index of $\Gamma = 2.05\pm0.12$, with no obvious change in the index within the X-ray PWN. The measured X-ray absorption indicates that the distance to the source is much less than $10~{\rm kpc}$ inferred by radio data. Aside from the PWN, no extended emission was observed around PSR J2021+3651 even by Suzaku. Archival data from the XMM-Newton were also analyzed to complement the Suzaku observations, indicating that the eastern region of the X-ray PWN has a similar spectrum ($N(\rm H)=(7.5 \pm 0.9) \times 10^{21}~{\rm cm^{-2}}$ and $\Gamma=2.03 \pm 0.10$) and source extent up to at least $12'$ along the major axis. The lack of significant change in the photon index and the source extent in X-ray are used to constrain the advection velocity or the diffusion coefficient for accelerated X-ray-producing electrons. A mean magnetic field of ${\sim}3~\mu{\rm G}$ is required to account for the measured X-ray spectrum and reported TeV $\gamma$-ray spectrum. A model calculation of synchrotron radiation and inverse Compton scattering was able to explain ${\sim}80\%$ of the reported TeV flux, indicating that the X-ray PWN is a major contributor of VER J2019+368.

Tables

Table 1. Table 1: Summary of observations

Observatory	Region	Pointingaa Position of the center of the XIS (Suzaku) or MOS (XMM-Newton) FOV.		Observation date	Net exposure
		RA (deg)	DEC (deg)		(ks)
Suzaku	S1	305.064	36.873	2014 Nov. 09	35.0
	S2	304.792	36.828	2014 Nov. 10	35.7
XMM-Newton	X1	305.273	36.851	2012 Apr. 07	83.4

Table 2. Table 2: Summary of the spectral fits of the PWN-west by Suzaku -XIS

Regiona	$N\left(H\right)$	$\Gamma$	$f\left(0.5\text{–}2\mathrm{keV}\right)$	$f\left(2\text{–}10\mathrm{keV}\right)$	${\chi }^2$/DOF
	(${10}^{21}{\mathrm{cm}}^{-2}$)		(${10}^{-13}\mathrm{erg}{s}^{-1}{\mathrm{cm}}^{-2}$)	(${10}^{-13}\mathrm{erg}{s}^{-1}{\mathrm{cm}}^{-2}$)
all($0^{\prime }\text{–}{15}^{\prime }$)	${8.2}_{-1.1}^{+1.3}$	$2.05\pm 0.12$	${6.04}_{-0.40}^{+0.42}$	${26.1}_{-1.8}^{+1.6}$	211.1/188
$0^{\prime }\text{–}{3}^{\prime }$	${8.2}_{-2.0}^{+2.1}$	$2.07\pm 0.21$	${1.97}_{-0.22}^{+0.16}$	${8.27}_{-1.06}^{+0.89}$	65.5/61
$3^{\prime }\text{–}{6}^{\prime }$	${6.2}_{-1.5}^{+1.7}$	$1.96\pm 0.18$	${1.68}_{-0.20}^{+0.17}$	${6.52}_{-0.73}^{+0.76}$	60.5/67
$6^{\prime }\text{–}{9}^{\prime }$	${7.1}_{-1.5}^{+1.8}$	$2.06\pm 0.18$	$1.22\pm 0.13$	${4.71}_{-0.52}^{+0.51}$	128.9/99
$9^{\prime }\text{–}{12}^{\prime }$	${12.5}_{-3.2}^{+3.9}$	$2.30\pm 0.32$	${0.57}_{-0.12}^{+0.07}$	${2.77}_{-0.53}^{+0.39}$	85.3/73
${12}^{\prime }\text{–}{15}^{\prime }$	${13.7}_{-4.4}^{+5.8}$	$2.29\pm 0.42$	${0.30}_{-0.10}^{+0.07}$	${1.69}_{-0.61}^{+0.35}$	54.7/52
$0^{\prime }\text{–}{3}^{\prime }$	8.2(fixed)	$2.07\pm 0.10$	${1.96}_{-0.18}^{+0.20}$	${8.26}_{-0.90}^{+0.86}$	65.5/62
$3^{\prime }\text{–}{6}^{\prime }$	8.2(fixed)	$2.14\pm 0.11$	${1.61}_{-0.18}^{+0.15}$	${6.28}_{-0.61}^{+0.63}$	64.5/68
$6^{\prime }\text{–}{9}^{\prime }$	8.2(fixed)	$2.15\pm 0.10$	${1.20}_{-0.12}^{+0.11}$	${4.62}_{-0.46}^{+0.51}$	130.0/100
$9^{\prime }\text{–}{12}^{\prime }$	8.2(fixed)	$1.94\pm 0.14$	${0.60}_{-0.09}^{+0.08}$	${2.93}_{-0.46}^{+0.47}$	90.5/74
${12}^{\prime }\text{–}{15}^{\prime }$	8.2(fixed)	$1.84\pm 0.19$	${0.33}_{-0.08}^{+0.07}$	${1.80}_{-0.36}^{+0.71}$	59.3/53

Table 3. Table 3: Summary of the spectral fits of the PWN-west, PWN-east, and Arc by XMM-Newton MOS

Region	$N\left(H\right)$	$\Gamma$	$f\left(0.5\text{–}2\mathrm{keV}\right)$	$f\left(2\text{–}10\mathrm{keV}\right)$	${\chi }^2$/DOF
	${10}^{21}{\mathrm{cm}}^{-2}$		${10}^{-13}\mathrm{erg}{s}^{-1}{\mathrm{cm}}^{-2}$	${10}^{-13}\mathrm{erg}{s}^{-1}{\mathrm{cm}}^{-2}$
west	$8.1\pm 1.0$	$2.10\pm 0.12$	${4.77}_{-0.30}^{+0.22}$	${19.3}_{-1.5}^{+0.9}$	271.8/225
east	$7.5\pm 0.9$	$2.03\pm 0.10$	${5.01}_{-0.26}^{+0.25}$	${20.4}_{-1.3}^{+1.1}$	304.4/244
arc	${9.3}_{-1.4}^{+1.5}$	$2.13\pm 0.18$	${1.66}_{-0.13}^{+0.12}$	${7.31}_{-0.79}^{+0.67}$	138.8/116

Table 4. Table 4: Summary of the spectral fits of the PWN-west and PWN-east divided into segments

Region		$N\left(H\right)$	$\Gamma$	$f\left(0.5\text{–}2\mathrm{keV}\right)$	$f\left(2\text{–}10\mathrm{keV}\right)$	${\chi }^2$/DOF
		${10}^{21}{\mathrm{cm}}^{-2}$		${10}^{-13}\mathrm{erg}{s}^{-1}{\mathrm{cm}}^{-2}$	${10}^{-13}\mathrm{erg}{s}^{-1}{\mathrm{cm}}^{-2}$
west	$0^{\prime }\text{–}{3}^{\prime }$	$7.8\pm 1.1$	${1.99}_{-0.14}^{+0.13}$	$1.71\pm 0.09$	${7.60}_{-0.53}^{+0.50}$	123.0/104
	$3^{\prime }\text{–}{6}^{\prime }$	${8.1}_{-1.3}^{+1.5}$	${2.07}_{-0.18}^{+0.19}$	$1.33\pm 0.10$	${5.57}_{-0.52}^{+0.47}$	82.6/80
	$6^{\prime }\text{–}{9}^{\prime }$	${8.3}_{-2.0}^{+1.7}$	${2.44}_{-0.33}^{+0.27}$	${1.23}_{-0.17}^{+0.09}$	${3.52}_{-0.53}^{+0.65}$	110.5/73
	$9^{\prime }\text{–}{12}^{\prime }$	${10.8}_{-3.9}^{+4.7}$	${2.53}_{-0.59}^{+0.73}$	${0.61}_{-0.20}^{+0.10}$	${2.00}_{-0.85}^{+0.41}$	79.3/57
	$0^{\prime }\text{–}{3}^{\prime }$	8.1(fix)	$2.02\pm 0.08$	${1.70}_{-0.08}^{+0.10}$	${7.54}_{-0.43}^{+0.52}$	123.3/105
	$3^{\prime }\text{–}{6}^{\prime }$	8.1(fix)	$2.07\pm 0.10$	${1.33}_{-0.07}^{+0.08}$	${5.58}_{-0.47}^{+0.48}$	82.6/81
	$6^{\prime }\text{–}{9}^{\prime }$	8.1(fix)	$2.40\pm 0.16$	$1.24\pm 0.10$	${3.56}_{-0.46}^{+0.53}$	110.4/74
	$9^{\prime }\text{–}{12}^{\prime }$	8.1(fix)	$2.18\pm 0.30$	${0.61}_{-0.11}^{+0.12}$	${2.25}_{-0.58}^{+0.64}$	80.5/58
east	$0^{\prime }\text{–}{3}^{\prime }$	${8.3}_{-1.2}^{+1.4}$	$1.99\pm 0.15$	${1.30}_{-0.09}^{+0.08}$	${6.11}_{-0.40}^{+0.44}$	109.1/96
	$3^{\prime }\text{–}{6}^{\prime }$	${7.0}_{-1.2}^{+1.3}$	$2.11\pm 0.17$	$1.46\pm 0.10$	${5.21}_{-0.45}^{+0.51}$	106.1/92
	$6^{\prime }\text{–}{9}^{\prime }$	${7.8}_{-1.6}^{+1.8}$	${2.00}_{-0.21}^{+0.22}$	${1.08}_{-0.10}^{+0.07}$	${4.74}_{-0.54}^{+0.44}$	75.2/74
	$9^{\prime }\text{–}{12}^{\prime }$	${7.2}_{-1.8}^{+2.1}$	${1.97}_{-0.25}^{+0.27}$	${1.11}_{-0.14}^{+0.10}$	${4.75}_{-0.58}^{+0.76}$	75.5/72
	$0^{\prime }\text{–}{3}^{\prime }$	7.5(fix)	${1.91}_{-0.09}^{+0.08}$	$1.31\pm 0.09$	${6.22}_{-0.41}^{+0.47}$	110.4/97
	$3^{\prime }\text{–}{6}^{\prime }$	7.5(fix)	$2.17\pm 0.10$	$1.46\pm 0.10$	${5.13}_{-0.40}^{+0.39}$	106.6/93
	$6^{\prime }\text{–}{9}^{\prime }$	7.5(fix)	$1.97\pm 0.12$	${1.08}_{-0.08}^{+0.09}$	${4.78}_{-0.40}^{+0.53}$	75.3/75
	$9^{\prime }\text{–}{12}^{\prime }$	7.5(fix)	$2.00\pm 0.15$	${1.11}_{-0.11}^{+0.10}$	${4.71}_{-0.60}^{+0.51}$	75.6/73

Table 5. Table 5: Summary of the spectral fit of the entire FOV data

Parameter	Value
$N{\left(H\right)}_{\mathrm{CXB}}\left({10}^{21}{\mathrm{cm}}^{-2}\right)$	30(fixed)
${\Gamma }_{\mathrm{CXB}}$	1.46(fixed)
${\mathrm{Norm}}_{\mathrm{CXB}}$	$5.79\times {10}^{-4}$(fixed)
$N{\left(H\right)}_{\mathrm{high}}\left({10}^{21}{\mathrm{cm}}^{-2}\right)$	30(fixed)
$k{T}_{\mathrm{high}}\left(\mathrm{keV}\right)$	2.5(fixed)
$A_{\mathrm{high}}\left(Z_{☉}\right)$	0.3(fixed)
${\mathrm{EM}}_{\mathrm{high}}$	$2.5\times {10}^{-3}\left(\le 9.3\times {10}^{-3}\right)$
$N{\left(H\right)}_{\mathrm{in}}\left({10}^{21}{\mathrm{cm}}^{-2}\right)$	6.7(fixed)
${\Gamma }_{\mathrm{in}}$	1.45(fixed)
${\mathrm{Norm}}_{\mathrm{in}}$	$6.0\times {10}^{-5}$(fixed)
$N{\left(H\right)}_{\mathrm{out}}\left({10}^{21}{\mathrm{cm}}^{-2}\right)$	6.7(fixed)
${\Gamma }_{\mathrm{out}}$	1.82(fixed)
${\mathrm{Norm}}_{\mathrm{out}}$	$2.15\times {10}^{-4}$(fixed)
$N{\left(H\right)}_{\mathrm{src}}\left({10}^{21}{\mathrm{cm}}^{-2}\right)$	$1.5\left(\le 40.7\right)$
${\Gamma }_{\mathrm{src}}$	${1.84}_{-0.16}^{+0.87}$
${\mathrm{Norm}}_{\mathrm{src}}$	$\left({4.1}_{-1.7}^{+11.8}\right)\times {10}^{-3}$
${\Gamma }_{\mathrm{SP}}$	${0.21}_{-0.51}^{+1.16}$
${\mathrm{Norm}}_{\mathrm{SP}}$	$\left({1.1}_{-1.0}^{+9.0}\right)\times {10}^{-2}$
${\chi }^2$/DOF	174.1/142

Table 6. Table 6: Summary of the spectral fits of the background region obtained through a joint-fit with the PWN-west, PWN-east, and Arc

	PWN-west	PWN-east	Arc
$N{\left(H\right)}_{\mathrm{CXB}}\left({10}^{21}{\mathrm{cm}}^{-2}\right)$	30(fixed)	30(fixed)	30(fixed)
${\Gamma }_{\mathrm{CXB}}$	1.46(fixed)	1.46(fixed)	1.46(fixed)
${\mathrm{Norm}}_{\mathrm{CXB}}$	$4.87\times {10}^{-5}$(fixed)	$4.83\times {10}^{-5}$(fixed)	$1.99\times {10}^{-5}$(fixed)
$N{\left(H\right)}_{\mathrm{high}}\left({10}^{21}{\mathrm{cm}}^{-2}\right)$	30(fixed)	30(fixed)	30(fixed)
$k{T}_{\mathrm{high}}\left(\mathrm{keV}\right)$	2.5(fixed)	2.5(fixed)	2.5(fixed)
$A_{\mathrm{high}}\left(Z_{☉}\right)$	0.3(fixed)	0.3(fixed)	0.3(fixed)
${\mathrm{EM}}_{\mathrm{high}}$	$\left(6.4\pm 2.3\right)\times {10}^{-4}$	$\left(8.3\pm 2.2\right)\times {10}^{-4}$	$\left({2.7}_{-1.0}^{+0.9}\right)\times {10}^{-4}$
$N{\left(H\right)}_{\mathrm{mid}}\left({10}^{21}{\mathrm{cm}}^{-2}\right)$	${7.09}_{-0.74}^{+0.72}$	${6.19}_{-0.74}^{+0.70}$	${7.48}_{-0.82}^{+0.79}$
$k{T}_{\mathrm{mid}}\left(\mathrm{keV}\right)$	${0.643}_{-0.074}^{+0.072}$	${0.643}_{-0.041}^{+0.078}$	${0.639}_{-0.058}^{+0.069}$
$A_{\mathrm{mid}}\left(Z_{☉}\right)$	0.3(fixed)	0.3(fixed)	0.3(fixed)
${\mathrm{EM}}_{\mathrm{mid}}$	$\left({1.40}_{-0.31}^{+0.35}\right)\times {10}^{-3}$	$\left({1.13}_{-0.26}^{+0.30}\right)\times {10}^{-3}$	$\left({6.5}_{-1.6}^{+1.9}\right)\times {10}^{-4}$
$N{\left(H\right)}_{\mathrm{low}}\left({10}^{22}{\mathrm{cm}}^{-2}\right)$	0(fixed)	0(fixed)	0(fixed)
$k{T}_{\mathrm{low}}\left(\mathrm{keV}\right)$	0.1(fixed)	0.1(fixed)	0.1(fixed)
$A_{\mathrm{low}}\left(Z_{☉}\right)$	1.0(fixed)	1.0(fixed)	1.0(fixed)
${\mathrm{EM}}_{\mathrm{low}}$	$\left(2.56\pm 0.39\right)\times {10}^{-4}$	$\left(2.43\pm 0.38\right)\times {10}^{-4}$	$\left(8.7\pm 2.7\right)\times {10}^{-5}$
$E_1\left(\mathrm{keV}\right)$	1.49(fixed)	1.49(fixed)	1.49(fixed)
${\mathrm{Norm}}_1$	$\left(4.08\pm 0.29\right)\times {10}^{-5}$	$\left(3.83\pm 0.27\right)\times {10}^{-5}$	$\left(1.63\pm 0.17\right)\times {10}^{-5}$
$E_2\left(\mathrm{keV}\right)$	1.75(fixed)	1.75(fixed)	1.75(fixed)
${\mathrm{Norm}}_2$	$\left({6.3}_{-2.2}^{+2.3}\right)\times {10}^{-6}$	$\left(1.25\pm 0.23\right)\times {10}^{-5}$	$\left(4.5\pm 1.4\right)\times {10}^{-6}$
$E_3\left(\mathrm{keV}\right)$	0.65(fixed)	0.65(fixed)	0.65(fixed)
${\mathrm{Norm}}_3$	$\left(1.80\pm 0.43\right)\times {10}^{-5}$	$\left(1.68\pm 0.41\right)\times {10}^{-5}$	$\left({7.0}_{-1.6}^{+1.7}\right)\times {10}^{-6}$
${\Gamma }_{\mathrm{SP},\mathrm{bg}}$	0.21(fixed)	0.21(fixed)	0.21(fixed)
${\mathrm{Norm}}_{\mathrm{SP},\mathrm{bg}}$	$3.33\times {10}^{-4}$(fixed)	$3.33\times {10}^{-4}$(fixed)	$3.33\times {10}^{-4}$(fixed)
${\Gamma }_{\mathrm{SP},\mathrm{src}}$	0.21(fixed)	0.21(fixed)	0.21(fixed)
${\mathrm{Norm}}_{\mathrm{SP},\mathrm{src}}$	$1.11\times {10}^{-3}$(fixed)	$1.07\times {10}^{-3}$(fixed)	$4.51\times {10}^{-3}$(fixed)
${\chi }^2$/DOF	271.8/225	304.4/244	138.8/116

Equations

$$E_e\simeq 132\mathrm{TeV}{\left(\frac{{ϵ}_{\mathrm{X}}}{1\mathrm{keV}}\right)}^{0.5}{\left(\frac{B}{3\mu \mathrm{G}}\right)}^{-0.5}.$$

$$E_e\simeq 17.2\mathrm{TeV}{\left(\frac{{ϵ}_{\gamma }}{1\mathrm{TeV}}\right)}^{0.5}.$$

$$E_{\mathrm{bk}}\simeq 80\mathrm{TeV}{\left(\frac{t_0}{17.2\mathrm{kyr}}\right)}^{-1}{\left(\frac{B}{3\mu \mathrm{G}}\right)}^{-2},$$

$$\tau ({ϵ}_{\mathrm{X}})\simeq 10.5\mathrm{kyr}{\left(\frac{{ϵ}_{\mathrm{X}}}{1\mathrm{keV}}\right)}^{-0.5}{\left(\frac{B}{3\mu \mathrm{G}}\right)}^{-1.5}.$$

$${\beta }_{\mathrm{adv}}\ge 7.9\times 10^{-3}{\left(\frac{B}{3\mu \mathrm{G}}\right)}^{1.5}\left(\frac{d}{1.8\mathrm{kpc}}\right).$$

$$D\ge 2.9\times 10^{27}\mathrm{c}{\mathrm{m}}^2{\mathrm{s}}^{-1}{\left(\frac{B}{3\mu \mathrm{G}}\right)}^{1.5}{\left(\frac{d}{1.8\mathrm{kpc}}\right)}^2.$$

$$D$$

$$\eta \simeq 0.60{\left(\frac{B}{3\mu \mathrm{G}}\right)}^3{\left(\frac{d}{1.8\mathrm{kpc}}\right)}^2.$$

Full text

pulsars: general — cosmic rays — gamma rays: observations — X-rays: ISM

\AuthorCallLimit

=1 \fullcollaborationNameThe Friends of AASTeX Collaboration

1 Introduction

Star-forming regions host several possible cosmic-ray (CR) accelerators such as supernova remnants (SNRs), pulsars and pulsar wind nebulae (PWNs), Wolf–Rayet stars, and OB associations. Cygnus-X (Piddington & Minnett, 1952; Uyaniker et al., 2001) is one of such nearby star-forming region; it is located at approximately 1.5 kpc (Rygl et al., 2012) and has long been studied at various wavebands, although care must be taken to properly associate individual sources to Cygnus-X, as there are several spiral arms in the same direction. A survey of the Northern Hemisphere sky by the Milagro Gamma-Ray Observatory identified several bright and extended TeV $\gamma$-ray sources (Abdo et al., 2007). MGRO J2019+37 is the brightest Milagro source in the direction of Cygnus-X with a measured flux of approximately $80\%$ of the Crab Nebula flux at 20 TeV. Despite extensive studies at various wavebands, the nature of MGRO J2019+37 remained unsettled because of its large source extent ($\sigma=0\fdg 7$ when modeled with a two-dimensional Gaussian probability density function; Abdo et al., 2012). The imaging atmospheric Cherenkov telescope array VERITAS carried out a deep observation of the MGRO J2019+37 region and resolved it into two sources. The brighter source, VER J2019+368, is an extended source that accounts for the bulk of MGRO J2019+37 in terms of morphology and spectrum (Aliu et al., 2014). Its peak is located at right ascension RA (J2000) ${\rm 20^{h}19^{m}25^{s}}$ and declination DEC (J2000) ${\rm 36\arcdeg 48^{\prime}14^{\prime\prime}}$, and its angular extension is estimated to be ${\sim}0\fdg 34$ and ${\sim}0\fdg 13$ along its major and minor axis respectively, with the orientation of the major axis $71\arcdeg$ east of north. Its TeV spectrum is hard and represented by a single power-law with a photon index $\Gamma\sim 1.75$ and an integrated energy flux at 1–10 TeV of ${\sim}6.7\times 10^{-12}~{}{\rm erg~{}s^{-1}~{}cm^{-2}}$. The emission region contains the energetic pulsar PSR J2021+3651 111 The pulse period and its derivatives are $104~{}{\rm ms}$ and $9.57\times 10^{-14}$, respectively (Abdo et al., 2009), giving the surface magnetic field of $3.2\times 10^{12}~{}{\rm G}$.

with the characteristic age of 17.2 kyr and spin-down luminosity of $3.4\times 10^{36}~{}{\rm erg~{}s^{-1}}$ (Roberts et al., 2002), its PWN, the H II region Sharpless 104 (Sh 2-104), and the Wolf–Rayet star WR 141, which are potential counterparts of the observed TeV emission. Fermi-LAT (Abdo et al., 2009) detected emissions from PSR J2021+3651 and gave an upper limit on the extended GeV $\gamma$-ray emissions around it, consistent with the hard spectrum of VER J2019+368 in the TeV band. Although the distance to the pulsar is inferred from the radio data to be ${\geq}10~{}{\rm kpc}$ (Roberts et al., 2002), this conclusion remains controversial (e.g., Van Etten et al., 2008) and does not coincide with a detailed study of the distance by Kirichenko et al. (2015).

Among the source classes of possible counterparts mentioned above, only PWNs are the established extended TeV $\gamma$-ray sources. Nevertheless, the association of the X-ray PWN (named G75.2+0.1; Hessels et al., 2004) to VER J2019+368 is a matter of debate, as its position is offset from the peak in TeV by ${\sim}20^{\prime}$ and its reported extent is $\leq 15^{\prime}$ (Roberts et al., 2008), which is much smaller than the size of the TeV $\gamma$-ray emission region. To resolve this, we carried out deep X-ray observations using the X-ray Imaging Spectrometer (XIS) (Koyama et al., 2007) on board the Suzaku satellite (Mitsuda et al., 2007), which is very sensitive to extended X-ray emission. We aimed to accurately measure the spectral and morphological properties of the X-ray PWN and to observe unknown extended X-ray emissions in the region of VER J2019+368. We also analyzed archival XMM-Newton (Jansen et al., 2001) data in order to complement the Suzaku-XIS observations, which did not cover the entire PWN. This paper is organized as follows. We describe the X-ray observations and our data reduction in Section 2. The results of the data analysis are presented in Section 3, in which we provide the detailed spectral and morphological properties of the X-ray PWN to the west of the pulsar and make a comparison between the eastern and western regions of the PWN. Discussion of the PWN’s association with VER J2019+368 based on its X-ray properties and a multiwavelength spectrum is provided in Section 4. A summary of this study and future prospects are presented in Section 5.

2 Observations and Data Reduction

In November 2014, we carried out deep X-ray observations of VER J2019+368 region using Suzaku-XIS. In order to constrain the X-ray properties of the PWN around PSR J2021+3651 and search for unknown extended X-ray emissions, we conducted two observations. As shown in Figure 1, these covered the main region of the TeV emission. The objective of the first pointing (S1) was to characterize the X-ray properties of the western region of the PWN, while the second pointing (S2) had the objective of searching for unknown extended X-ray emissions.

The observations were carried out using the XIS on the focal plane of the X-Ray Telescope (XRT; Serlemitsons et al., 2007) on Suzaku. The XIS consists of two front-illuminated (FI) X-ray charge coupling devices (CCDs) (XIS0 and 3) 222Because of an anomaly that occurred in November 2006, the operation of another FI sensor, XIS2, has been terminated. and one backside-illuminated (BI) X-ray CCD (XIS1). The combined XIS and XRT system are sensitive within the energy range of 0.3–12 keV. Although its angular resolution is moderate (half-power diameter $\sim 2^{\prime}$), the XIS+XRT system provides a low and stable instrumental background (Mitsuda et al., 2007; Tawa et al., 2008) and is, therefore, suitable for the detailed study of extended emissions with low surface brightness. Data were analyzed using the HEASOFT 333http://heasarc.nasa.gov/lheasoft/ 6.15.1 software package with the calibration database released on October 10, 2015. We analyzed so-called cleaned events that had passed the following standard event selection criteria: (a) Only ASCA-grade 0, 2, 3, 4, and 6 events were accumulated with hot and flickering pixels removed. (b) More than 436 s had elapsed since passing through the South Atlantic Anomaly. (c) The pointing directions were at least $5^{\circ}$ and $20^{\circ}$ above the rim of the Earth during the night- and daytime, respectively. To further reduce the non-X-ray background (NXB), we also required that (d) the geomagnetic cutoff rigidity exceeded 6 GV. Details concerning the observation and net exposures of the screened events are summarized in Table 1.

In order to complement the Suzaku observations which did not cover the entire PWN (see Figure 1), we also analyzed archival XMM-Newton data pointing at the position of PSR J2021+3651. XMM-Newton is equipped with two types of X-ray CCD, PN and MOS — both of which have sensitivities within 0.15–12 keV when combined with the X-ray telescope. Although the background is rather high and unstable owing to satellite’s highly elliptical orbit, XMM-Newton has a larger field of view (FOV) and effective area and is, therefore, complementary to Suzaku-XIS. Because PN was operated in timing mode to study the PSR J2021+3651, which is not suitable for studying PWNs, we used only the MOS data. Among two MOS CCD cameras, we focused on MOS2 data since one CCD chip of MOS1 that covers part of the PWN was not functional in this observation. The SAS 444https://www.cosmos.esa.int/web/xmm-newton/download-and-install-sas 15.0.0 and ESAS 555http://heasarc.gsfc.nasa.gov/docs/xmm/xmmhp_xmmesas.html 13 software packages were used in analyzing the data. Details of the procedure for reducing and estimating the particle-induced background is given in Section 3.2. A summary of the XMM-Newton observation and net exposure is shown in Table 1.

3 Data Analysis and Results

3.1 Suzaku Data

3.1.1 X-ray Images and PWN-West Morphology

We extracted X-ray images from XIS3 (FI CCD), which has better imaging quality than BI CCD (XIS1) owing to its lower instrumental background (Mitsuda et al., 2007; Tawa et al., 2008). Although XIS0 also has good imaging quality, we did not use it to construct an image in order to avoid artifacts resulting from its unusable area ($\sim$1/4 of the CCD chip). We defined the soft and hard bands as 0.7–2 and 2–10 keV, respectively, and excluded the corners of the CCD chips illuminated by the ${\rm{}^{55}Fe}$ calibration sources. We then estimated the NXB contribution from the nighttime Earth data and subtracted it from the images using xisnxbgen (Tawa et al., 2008). Vignetting was then corrected by dividing the soft- and hard-band images by flat sky images simulated at 1 and 4 keV, respectively, using the XRT+XIS simulator xissim (Ishisaki et al., 2007). In the flat image simulations, we assumed a uniform intensity of ${\rm 1~{}photon~{}s^{-1}~{}cm^{-2}~{}sr^{-1}}$, and therefore, the approximate unit of the obtained vignetting-corrected image is ${\rm photons~{}s^{-1}~{}cm^{-2}~{}sr^{-1}}$. The obtained images are shown in Figure 2, in which smoothing with a Gaussian kernel of $\sigma=0.^{\prime}28$ is applied for visualization.

Extended emission from the western part of the PWN is apparent in our first observation (S1, Figure 1), but no obvious extended emission is seen in our second observation (S2, Figure 1). We also identified two bright sources in S1: PSR J2021+3651, located at the east edge of the CCD chip, and a bright field star, USNO-B1.0 1268-044892 (already reported by Van Etten et al., 2008), which is seen mainly in the soft band. The PWN emission is roughly along the south side of the CCD that was tilted by $71\fdg 4$ east from the north, suggesting that the major axis of the X-ray PWN is almost parallel to that of VER J2019+368. To examine the source extent quantitatively, we defined $5^{\prime}\times 1^{\prime}$ rectangles 666In 15th to 17th regions (from the pulsar) in S1 observation, we used $4^{\prime}\times 1^{\prime}$ (instead of $5^{\prime}\times 1^{\prime}$) rectangles to avoid the corners illuminated by the ${\rm{}^{55}Fe}$ calibration sources. as shown in Figure 2 and calculated the PWN count rate profile. We used both XIS0 and XIS3 to conduct a morphology analysis of the PWN emission along its major axis. The size and position of the rectangles were chosen to avoid the unusable area of XIS0. We also removed two sources (PSR J2021+3651 and USNO-B1.0 1268-044892) in S1 and one hard source (presumably the background active galactic nucleus) in S2: the radius of the circles for exclusion was $120^{\prime\prime}$, $90^{\prime\prime}$, and $90^{\prime\prime}$ for PSR J2021+3651, USNO-B1.0 1268-044892, and the hard source seen in S2, respectively. As PSR J2021+3651 was located near the edge of the CCD chip, its position was not discernible from the XIS image. As the position accuracy of the XIS is known to be ${\sim}20^{\prime\prime}$ (Uchiyama et al., 2008), we did not use the reported position of the pulsar; instead, we referred to USNO-B1.0 1268-044892 and obtained the shifts as $-22^{\prime\prime}$ and $-9^{\prime\prime}$ in RA and DEC, respectively, and determined the position of PSR J2021+3651 in our image. From the position of PSR J2021+3651 toward the southwest (with the orientation of $108\fdg 6$ west from the north), we defined 17 rectangles in S1 and 11 rectangles in S2, with three rectangles overlapped. We then examined the morphology of the PWN up to $25^{\prime}$ from the position of the pulsar after first subtracting the NXB estimated by xisnxbgen (Tawa et al., 2008). The remaining X-ray background — presumably the cosmic X-ray background (CXB) and the Galactic ridge X-ray emission (GRXE) (e.g., Worrall et al., 1982; Warwick et al., 1985; Koyama et al., 1986), which are expected to be almost uniform within the XIS FOV — was estimated using $10^{\prime}\times 4^{\prime}$ rectangles in each observation as shown in Figure 2. The NXB-subtracted background count rate was subtracted from the NXB-subtracted source count rate with vignetting taken into account. Finally, the obtained count rate of each bin was corrected for vignetting and the region size (normalized to the count rate of the ninth rectangle from the pulsar which is the closest to the FOV center of the S1 observation), as summarized in Figure 3. It is seen from the figure that the PWN emission extends from the pulsar in the southwest direction by up to $15^{\prime}$ in the soft band and $18^{\prime}$ in the hard band. Although we corrected for the vignetting effect, it is severe at high energies and near the edge of the XIS (Serlemitsons et al., 2007). Therefore, we concluded rather conservatively that the PWN emission extends from the pulsar in the southwest direction at least up to $15^{\prime}$ in both the soft and hard bands.

We also examined the PWN morphology in the minor axis direction, and in the region between $3^{\prime}$ and $6^{\prime}$ from the pulsar along the major axis to avoid contamination from the pulsar, as shown in Figure 4(a). Because the XIS0 has an unusable area on the south side of the FOV, we only used XIS3 and the hard band in order to avoid emissions from USNO-B1.0 1268-044892. Background subtraction and vignetting correction were conducted in the same manner as in the morphology study along the major axis. The obtained count rate profile is shown in Figure 4(b), in which distance is measured from the south edge of the XIS toward the north, and the bins from $3^{\prime}$ to $8^{\prime}$ are located within rectangles used to study the morphology along the major axis. It is seen from the figure that the PWN emission has a source extent of at least $10^{\prime}$ along the minor axis; thus, these results show for the first time that the western region of the PWN has a source extent of at least $15^{\prime}$ and $10^{\prime}$ along the major and minor axis, respectively. The count rates in the $3^{\prime}$–$8^{\prime}$ bins (i.e., those within the area of study along the major axis) and within $0^{\prime}$–$10^{\prime}$ bins (the entire PWN emission) are $12.64\pm 0.74$ and $18.60\pm 1.02~{}{\rm c~{}s^{-1}}$, respectively, giving a ratio of $1.47\pm 0.12$ to convert the flux within a region of $5^{\prime}$ width to that of the entire PWN-West emission (see also Section 3.3).

3.1.2 Spectrum of the PWN-West

As described in Section  3.1.1, we confirmed that the PWN-west region extends up to (at least) $15^{\prime}$ westward from the pulsar. We then extracted spectra obtained by three CCD cameras (XIS0, XIS1, and XIS3) for 15 rectangles from S1, starting with the rectangle closest to the pulsar and with two point sources excluded, in order to maximize the photon statistics while avoiding the unusable area of XIS0. On the basis of the morphology of the major axis shown in Figure 3, we assumed a linear decrease in intensity (from 1 to 0 in relative) from $0^{\prime}$ to $15^{\prime}$ in calculating the ancillary response files (ARFs) using xissimarfgen (Ishisaki et al., 2007). The losses of effective area owing to the exclusion of point sources and area illuminated by calibration sources were taken into account in calculating the ARFs. In the spectral analysis, the response matrix files (RMFs) were calculated using xisrmfgen, and the integrated NXB spectrum over the source spectrum region was estimated using xisnxbgen (Tawa et al., 2008) and subtracted from the source spectrum. As the NXB-subtracted X-ray spectrum was expected to suffer from the (X-ray) background owing to the CXB and GRXE, the background was estimated again using the $10^{\prime}\times 4^{\prime}$ source-free region (see Figure 2). We first subtracted the NXB contribution from the background spectrum and then subtracted the NXB-subtracted background spectrum from the NXB-subtracted source spectrum with vignetting at 2 keV taken into account. The vignetting correction factors at 1 and 4 keV differ from that at 2 keV by only ${\sim}3\%$ and ${\sim}2\%$, respectively. The obtained spectrum was well fitted [reduced chi-square $\chi^{2}/{\rm degrees~{}of~{}freedom~{}(DOF)}=211.1/188$] by an absorbed power-law model (${\tt wabs}\times{\tt pow}$ in XSPEC), as shown in Figure 5 and Table 2. The best-fit hydrogen column density of the photoelectric absorption, $N({\rm H})=(8.2^{+1.3}_{-1.1})\times 10^{21}~{}{\rm cm^{-2}}$, 777Here and hereafter, errors are calculated for single-parameter 90% confidence limit. was consistent with the absorption toward the vicinity of the pulsar measured by the Chandra X-ray Observatory as reported by Van Etten et al. (2008) ($6.7^{+0.8}_{-0.7}\times 10^{21}~{}{\rm cm^{-2}}$). Therefore, we confirmed that the extended X-ray emission comes from the PWN around PSR J2021+3651. In order to examine the possible spectral change along the major axis, we divided the source region into five segments, each $3^{\prime}$ length. We repeated the same analysis procedure [response calculation, NXB subtraction, and X-ray background subtraction with vignetting correction] as described above and fit each of the spectra with an absorbed power-law model. We first let the absorption free to vary in each region and obtained the parameters as summarized in Table 2. Although there seems to be a slight softening of the spectra in outer regions (${\rm distance\geq 9^{\prime}}$), we observe a correlated increase of the absorption and the photon indices of all five subregions are consistent with that of the whole spectrum within statistical errors. We also fixed the absorption at $8.2\times 10^{21}~{}{\rm cm^{-2}}$, which was the best-fit value of the whole spectrum. As is seen from Table 2, the photon index does not change significantly over the entire western region of the PWN. We can also see that the intensity gradually decreases in a manner approximately proportional to the distance from the pulsar.

3.2 XMM-Newton Data

3.2.1 X-ray Images and PWN Morphology

In order to examine the overall properties of the PWN, we also analyzed archival data produced when the XMM-Newton was aimed at the position of PSR J2021+3651 (see Table 1). To reduce the particle-induced background and estimate the residual background as accurately as possible, we processed data using the ESAS software package. Details of the XMM-Newton CCD-camera background and analysis procedures for extended objects can be found in Kuntz & Snowden (2008) and Snowden & Kuntz (2014). We first excluded a period of time in which the data was severely contaminated by highly fluctuating background induced by soft protons, by making a light curve in 2.5–8.5 keV from the entire FOV, and created a count map. A cut was made by setting the threshold at $\pm 1.5\sigma$ from the average, resulting in the net exposure of 83.4 ks. We next created a count map due to the quiescent particle background (QPB) based on data obtained when the filter wheel was in the closed position (FWC data), subtracted it from the cleaned count map, and then divided the subtracted map by an exposure map to correct the vignetting. The procedures described above were made by running mos-filter and adapt commands. The obtained background-excluded/subtracted and exposure-corrected images in the soft-band (0.7–2 keV) and hard-band (2–10 keV) are shown in Figure 6, showing that the western and eastern regions of the PWN have a similar source extent. Three bright sources are also identifiable: PWN J2021+3651 in the middle of the image, WR 141 to the northeast of the pulsar, and USNO-B1.0 1268-044892 to the southwest of the pulsar.

To examine the morphology of the western and eastern regions of the PWN, we defined 12 rectangles of $5^{\prime}\times 1^{\prime}$ in each region as shown in the figure, similar to those used for the Suzaku-XIS data analysis. We excluded PSR J2021+3651 using a circular masking region with a radius of $60^{\prime\prime}$, and WR 141 and USNO-B1.0 1268-044892 using circular regions with a radius of $45^{\prime\prime}$. We also excluded several less-bright sources as using circular regions of $30^{\prime\prime}$ radius indicated by the figure. Although we applied the temporal filtering based on the light curve to reduce the soft proton contamination and subtracted the QPB, there remains non-negligible residual background induced by soft protons. To estimate the residual background, we extracted the spectrum from the entire FOV and fit it with a model to represent the X-ray emission plus residual soft-proton background modeled as a simple power-law convolved with the response matrix of diagonal unity elements (Kuntz & Snowden, 2008; Snowden & Kuntz, 2014). The details of the procedure and best-fit model parameters are given in Appendix A. In the following analysis of the morphology and spectrum of the PWN, the spectral index of the soft proton contamination is fixed to the best fit value for the entire FOV, and the normalization is scaled by using the proton-scale command.

The X-ray background (presumably the CXB and GRXE) for the PWN was estimated by calculating the count rate of the background region ($5^{\prime}\times 4^{\prime}$ rectangle located in northwest of the pulsar shown in Figure 6) after the QPB and the residual soft-proton background estimated from the entire FOV was subtracted (see above). We thus obtained the count rate profile of the western and eastern part of the PWN, with the QPB contribution estimated based on FWC data and subtracted, the residual soft-proton contamination estimated from the entire FOV and subtracted, and the X-ray background estimated from the background region and subtracted with the vignetting taken into account, as summarized in Figure 7. As was done for the morphology analysis by Suzaku data (Section 3.1.1), the count rate of each bin was corrected for vignetting and the region size (normalized to the entire rectangle of $5^{\prime}\times 1^{\prime}$ closest to the pulsar in the PWN-west). It is seen from the figure that the PWN emission extends from the pulsar in the southwest and northeast directions by up to at least $12^{\prime}$ in both the soft (0.7–2 keV) and hard (2–10 keV) bands. We can also see that the count rate profile of the PWN-west is roughly the same as that seen by the Suzaku morphology analysis, while the decrease of the count rate is less pronounced in the PWN-east. We also note that the estimated ratio of the X-ray count rate (PWN, CXB, and GRXE) to the NXB count rate (QPB and residual soft proton background) at the western edge of the PWN in XMM-Newton image ($11^{\prime}\mbox{--}12^{\prime}$ away from the pulsar) is about 0.29 in 2–10 keV, while that seen in Suzaku data is about 4.3. Therefore the XMM-Newton data might suffer from the larger systematic uncertainty of the NXB. This could be why we observe an enhancement of the intensity in western and eastern edges of the XMM-Newton image in the hard band (Figure 6b), the former of which was absent in Suzaku image (Figure 2b).

3.2.2 Spectrum of the PWN

We then proceed to the spectral analysis. We first extracted spectra for the western and eastern parts of the PWN for 12 rectangles starting with those closest to the pulsar with bright spots excluded (see Figure 6). On the basis of the morphology along the major axis shown in Figure 3 and 7, we assumed a linear decrease in intensity (from 1 to 0 in relative) from $0^{\prime}$ to $15^{\prime}$ in calculating the ARFs using the arfgen command for the PWN-west, and another linear decrease in intensity (from 1 to 0.5 in relative) from $0^{\prime}$ to $15^{\prime}$ for the PWN-east. The losses of the effective area owing to the exclusion of point sources were taken into account in calculating the ARFs. In the spectral analysis, the QPB contribution was estimated using mos-filter and subtracted from the source spectrum, the residual soft-proton contamination was estimated using the data of the entire FOV and subtracted with the scale factor calculated by proton-scale, and the contribution from the X-ray background were estimated by simultaneously fitting the spectra of the source and background regions shown in Figure 6. The extended PWN emission was modeled by an absorbed power-law model (${\tt wabs}\times{\tt pow}$ in XSPEC) in the source spectrum. We also analyzed the spectrum of the so-called ”Arc” (Van Etten et al., 2008) by the same procedure with a flat intensity profile assumed in calculating the ARF. The obtained spectra are shown in Figure 8 and parameters of the source spectra are summarized in Table 3. Detailed descriptions of the spectral modeling and obtained parameters of the background region are given in Appendix A. It is seen from Table 3 that the absorption ($N({\rm H})$), photon index ($\Gamma$), and flux are similar between the PWN-west and PWN-east, and similar $N({\rm H})$ and $\Gamma$ are obtained for the Arc, supporting the same physical origin of three regions. $N({\rm H})$ and $\Gamma$ of the PWN-west are similar to those measured by Suzaku (Table 2), and the obtained flux in 2–10 keV agrees with that integrated over $0^{\prime}\mbox{--}12^{\prime}$ by Suzaku in ${\sim}10\%$. We also divided the western/eastern regions into four segments (each $3^{\prime}$ length) and fitted each of the spectra as we did for Suzaku data and summarize the results in Table 4, in which the significant change of the spectral index was not seen.

3.3 Summary of X-ray Data Analysis Results

Before proceeding to the discussion (Section 4), let us summarize the results of the X-ray data analysis.

Even with Suzaku-XIS, no extended emission was found in the western region of TeV emission (Section 3.1.1). 2. 2.

The source extent of the PWN-west was measured to be $15^{\prime}\times 10^{\prime}$ by Suzaku-XIS, with a linear decrease of the intensity from $0^{\prime}$ to $15^{\prime}$ (Section 3.1.1). The XMM-Newton data indicate that the PWN-east has a flatter intensity profile up to $12^{\prime}$, beyond which the source extent is not constrained. The Arc has a source extent of ${\sim}9^{\prime}\times 2.^{\prime}5$. (Section 3.2.1) 3. 3.

The orientation of the PWN major axis is ${\sim}71\fdg 4$ east from the north (Section 3.1.1). 4. 4.

The PWN-west spectrum is represented by an absorbed power-law with $N({\rm H})=(8.2^{+1.3}_{-1.1})\times 10^{21}~{}{\rm cm^{-2}}$ and $\Gamma=2.05\pm 0.12$. No significant change of $\Gamma$ was found inside the region (Section 3.1.2). With the results of XMM-Newton spectral analysis (Section 3.2), we confirm that the PWN-east and the Arc have similar spectral parameters to those of the PWN-west. (Section 3.2.2) 5. 5.

The 2–10 keV observed flux $f({\rm 2\mbox{--}10~{}keV})$ in the region of $15^{\prime}\times 5^{\prime}$ to the west of the PWN was measured to be $2.6\times 10^{-12}~{}{\rm erg~{}s^{-1}~{}cm^{-2}}$ (Section 3.1.2), giving the absorption-corrected flux $F({\rm 2\mbox{--}10~{}keV})$ of $2.8\times 10^{-12}~{}{\rm erg~{}s^{-1}~{}cm^{-2}}$. On the basis of the morphology along the minor axis (Section 3.1.1), we obtained $F({\rm 2\mbox{--}10~{}keV})\sim 4.1\times 10^{-12}~{}{\rm erg~{}s^{-1}~{}cm^{-2}}$ for the entire PWN-west. Although it was not possible to constrain the extent of the PWN-east beyond the XMM-Newton FOV, the results of XMM-Newton spectral analysis (Section 3.2.2) suggest that the PWN-east plus Arc has $F({\rm 2\mbox{--}10~{}keV})\sim 3.0\times 10^{-12}~{}{\rm erg~{}s^{-1}~{}cm^{-2}}$ at least, giving the lower limit of $F({\rm 2\mbox{--}10~{}keV})$ for the overall PWN (the PWN-west, PWN-east, and Arc) to be ${\sim}7\times 10^{-12}~{}{\rm erg~{}s^{-1}~{}cm^{-2}}$.

4 Discussion

4.1 Properties of the X-ray PWN

Here we describe the properties of the X-ray PWN, and the implications of these properties, without discussing its relation with VER J2019+368.

The photon index of the PWN-west we obtained, $\Gamma=2.05\pm 0.12$, is significantly larger than that of the PWN emission close to the pulsar measured by Chandra as reported by Van Etten et al. (2008); they obtained $\Gamma=1.0\mbox{--}1.5$ within ${\sim}10^{\prime\prime}$ from the pulsar (“Inner nebula”), and $\Gamma\sim 1.7$ for their “Jet” and “Outer nebula-east” (which is within ${\sim}1.5^{\prime}$ of the pulsar). The photon index we measured is consistent, however, with $\Gamma=1.93\pm 0.13$ for their “Outer nebula-west” (which is 1.5–3*′* from the pulsar). As no significant change of $\Gamma$ is observed up to $15^{\prime}$ from the pulsar toward the southwest direction in Suzaku data, we can conclude that the CR electrons accelerated at the PWN termination shock (${\sim}10^{\prime\prime}$ from the pulsar; Van Etten et al., 2008) suffer from synchrotron cooling close to the pulsar (within ${\sim}1.5^{\prime}$) but propagate outward without significant cooling. The XMM-Newton data indicate similar conclusions on the PWN-east; the photon index ($\Gamma=2.03\pm 0.10$) is larger than that of the pulsar, Jet, and Outer nebula-east and no significant change of $\Gamma$ is observed up to $12^{\prime}$ toward the northeast.

Through observations by Suzaku (Section 3.1), XMM-Newton (Section 3.2), and Chandra (Van Etten et al., 2008), the absorption of the X-ray PWN was found to be $(6\mbox{--}9)\times 10^{21}~{}{\rm cm^{-2}}$, significantly lower than the total Galactic absorption in the direction of Cygnus-X of ${\geq}2.5\times 10^{22}~{}{\rm cm^{-2}}$ estimated by Mizuno et al. (2015) using X-ray source spectra and $\gamma$-ray data. Therefore, the pulsar and its PWN are unlikely to be located at a distance ${\geq}10~{}{\rm kpc}$, as was inferred from radio data (see Section 1). Instead, we adopt the distance $d=1.8^{+1.7}_{-1.4}~{}{\rm kpc}$ estimated by Kirichenko et al. (2015) based on the absorption-distance relation using red-clump stars in the direction of the pulsar.

4.2 Relation to VER J2019+368

We first discuss implications of the determined properties in the X-ray and TeV $\gamma$-ray regimes. We then examine particle transport (and magnetic fields), primarily within the X-ray PWN, and implications. We finally present a possible model to explain the multiwavelength data.

First, the fact that the major axes of the X-ray PWN and VER J2019+368 are almost parallel (Section 3.1.1) strongly support that the X-ray PWN is physically associated with VER J2019+368. If the X-ray PWN is a counterpart of the TeV emission, TeV $\gamma$ rays are likely to be produced by the inverse Compton (IC) scattering by X-ray producing CR electrons. In the case of the PWN synchrotron/IC scenario, temporarily neglecting the details of the electron spectrum and the Klein-Nishina effect produces a ratio of X-ray to TeV $\gamma$-ray luminosities given by the ratio of magnetic field energy density to photon field energy density, $U_{\rm mag}/U_{\rm ph}$. The absorption-corrected PWN flux in the X-ray regime, $F({\rm 2\mbox{--}10~{}keV})\sim 7\times 10^{-12}~{}{\rm erg~{}s^{-1}~{}cm^{-2}}$, is close to the TeV $\gamma$-ray flux of VER J2019+368 [$F({\rm 1\mbox{--}10~{}TeV})\sim 6.7\times 10^{-12}~{}{\rm erg~{}s^{-1}~{}cm^{-2}}$], indicating either that $U_{\rm mag}$ is close to the energy density of the cosmic microwave background (CMB) or that the average magnetic field $B$ of the PWN is rather low and close to the typical interstellar magnetic field, $B\sim 3~{}\mu{\rm G}$. This value should be taken as a lower limit since the X-ray observations did not cover the whole TeV-emitting region (see Figure 1). Nevertheless, a much larger value of the magnetic field is unlikely, since the two Suzaku observations covered the central region of the TeV emission. Hereafter we express physical quantities with $B$ normalized by $3~{}\mu{\rm G}$.

From the spectra and morphologies in the X-ray and TeV $\gamma$-ray regimes, we can constrain the properties of the accelerated CR electrons. Hereafter, we assume a constant injection of accelerated CR electrons into uniform magnetic and radiation fields over the lifetime of the pulsar for simplicity.

From discussions in, e.g., Longair (2011) and de Jager & Djannati-Ataï (2008), the characteristic energy of X-rays owing to synchrotron radiation ($\epsilon_{\rm X}$) in the magnetic field $B$ is related to the electron energy, $E_{e}$, as

[TABLE]

On the other hand, the typical energy of $\gamma$ rays ($\epsilon_{\gamma}$) generated by IC scattering of the CMB photons is related to $E_{e}$ as,

[TABLE]

Equations (1) and (2) indicate that electrons with $E_{e}\geq 100~{}{\rm TeV}$ are required to generate synchrotron X-rays above 1 keV with $B\sim 3~{}\mu{\rm G}$, while electrons with $E_{e}\leq 50~{}{\rm TeV}$ produce $\gamma$ rays below 10 TeV. As we obtained $\Gamma=2.05$ for the X-ray PWN and $\Gamma=1.75$ was reported for TeV $\gamma$-ray emission, there must be a spectral break of accelerated CR electrons at around 50–100 TeV. Although the TeV $\gamma$-ray photon index has a rather large uncertainty of ${\sim}0.2$ (Aliu et al., 2014), the Klein-Nishina effect softens the electromagnetic spectrum more than in the Thomson regime, in which the photon index of IC emission is the same as that of the synchrotron radiation with the same CR electron spectral index. Therefore, our hypothesis of a spectral break is robust. The process causing this spectral break is likely to be a synchrotron cooling, in which the break energy $E_{\rm bk}$ is related to the injection time $t_{0}$ and magnetic field $B$ as

[TABLE]

where $t_{0}$ is normalized to the characteristic age of the pulsar. The electron spectral index changes by 1, and the synchrotron radiation and IC emission each change by 0.5, indicating that the difference in spectral slopes between X-rays and TeV $\gamma$ rays can be naturally explained by the canonical age of the pulsar (characteristic age of 17.2 kyr) and $B=3~{}\mu{\rm G}$. We should also take into account the cooling of CR electrons during propagation; for electrons producing X-rays above 1 keV, the main mechanism for this is synchrotron cooling. Then, using Equation (1), the cooling time $\tau$ can be expressed as

[TABLE]

This indicates that the lifetimes of CR electrons producing X-rays at 1 and 10 keV (under $B=3~{}\mu{\rm G}$) are 10.5 and 3.3 kyr, respectively. If we also take into account the cooling owing to IC scattering of the CMB and infrared background (based on the blackbody radiation at a temperature of 30 K and energy density of $0.3~{}{\rm eV~{}cm^{-3}}$; see below) using the procedure described in Moderski et al. (2005), the true lifetimes of electrons producing 1 and 10 keV X-rays are found to be 7.9 and 3.0 kyr, respectively. Therefore, Equation (4) is valid to within 25%.

Let us then discuss particle transport and its implications. The CR electrons are transported via either diffusion or advection caused by the pulsar wind. If advection is the dominant process, high-energy electrons with shorter lifetimes [Equation (4)] will make it closer from the pulsar. This implies a spectral softening not seen in our detailed study of the X-ray spectrum (Section 3.1.2). Therefore, the highest energy electrons we consider propagate over a distance ${\geq}15^{\prime}$ during their lifetime. Since the angular extent of $15^{\prime}$ corresponds to $8~{}{\rm pc}\left(\frac{d}{1.8~{}{\rm kpc}}\right)$, and the lifetime of electrons producing 10 keV X-rays due to synchrotron radiation is $3.3~{}{\rm kyr}\left(\frac{B}{3~{}\mu{\rm G}}\right)^{-1.5}$, the advection velocity divided by the speed of light ($\beta_{\rm adv}$) should satisfy

[TABLE]

In this scenario, the absence of X-ray emission beyond the peak position of the TeV emission is due to the lower surface brightness of synchrotron X-rays caused by the lower magnetic field or lower CR electron density. The scenario can naturally explain the larger size of TeV emission produced by electrons of lower energy (longer lifetime). In the case of diffusion-dominated scenario, the electrons propagate the diffusion length of $\sqrt{2D\tau}$, where $D$ and $\tau$ are the diffusion coefficient and the electron lifetime, respectively. Then we can constrain $D$ as we did to constrain $\beta_{\rm adv}$. Let us first examine the case of energy-independent diffusion, as predicted by, e.g., Porth et al. (2016) through three-dimensional magnetohydrodynamics simulations. In order for the diffusion length to exceed the length of the X-ray PWN, even for electrons producing 10 keV X-rays, we obtain

[TABLE]

Like the advection-dominated scenario, the absence of X-ray emission beyond the TeV emission peak is due to the lower magnetic field or lower CR electron density, and the larger size of TeV emission is due to the cooling of electrons producing X-rays. Alternatively, diffusion can naturally explain the apparent lack of spectral softening, if $D$ depends on the particle energy as $D\propto E_{\rm e}^{\delta}$ with $\delta\sim 1$ (e.g., Van Etten et al., 2011). If $\delta=1$ the diffusion coefficient $D$ can be expressed as $D=\frac{1}{3}\lambda_{g}c\eta$, where $\lambda_{g}$ is the electron gyroradius, $c$ is the speed of light, and the parameter $\eta$ is related to the degree of magnetic turbulence. By substituting the physical constants and also using Equation (1), we obtain

[TABLE]

Then, in order for the source extent not to exceed the diffusion length in electron lifetime ($\sqrt{2D\tau}$), we obtain [by substituting Equation (4)]

[TABLE]

Therefore, under the condition of $B\sim 3~{}\mu{\rm G}$ and $d\sim 1.8~{}{\rm kpc}$, $\eta\sim 1$ (i.e., close to Bohm limit) is required, suggesting that the magnetic field is highly turbulent. Energy-dependent diffusion alone, however, is not able to explain the larger size of the TeV emission. We thus constrain the advection velocity $\beta_{\rm adv}$ or the diffusion coefficient $D$ from the morphology of the X-ray PWN.

On the basis of the discussions above (in particular, regarding $E_{\rm bk}$ and $B$), under the assumption of the constant injection and uniform magnetic field, we present a possible multiwavelength spectral model in Figure 9 in which the electron spectrum is assumed to be a power-law with a photon index of 2.1 below 80 TeV and 3.1 above 80 TeV and with an exponential cutoff at 1 PeV. Contributions from synchrotron radiation and IC scattering are computed based on Crusius & Schlickeiser (1986) and Blumenthal & Gould (1970) respectively. We adopted $B=3~{}\mu{\rm G}$ and adjusted the model normalization to explain the entire X-ray PWN flux in the region 2–10 keV (${\sim}7\times 10^{-12}~{}{\rm erg~{}s^{-1}~{}cm^{-2}}$; see Section 3.3). To calculate the IC scattering emission, we referred to the radiation field model of Porter et al. (2008). Because the source distance from the Galactic center is estimated to be 8.2 kpc in our case where $l=75\arcdeg$ and $d=1.8~{}{\rm kpc}$, we adopted their model at the solar circle and assumed a CMB and infrared background with temperature 30 K and energy density $0.3~{}{\rm eV~{}cm^{-3}}$. We also overlaid the $\gamma$-ray emission from the pulsar and an upper limit of the PWN in the GeV band taken from Abdo et al. (2009). Paredes et al. (2009) reported extended radio emission of ${\sim}700~{}{\rm mJy}$ at 1.4 GHz in the vicinity of VER J2019+368. Although they did not provide information on the position and spatial extent, we plot their flux for reference. As is seen from Figure 9, the model explains about 80% of the TeV emission, indicating that the X-ray PWN is a major contributor to VER J2019+368. Considering the assumptions we made (constant injection of CR electrons into uniform magnetic and radiation fields over the pulsar lifetime), limited coverage of X-ray observations (see Figure 1), and the apparent offset of the pulsar from the peak of the TeV emission (which cannot be explained by our simplified scenario), we do not rule out X-ray emission from the nebula further out and/or the confusion of TeV source(s) physically unrelated to the X-ray PWN. Further observations in X-rays and TeV $\gamma$ ray are worthwhile to fully understand the system. In particular, TeV $\gamma$-ray observations at the better sensitivity and angular resolution by the Cherenkov Telescope Array (CTA) (Actis et al., 2011) are anticipated to reveal the TeV $\gamma$-ray properties in more detail.

In the discussions above, we have assumed constant injection of electrons into uniform magnetic fields for simplicity. If the pulsar has already experienced significent energy losses, it injected more electrons in the past, therefore the ratio of TeV to X-ray flux is increased. In order for the predicted TeV flux not to exceed the observed value, a magnetic field larger than $B=3~{}\mu{\rm G}$ is required. In this case, $E_{\rm bk}\sim 80~{}{\rm TeV}$ (which is required to explain the different spectral indices between X-ray and $\gamma$ ray) can be achieved if the true age of the pulsar is younger (see also Section 4.3).

So far we have assumed only the PWN since no evidence of a host SNR is found (Van Etten et al., 2008). If the parent SNR is found in future, the discussion on particle transport and the relation between X-ray and TeV $\gamma$ ray might be affected.

4.3 Comparison with Other PWNe

Finally, we compare the properties of the X-ray PWN and VER J2019+368 with other X-ray PWNe associated with TeV $\gamma$ rays. According to Mattana et al. (2009) who compiled the properties of 14 PWNe, the $\gamma$ ray to X-ray energy flux ratio is approximately proportional to the pulsar characteristic age, owing to the effect of severe cooling on X-ray production. The energy flux ratio at 1–30 TeV and 2–10 keV in our case is $(11.6\times 10^{-12})/(7\times 10^{-12})\sim 1.6$, which is roughly consistent with their Figure 1. Therefore, the smaller flux and size of the X-ray region can be understood, as with other TeV-emitting PWNs, to be caused by faster cooling of X-ray-producing electrons. Bamba et al. (2010) studied eight PWNs with various characteristic ages, which are associated with TeV $\gamma$-ray sources; in particular, they studied the size of an X-ray PWN as a function of the pulsar characteristic age. They found a rather constant size up to $10~{}{\rm kyr}$ and then a gradual increase in size thereafter, possibly caused by an increase in advection speed or a decrease in the magnetic field turbulence as the pulsar/PWN grows older. Our measured size of the X-ray PWN, $8~{}{\rm pc}\left(\frac{d}{1.8~{}{\rm kpc}}\right)$, is consistent with those of their samples of similar characteristic age.

We should also compare with the archetype evolved PWN HESS J1825-137 and its extended X-ray PWN around PSR J1826-1334. The properties of the pulsar are similar to those of PSR J2021+3651. The pulse period and its derivatives are $101~{}{\rm ms}$ and $7.5\times 10^{-14}$, respectively (Clifton et al., 1992), giving the surface magnetic field of $2.8\times 10^{12}~{}{\rm G}$, the characteristic age of 21.4 kyr, and spin-down luminosity of $2.8\times 10^{36}~{}{\rm erg~{}s^{-1}}$. The PWN has also similar properties. The TeV PWN extends more than $1\arcdeg$ with the $\gamma$-ray peak position being offset from the pulsar and X-ray peak position by ${\sim}10^{\prime}$ (Aharonian et al., 2006). The source extent of the X-ray PWN was measured by Suzaku (Uchiyama et al., 2009; Van Etten et al., 2011) to be ${\sim}15^{\prime}$ towards the south (the northern part of the pulsar has not been observed; see Figure 1 of Van Etten et al., 2011). While the compact core of the X-ray PWN has a hard photon index of ${\sim}1.6$ (Gaensler et al., 2003), the outer part of the X-ray PWN does not exhibit significant spectral softening with the photon index of ${\sim}2$ (Uchiyama et al., 2009). There are, however, two distinct properties between two systems. Firstly, the energy fluxes in X-ray (2–10 keV) and TeV $\gamma$ ray (1–10 TeV) of HESS J1825-137 are $4.5\times 10^{-12}$ and $51\times 10^{-12}~{}{\rm erg~{}s^{-1}~{}cm^{-2}}$, respectively (Uchiyama et al., 2009; Aharonian et al., 2006), giving ${\sim}10$ times larger $F(1\mbox{--}10~{}{\rm TeV})/F(2\mbox{--}10~{}{\rm keV})$ ratio than that we obtained for VER J2019+368. Secondly, the photon index in TeV $\gamma$ ray of HESS J1825-137 is ${\sim}2.4$ on average and is significantly softer than that of VER J2019+368. These two facts can be explained naturally by assuming more severe cooling of electrons in HESS J1825-137. If the CR electrons of energies less than ${\sim}50~{}{\rm TeV}$ have already suffered from cooling, a softer spectral index in TeV and a larger ratio of the TeV $\gamma$-ray flux to the X-ray flux are expected. Probably VER J2019+368 has weaker mean magnetic field and/or the true age of the pulsar is younger. We also note that Van Etten et al. (2011) carried out a modeling of X-ray and TeV $\gamma$-ray data of HESS J1825-127 and constrained the electron injection history, profile of the magnetic field, advection velocity, and diffusion coefficient. Although such an extensive modeling is beyond the scope of our study, detailed study of the morphology in TeV $\gamma$ ray by future observations by CTA is anticipated to better understand the VER J2019+368 system.

5 Summary

We conducted deep X-ray observations of the VER J2019+368 region using Suzaku-XIS to examine the properties of the X-ray PWN around PSR J2021+3651 and to search for previously unknown extended X-ray emissions. We also analyzed archival XMM-Newton data to complement the Suzaku observations, which did not cover the entire region of VER J2019+368. We found that the total size of the X-ray PWN along the major axis is more than $27^{\prime}$: the PWN-west have a source extent of approximately $15^{\prime}\times 10^{\prime}$ with an orientation of its major axis nearly parallel to that of TeV emission, and the PWN-east extends up to at least $12^{\prime}$ from the pulsar. The PWN spectra were well fitted by an absorbed power-law for absorption at ${\sim}8\times 10^{21}~{}{\rm cm^{-2}}$ and a photon index of ${\sim}2$, with no obvious change in the index occurring within the X-ray PWN. The measured X-ray absorption favors the distance to the source to be much smaller than $10~{}{\rm kpc}$ inferred from radio data. Aside from the PWN around PSR J2021+3651, no extended emission was found by even Suzaku-XIS. The uniformity of the X-ray photon index constrains the advection velocity or the diffusion coefficient depending on the primary process of particle transport for X-ray-producing CR electrons. From the measured X-ray spectrum, reported TeV $\gamma$-ray spectrum and X-ray source extent, under the assumption of the constant injection of CR electrons into the uniform magnetic and radiation fields over the characteristic age of the pulsar, we obtained a rather low magnetic field of ${\sim}3~{}\mu{\rm G}$. Our synchrotron/IC model is able to explain ${\sim}80\%$ of the TeV emission, indicating that the X-ray PWN is a major contributor to VER J2019+368. To fully understand the nature of the extended TeV emission, higher-sensitivity and higher-resolution observations by facilities such as CTA will be useful.

We would like to thank S. Kisaka for valuable comments, and H. Katagiri for helping the calculation of multiwavelength spectral model. We also thank the referee for his/her valuable comments, and the Suzaku team members and XMM-Newton team members for their dedicated support of the satellite operation and calibration. This work was partially supported by JSPS Grant-in-Aid for Scientific Research Grant Numbers JP25287059 (T.M.), JP15K05088 (R.Y.), and JP26800160 (K. H.).

Appendix A Detailed Descriptions and Parameters of the XMM-Newton Data Analysis

In order to estimate the residual soft-proton background (Section 3.2), we accumulated the spectrum from the entire FOV and modeled it with a model which consists of a simple power-law (pow in XSPEC to model the residual soft-proton contamination), an absorbed power-law (${\tt wabs}\times{\tt pow}$ to model the CXB), thin-thermal plasma emission (${\tt wabs}\times{\tt apec}$ to model the hard-temperature emission of the GRXE), two absorbed power-laws (${\tt wabs}\times{\tt pow}$ to reproduce the spectra of inner nebula and outer nebula reported by Van Etten et al. (2008)), and another absorbed power-law (${\tt wabs}\times{\tt pow}$ to approximate the sum of point sources, emission from the pulsar, and rest of the PWN emission). The response matrix of diagonal unity elements is assumed in the first component. For simplicity, a flat intensity profile is assumed for the others. Some parameters were fixed to typical values (parameter values of the CXB were taken from Kuntz & Snowden (2008), those of the GRXE were referred to Mizuno et al. (2015), and those of the inner/outer nebulae were taken from Van Etten et al. (2008)). Since we aim to constrain the residual soft-proton background which is prominent in high energy, we focused on data in 3–12 keV. The obtained best-fit parameters are summarized in Table A1. In the rest of the XMM-Newton data analysis, the spectral index of the residual soft-proton background is fixed to what obtained here ($\Gamma=0.21$) with the normalization scaled using the proton-scale command.

We then estimated the X-ray background to examine the count rate profile by calculating the count rate of the background region ($5^{\prime}\times 4^{\prime}$ rectangle located in the northwest of the pulsar shown in Figure 6) after the QPB and the residual soft-proton background (estimated from the entire FOV as above) was subtracted. The background count rate was then subtracted from the count rate in the source region with the vignetting taken into account. The obtained count rate profiles of the PWN are summarized in Figure 7, in which data in 1.4–1.6 and 1.7–1.8 keV were discarded to reduce the contamination from the instrumental background due to Al ${\rm K_{\alpha}}$ and Si ${\rm K_{\alpha}}$ fluorescent lines (Kuntz & Snowden, 2008; Snowden & Kuntz, 2014).

We finally analyzed the spectrum of the PWN-west, PWN-east, and Arc. In addition to the CXB (${\tt wabs}\times{\tt pow}$) and high-temperature emission of the GRXE (${\tt wabs}\times{\tt apec}$), we added two more thin-thermal plasma models (${\tt wabs}\times{\tt apec}$ to model the soft-temperature emission of the GRXE and local diffuse X-ray emission), gaussian lines (gauss) at 1.49 and 1.75 keV to model the fluorescent background of Al and Si, and a line (gauss) at 0.65 keV to model the contribution of the Solar-Wind Charge eXchange (SWCX, Kuntz & Snowden, 2008; Snowden & Kuntz, 2014). Those X-ray backgrounds were estimated by carrying out the joint spectral-fitting over the background region and the source region, with the parameters coupled with the vignetting taken into account. Again some parameters were fixed to typical values. The obtained best-fit parameters of the background region are summarized in Table A2, and those of the sources are given in Table 3.