SEARCHING FOR OVERIONIZED PLASMA IN THE GAMMA-RAY-EMITTING SUPERNOVA REMNANT G349.7+0.2

G349.7+0.2 was observed with Suzaku on 2011 September 19 with XIS for a total exposure time of 160 ks (observation ID: 506064010). The XIS consists of three active X-ray CCD detectors (XIS 0, 1, and 3). XIS0 and XIS3 have front-illuminated (FI) CCDs, while XIS1 contains a back-illuminated (BI) CCD.

For data reduction, we used the High Energy Astrophysics Software (HEASoft) package version 6.16 and the calibration database (CALDB) version 20140520. We used xspec version 12.8.2 (Arnaud 1996) for spectral analysis. The redistribution matrix and the auxiliary response functions were generated by xisrmfgen and xissimarfgen (Ishisaki et al. 2007), respectively.

Since G349.7+0.2 is located on the Galactic plane, the background of this SNR consists of the instrumental non-X-ray background (NXB), the Cosmic X-ray background, and the Galactic ridge X-ray emission (GRXE). Considering the position dependence of the GRXE, we constructed the background spectrum from the nearby blank sky (observation ID: 100028030) because the background is dominated by the GRXE (Uchiyama et al. 2013). The NXB spectra for the source and background were estimated from the Suzaku database of dark-Earth observations using the procedure of Tawa et al. (2008). The source and background spectra were made by subtracting the NXB spectra using mathpha. The NXB-subtracted background spectrum was subsequently subtracted from the source spectrum using xspec.

Figure 1 shows the Suzaku combined XIS0, XIS1, and XIS3 image of G349.7+0.2 in the 0.3–10.0 keV energy band. In order to study the spectral properties, we extracted the data from a $1\buildrel{\,\prime}\over{.} 7\times 1\buildrel{\,\prime}\over{.} 3$ size box centered at ${\rm R}{\rm .A}{\rm .}(2000)\;={{17}^{{\rm h}}}{{18}^{{\rm m}}}03\buildrel{\rm{s}}\over{.}7$ and ${\rm decl}.(2000)=-37{}^\circ 26^{\prime} 56\buildrel{\prime\prime}\over{.} 2$ for the E region and from a same-size box centered at ${\rm R}{\rm .A}{\rm .}(2000)\;={{17}^{{\rm h}}}{{17}^{{\rm m}}}57\buildrel{\prime\prime}\over{.} 8$ and ${\rm decl}.(2000)\;=$ $-37{}^\circ 25^{\prime} 37\buildrel{\prime\prime}\over{.} 9$ for the W region. These regions are shown with black rectangles in Figure 1. The E region corresponds to regions 1, 2, and 3 in Figure 1 of Lazendic et al. (2005) and the W region includes the regions 4, 5, and 6 in the same figure.

Zoom In Zoom Out Reset image size

For the E region, we first tried to fit the spectrum to a single-temperature plane-parallel shock model with variable abundances for all elements with $Z\leqslant 30$(VVPSHOCK in xspec) modified by the TBABS absorption model (Wilms et al. 2000) in xspec. The abundances are set to wilm (Wilms et al. 2000). In this fitting, the absorbing column density (${{N}_{{\rm H}}}$), upper limit of the ionization timescale (${{\tau }_{{\rm u}}}$), electron temperature ($k{{T}_{{\rm e}}}$), normalization, and abundances of Si, Ca, and Fe were left as free parameters. The other elements were fixed at their solar values (Wilms et al. 2000). Although this fit gave an acceptable ${{\chi }^{2}}$/degrees of freedom (dof) value of 448/429, residuals at ∼3.3 keV and ∼4.1–4.5 keV remained. The first residual around 3.3 keV corresponds to the Ar Lyα line. The other residual could be caused by the RRC corresponding to the H-like Ar.

Thus, we added a Gaussian and RRC (REDGE in xspec) components to the model. The line energy of the Ar Lyα was fixed at 3.3 keV, while the normalization was a free parameter. We note that the Gaussian line width parameter was fixed to zero eV. For the RRC component, the energy was fixed at 4.4 keV. After adding the Gaussian and RRC components to the model, the reduced ${{\chi }^{2}}$ value improved to 1.01 for 427 dof, which yields an F-test probability of 4.2 × 10⁻⁴.

We also applied the RP model VVRNEI in xspec for the E region. VVRNEI is a non-equilibrium recombining collisional plasma model with variable abundances for all elements ($Z\leqslant 30$). It is characterized by a constant electron temperature ($k{{T}_{{\rm e}}}$), initial temperature ($k{{T}_{{\rm init}}}$), and a single ionization parameter. N$_{{\rm H}}$, $k{{T}_{{\rm e}}}$, the abundances of Si, Ca, Fe, τ, and the normalization were free parameters during the fitting. We tried fitting by fixing $k{{T}_{{\rm init}}}$ to a range of values (2.0, 5.0, 10.0, and 50.0 keV) and for each fit we obtained similar ${{\chi }^{2}}$/dof and τ values for each of the $k{{T}_{{\rm init}}}$ values. This model gave a ${{\chi }^{2}}$/dof of 520.6/431, which required a large ionization timescale ($\tau \;\sim$ 2.8 × 10¹³ cm⁻³ s) showing that G349.7+0.2 is in ionization equilibrium and not in the RP phase.

As a next step, we applied an absorbed VVPSHOCK model for the W region. During the fitting, the absorption, electron temperature, ionization timescale, abundances of Si, Ar, Ca, and Fe, and normalization parameters were allowed to vary freely, while the other elements were fixed at their solar values (Wilms et al. 2000). We obtained an acceptable fit with a ${{\chi }^{2}}{\rm /dof}$ value of 539.7/529 and found no significant residuals. Thus, the spectral analysis of the E and W regions showed that there is no overionized plasma in G349.7+0.2. The spectral fits for the E and W regions are given in Figure 2 and the best-fit parameters are summarized in Table 1, where all of the errors are given at the 90% confidence level.

Zoom In Zoom Out Reset image size

Table 1.  Best-fit Parameters with Corresponding Errors at the 90% Confidence Level for the E and W Spectra of G349.7+0.2

Component	Parameters	E	W
TBABS	N$_{{\rm H}}$ (1022 cm−2)	8.7$_{-0.5}^{+0.3}$	9.5 ± 0.3
VVPSHOCK	kTe (keV)	1.2 ± 0.1	1.3 ± 0.1
	Si (solar)	1.4$_{-0.2}^{+0.3}$	1.2 ± 0.1
	Ar (solar)	1 (fixed)	1.2 ± 0.2
	Ca (solar)	2.0$_{-0.5}^{+0.7}$	1.3$_{-0.3}^{+0.4}$
	Fe (solar)	3.3$_{-1.5}^{+2.1}$	1.5 ± 0.4
	${{\tau }_{{\rm u}}}$ (1011 cm−3 s)	13.9$_{-6.3}^{+1.7}$	4.8$_{-0.9}^{+1.2}$
	Norm (10−2 photon cm−2 s−1)	5.2$_{-0.7}^{+0.4}$	5.8 ± 0.3
Ar Lyα	E (keV)	3.3 (fixed)	⋯
	σ (keV)	0.0 (fixed)	⋯
	Norm (10−5 photon cm−2 s−1)	1.4$_{-0.8}^{+0.7}$	⋯
RRC H-like Ar	E (keV)	4.4 (fixed)	⋯
	Norm (10−5 photon cm−2 s−1)	3.0$_{-1.5}^{+0.9}$	⋯
	${{\chi }^{2}}$/dof	432/427	539.7/529
	reduced ${{\chi }^{2}}$	1.01	1.02

Note. The Spectral Fits were Carried out in the 1.0–8.0 keV Energy Band.

Download table as:  ASCIITypeset image

We searched for the best-fit location of G349.7+0.2 within the ROI, which was found as R.A.(J2000) = 259509 ± 0013 and decl.(J2000) = −37455 ± 0013. The TS value of the best-fit position⁷ was enhanced by 2.7σ over the position of 2FGL J1718.1–3725 in the second Fermi-LAT catalog. Then, the model was refitted using the best-fit position to compute the TS map (Section 2.2.2 and Figure 3) and the spectrum.

To check the functional form of the spectral energy distribution (SED), we considered G349.7+0.2 as a point-like source. First, the PL function was fit to the data between 200 MeV and 300 GeV. Because the spectrum deviates from a PL function, we tested if the gamma-ray emission is better described by a log-parabola (LP) or a broken power-law (BPL) function, where their functional forms are as follows.

• 1.  
  Log-parabola:

  $$\begin{eqnarray*}&&F{{(E)}^{{\rm LP}}}={{N}_{{}^\circ }}{{(E/{{E}_{b}})}^{({{{\Gamma }}_{1}}+{{{\Gamma }}_{2}}{\rm ln} (E/{{E}_{b}}))}}\end{eqnarray*}$$
• 2.  
  Broken Power-law:

  $$\begin{eqnarray*}\begin{array}{ccccccccccccccc} F{{(E)}^{{\rm BPL}}} & = & {{N}_{{}^\circ }}{{(E/{{E}_{b}})}^{-{{{\Gamma }}_{1}}}}\;{\rm for}\;E\;\lt \;{{E}_{b}} \\ {} & = & {{N}_{{}^\circ }}{{(E/{{E}_{b}})}^{-{{{\Gamma }}_{2}}}}\;{\rm for}\;E\;\gt \;{{E}_{b}} \\ \end{array}\end{eqnarray*}$$

Here, F(E)^LP and F(E)^BPL are representing the flux for the LP- and BPL-type spectrum, respectively. The spectrum breaks or turns at a certain energy of E_b, where the spectral index changes from ${{{\Gamma }}_{1}}$ to ${{{\Gamma }}_{2}}$. N$_{{}^\circ }$ is the normalization parameter.

Using different functions in fitting the spectrum of ${\rm G349}{\rm .7+0}{\rm .2},$ the likelihood ratio, TS, was used as a measure of the improvement of the likelihood fit with respect to the simple PL. The results to the spectral fits are summarized in Table 2, where the highest TS value (168) was found for the BPL fit.

Table 2.  Spectral Fit Parameters for PL, LP, and BPL between 200 MeV and 300 GeV, Assuming G349.7+0.2 to be a Point-like Source

Spectral	Photon Flux	${{{\Gamma }}_{1}}$	${{{\Gamma }}_{2}}$	Eb	TS
Model	(10−8 ph cm−2 s−1)			(MeV)
PL	3.59 ± 0.21	2.00 ± 0.03	⋯	⋯	137
LP	1.66 ± 0.13	1.67	0.19 ± 0.31	1310	151
BPL	1.24 ± 0.14	4.98	2.28 ± 0.04	568	168

Download table as:  ASCIITypeset image

The PL resulted in a spectral index of ${{{\Gamma }}_{1}}$ = 2.00 ± 0.03, which is in agreement with the best-fit PL index value given for 2FGL J1718.1–3725 in the second Fermi-LAT catalog (∼1.98), Nolan et al. (2012). This result also matches the results obtained by Castro & Slane (2010). The best-fit parameters for the BPL fit are N$_{{}^\circ }$ = (2.30 ± 0.21) × 10⁻¹¹ MeV⁻¹ cm⁻² s⁻¹, ${{{\Gamma }}_{1}}\sim$ 4.98, and ${{{\Gamma }}_{2}}$ = 2.28 ± 0.04, where the given uncertainties are statistical. The total energy flux was found to be (3.70 ± 0.10) × 10⁻¹¹ erg cm⁻² s⁻¹.

We also calculated the systematic errors that stem from the uncertainties in the Galactic diffuse background intensity, where we followed the methods described in Abdo et al. (2009) and Castro et al. (2013). We varied the normalization value of the Galactic background by ±6% from the best-fit value and used these new frozen values of the normalization parameter to recalculate the SED of G349.7+0.2. The systematic errors on the SED are shown in Figure 4 with red solid lines on top of the statistical errors given in black lines.

To investigate the morphology of G349.7+0.2, we created a 2° × 2° TS map of G349.7+0.2 and its neighborhood with a bin size of 001 × 001. The TS map shown in Figure 3 was produced with pointlike using a background model file, which contained all of the point-like sources and diffuse sources, but excluded G349.7+0.2 from the model. Thus, it shows the TS distribution of gamma-rays originating dominantly from G349.7+0.2. In Figure 3, the blue contours represent the Suzaku combined XIS image in the 0.3–10.0 keV energy band (from Figure 1) and the black cross and circle represent the location and its statistical error of the GeV source from the second Fermi-LAT catalog corresponding to G349.7+0.2, respectively. The best-fit location of the gamma-ray emission is slightly offset from the center of the X-ray remnant, as shown in Figure 3 by the red cross and circle, and it is closer to the southeast region of the SNR.

Zoom In Zoom Out Reset image size

To search for the energy-dependent morphology, we split the data set into two energy ranges (200 MeV−1 GeV and 1–300 GeV) and analyzed the data in each range. We found no significant gamma-ray excess at the location of G349.7+0.2 for the energy range between 200 MeV and 1 GeV, but ${\rm G349}{\rm .7+0}{\rm .2}$ was detected in the higher-energy range of 1–300 GeV with a significance of ∼11σ using the BPL spectral model.

Assuming that G349.7+0.2 has a PL/BPL-type spectrum, we checked the extension of G349.7+0.2 using the TS of the extension (TS_ext) parameter. TS_ext is the the likelihood ratio comparing the likelihood for being a point-like source to a likelihood for an existing extension. The pointlike code calculates the TS_ext by fitting a source first with a disk template and then as a point-like source. In pointlike, the extended source detection threshold is defined as the source flux, where the TS_ext value averaged over many statistical realizations is 16 (Lande et al. 2012). Simulation studies showed that to resolve a disk-like extension with a radius (r) of 01 at the detection threshold, the source with a spectral index value of 2.0 and 2.5 must have a minimum flux of 3 × 10⁻⁷ ph cm⁻² s⁻¹ and 2 × 10⁻⁶ ph cm⁻² s⁻¹, respectively (Lande et al. 2012).

For G349.7+0.2 having either a PL- or a BPL-type spectrum, the TS_ext value was found to be about zero, after a disk template fitting, where the fit parameters are shown in Table 3. This indicates that a disk-like extension with r ∼ 01 could not be resolved at the integrated flux level of 3.59 × 10⁻⁸ ph cm⁻² s⁻¹ for a PL-type spectrum with a spectral index of 2.0 and of 1.24 × 10⁻⁸ ph cm⁻² s⁻¹ for the BPL-type spectrum with a spectral index of 2.3 (${{{\Gamma }}_{2}}$), respectively.

Zoom In Zoom Out Reset image size

We first looked for long-term variability in the light curve of G349.7+0.2 by taking data from the circular region of 1° around the best-fit position. Figure 5 shows the one- month binned light curve obtained after applying Fermi-LAT aperture photometry, where we checked for possible variations in the flux levels. If any or some of the flux data points are above 3σ, this would be an indication of a significant variability (e.g., a flare) in the circular ROI. In Figure 5, most of the flux data points remain within the 1σ and 3σ bands. Although one flux data point is above 3σ, we note that its error bars are large, showing a large statistical uncertainty in this time interval. Thus, we conclude that there is no long-term variability observed in the close neighborhood of G349.7+0.2.

Zoom In Zoom Out Reset image size

We have also checked if the spectral shape of G349.7+0.2 fits to the standard spectrum of a pulsar, Power Law with Exponential Cutoff (PLEC). The best-fit cutoff energy was found to be 24.60 ± 4.95 GeV, which is an order of magnitude away from the range of typical pulsar cutoff energies (Abdo et al. 2010d). The PLEC fit did not show a significant improvement over the PL, BPL, and LP spectral fits.

We used the CO data taken by the 12 m telescope of the NRAO (Reynoso & Magnum 2001) with a good angular resolution (54'') to understand the molecular environment around G349.7+0.2. We integrated the CO spectrum over the whole range of velocities from −137 to +57 km s⁻¹ to obtain the velocity integrated CO intensity (W_CO). The W_CO averaged over the region with a radius of 0025 covering the X-ray remnant was found to be ∼150 K km s⁻¹. Using the CO-to-H₂ conversion factor of X = 1.8 × 10²⁰ cm⁻² K⁻¹ km s⁻¹ (Dame et al. 2001), we found N(H₂) = 2.7 × 10²² cm⁻². Using the N(H i) value found by Tian & Leahy (2014), we calculated the total column density as N(total) = N(H i) + 2N(H₂) = 7.12 × 10²² cm⁻². Lazendic et al. (2005) calculated the total column density as 7.3 × 10²² cm⁻² adding the contributions from H i, 2.8 × 10²² cm⁻² assuming T_s = 140 K, and H₂, 4.6 × 10²² cm⁻².

Figure 6 shows the CO map produced in the velocity range of [+14, +20] km s⁻¹, where white contours represent the TS values shown in Figure 3. Also shown with black plus markers are the five masers measured by Frail et al. (1996), which are located inside of the blue Suzaku X-ray contours of 540 cts.

Zoom In Zoom Out Reset image size

The progenitor of the G349.7+0.2 has exploded within a large ringlike CO shell (seen as the color-coded arc structure in Figure 6), a remnant of an older SN explosion containing denser clumps of MCs. The explosion probably happened within the inter-clump gas (Reynoso & Magnum 2001) of the CO shell, at a location between the edge of the CO shell and an MC clump inside this shell. In our interpretation, we assume that while the SNR's shell was expanding, a small section of it broke out of the CO shell boundary and entered the rarefied interstellar medium (ISM) moving in a direction toward the observer.

The X-ray morphology of G349.7+0.2 is similar to the radio image consisting of a two-ring structure, where one ring is smaller and brighter than the other. The smaller ring might be the part of the SNR shock front entering into a higher density region, an MC clump called "cloud1" by Reynoso & Magnum (2001), and interacting with it. The shock of ${\rm G349}{\rm .7+0}{\rm .2}$ entering "cloud1" was confirmed to be of type C, which allows OH 1720 MHz masers to form (Lockett et al. 1999). The larger ring structure represents the SNR shock expanding into the relatively less dense medium of the CO shell.

Therefore, we calculated the proton density for two regions: one having a radius of 0015 corresponding to the inner smaller ring, and another one with a radius of 0022 representing the larger outer ring. The smaller ring radius was chosen such that it covers the highest number of X-ray counts (shown with the innermost blue contour for $\gt$540 cts in Figure 6). This region coincides mostly with the E region of the SNR selected in the X-ray analysis (Section 2.1.3). The larger ring was chosen to include the X-ray contour representing the X-ray counts higher than 365, where the W region is also partially covered. Both rings are within the 10σ contours of GeV gamma-ray TS distribution. To calculate the average proton density within the regions surrounded by these two rings, we first estimated the size of the emission regions to be ∼3 and ∼4.5 pc using the distance to the SNR to be 11.5 kpc. Then assuming a spherical geometry of the emission region, the proton density (n) corresponding to the inner and outer ring region was found to be 81 and 35 cm⁻³, respectively.

To understand whether the leptonic or hadronic scenario is dominating the SED of G349.7+0.2, we started the multi-wavelength modeling assuming a BPL-type of spectra for both the electron and proton distributions inside the source:

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \frac{d{{N}_{e,p}}}{dE} & = & N_{e,p}^{1}\;E_{e,p}^{-\alpha }\;{\rm for}\;{{E}_{e,p}}\;\lt \;{{E}_{{\rm br}}} \\ {} & = & N_{e,p}^{2}E_{e,p}^{-\beta }{\rm exp} \left( -\frac{{{E}_{e,p}}}{{{E}_{c}}} \right) \\ {} & {} & \quad {\rm for}\;{{E}_{{\rm br}}}\leqslant {{E}_{e,p}}\leqslant {{E}_{c}}, \\ \end{array}\end{eqnarray} \tag{ 1 }$$

where ${{E}_{e,p}}$ is the electron or proton energy and E_br is the spectral break energy with the spectral index changing from α to β. E_c is the maximum energy of electrons or protons. $N_{e,p}^{1}$ and $N_{e,p}^{2}$ are normalization constants.

For the multi-wavelength modeling, we used the radio data (Shaver et al. 1985; Whiteoak & Green 1996; Green 2009) and the gamma-ray data from Fermi-LAT and H.E.S.S. The H.E.S.S. flux data points (Abramowski et al. 2014) are shown in magenta in Figure 7. During the fitting procedure, we took E_c as 10 TeV for protons and 2 TeV for electrons. The value of E_c was chosen such that it would not cause the spectra of bremsstrahlung and inverse-Compton (IC)  to overestimate the TeV fluxes, when they are considered the dominant models for explaining the data. Similarly, E_c for protons was chosen such that the hadronic model does not overestimate the TeV fluxes.

Zoom In Zoom Out Reset image size

We first considered a purely leptonic model to explain the observed fluxes from radio wavelengths to gamma-rays. In the case of the leptonic model, bremsstrahlung and IC emission processes (Blumenthal & Gould 1970) were used to explain the GeV–TeV data. The total energetics of the electrons and protons need to be known to estimate the electron to proton ratio, defined as the ratio of dNe/dE and dNp/dE at E = 10 GeV, and the different model contributions. We first estimated the bremsstrahlung spectrum using a fixed value of ambient proton density. Parameters of the leptonic model were then adjusted in such a way that the bremsstrahlung spectrum could explain the observed fluxes at GeV–TeV energies. These parameters were then fixed and they were used in computing the IC and synchrotron spectra of the leptonic model. To explain the observed radio fluxes, the magnetic field was adjusted to get the radio synchrotron spectrum. Using the same ambient proton density, we estimated the hadronic contribution to the gamma-ray spectrum resulting from the decay of ${{\pi }^{{}^\circ }}$ (Kelner et al. 2006). Second, we repeated a similar procedure for the case of the IC process as the dominant model explaining the observed fluxes at GeV–TeV energies. At last, we considered the ${{\pi }^{{}^\circ }}$-decay model dominating among all of the processes, which is a purely hadronic model (Kelner et al. 2006). The parameters used for all three dominating processes are shown in Table 4 for the two ambient proton densities of 35 and 81 cm⁻³ derived in Section 2.2.4.

Table 4.  Parameters of the Multi-wavelength Model

Model	Parameters					Energetics
	α	β	${{k}_{e}}/{{k}_{p}}$	n	B	Leptonic	Hadronic
				(cm−3)	(μG)	(1049 ergs)	(1049 ergs)
For n = 35 cm−3
Bremsstrahlung	1.8	2.4	1.0	35	20	1.25	0.74
IC	1.8	2.6	1.0	0.1	3.0	27.5	13.1
${{\pi }^{{}^\circ }}$-decay	1.9	2.4	0.01	35	100	0.13	7.44
For n = 81 cm−3
Bremsstrahlung	1.8	2.4	1.0	81	40	0.55	0.31
IC	1.8	2.6	1.0	0.1	3	27.5	13.1
${{\pi }^{{}^\circ }}$-decay	1.9	2.4	0.01	81	150	0.07	3.13

Note. For bremsstrahlung and ${{\pi }^{{}^\circ }}$-decay dominated models, the ambient proton density is considered to be 35 and 81 cm⁻³.

Download table as:  ASCIITypeset image

Since the ambient proton density of 81 cm⁻³ represents the interaction region corresponding to the smaller and brighter SNR shell, including the five OH masers, the modeling results for this density will be important in terms of giving us clues about the proton spectrum. The best-fit parameters for the proton spectrum were obtained by a ${{\chi }^{2}}$-fitting procedure to the flux points. The dominating models of bremsstrahlung, IC, and ${{\pi }^{{}^\circ }}$-decay for the ambient proton density of 81 cm⁻³ are shown in the top, middle, and bottom panels of Figure 7, respectively. The estimated parameters are $\alpha =1.89\;\pm \;0.14$, $\beta =2.42\;\pm \;0.03$, and E_br=12 GeV. The ${{\chi }^{2}}$/dof is estimated to be $\simeq$0.53. The best-fit gamma-ray spectrum resulting from the decay of ${{\pi }^{{}^\circ }}$ is shown in the bottom panel of Figure 7 with the red solid line. The estimated total energy can be written as W_p $\simeq$ 3.13 × 10$^{49}\left( 81\;{\rm c}{{{\rm m}}^{-3}}/n \right)$ erg, where n is the effective gas number density for the p-p collision.