FERMI LARGE AREA TELESCOPE OBSERVATIONS OF THE SUPERNOVA REMNANT W28 (G6.4−0.1)

Figure 1 shows a smoothed count map in the 2–10 GeV energy band in a 10° × 10° region around W28. Figure 2 gives comparisons with images of other wavebands in close-up view. Black contours indicate the HESS significance map in panel (a), the CO (J = 1–0) line intensity taken by NANTEN for the velocity range from 0 km s⁻¹ to 20 km s⁻¹ (corresponding to kinematic distances of approximately 0 to 4 kpc) in panel (b), and a Very Large Array (VLA) image in panel (c). Correlations between GeV gamma rays observed by LAT and some of the HESS sources are evident. Gamma rays are also bright in the brightest spots in the CO contours, and bright regions in the gamma-ray image extensively overlap with bright regions in the VLA contours. There are two LAT sources in the vicinity of W28 in the 1FGL catalog (Abdo et al. 2010d): 1FGL J1801.3−2322c and 1FGL J1800.5−2359c. No obvious pulsations of gamma rays are found in these sources. Hereafter, we refer to 1FGL J1801.3−2322c as Source N and to 1FGL J1800.5−2359c as Source S.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to quantitatively evaluate the extension and location of these two sources, we apply the maximum likelihood tool, gtlike, which is publicly available as part of the Fermi ScienceTools. The likelihood is the product of the probability of observing the gamma-ray counts of each spatial and energy bin given the emission model, and the best parameter values are estimated by maximizing the likelihood of the data given the model (Mattox et al. 1996). The probability density function for the likelihood analysis included (1) individual sources detected in the preliminary LAT one-year catalog, (2) the Galactic diffuse emission resulting from cosmic-ray interactions with the interstellar medium and radiation based on the LAT standard diffuse background model gll_iem_v02 available from FSSC,⁵² and (3) the isotropic component to represent extragalactic and residual cosmic-ray backgrounds using the isotropic spectral template isotropic_iem_v02 from FSSC. Note that we make energy-dependent corrections of the Galactic diffuse model by multiplying a power-law function with the spectral index free to vary in the fit. This correction gives better spectral fits by taking into account local systematic discrepancies between the data and the Galactic diffuse model. The region of interest for the binned maximum likelihood analysis was a square region of 20°×20° centered on W28 with a pixel size of 01. The instrument response functions (IRFs) used in our work were the "Pass 6 v3" (P6_V3) IRFs, which were developed following launch to address gamma-ray detection inefficiencies that are correlated with background rates.

Using the tools described above, we investigated spatial extensions of the two LAT W28 sources. Here, we used only events above 2.15 GeV to take advantage of the narrower PSF in the higher energy band. Before investigating the detailed extension, we first determined the flux and spectral index for all components except for Source N and Source S. In this process, the normalization of the Galactic diffuse emission and the flux and spectral index of the power-law model for the sources within 5°of Source N were set free to account for the effects of sources around W28 on the background flux in the fit. The flux and spectral index, except for those of Source N and Source S, are fixed hereafter. We used a radially symmetric uniform disk to evaluate the extension. We varied the radius and location of the disk while holding the position of Source S fixed at the catalog position and evaluated the resulting maximum likelihood value (L_ex) with respect to the maximum likelihood for no source hypothesis (L₀) and with the point source hypothesis (L_ps). A two-source hypothesis (L_2s) is also evaluated by scanning the positions of the two point sources.

The best likelihood ratio −2ln(L_ps/L_ex) ≈ 76 is obtained for Source N with a disk radius of σ = 039, which rejects a point source hypothesis at more than 8σ. Note that we validated this procedure by applying this method on a nearby gamma-ray pulsar, PSR J1809−2332 (Abdo et al. 2010c; ∼2° away from W28), where we find the extension to be consistent with a point source. In addition, the best likelihood ratio −2ln(L₀/L_ex) ≈ 336 for the disk shape is substantially better than that for the two point sources, −2ln(L₀/L_2s) ≈ 276, where the positions of the two sources were free in the optimization. Therefore, we use the best-fit disk model for further analysis of Source N. The best location of the disk model is found at (R.A., decl.) = (18^h01^m21^s, −23°26'26'') with an error radius of 003 at 68% confidence level.

The extension of Source S is also investigated using the same procedure as above. We did not find significant extension. An upper limit on the size of the gamma-ray emission was obtained by investigating the decrease of the likelihood with increasing source size. Under the assumption of a disk shape, the upper limit amounts to 16' at the 68% confidence level. We treat Source S as a point source hereafter. The best-fit location of Source S was estimated to be (R.A., decl.) = (18^h00^m59^s, −24°11'31'') using gamma rays above 2.15 GeV. This location is 13' away from the location of 1FGL J1800.5−2359c using gamma rays above 0.1 GeV and assuming that 1FGL J1801.3−2322c is a point source. In order to investigate the origin of the position difference, we performed localization with the same energy cut as the catalog analysis assuming a point source hypothesis and the extended source hypothesis. We found that the localization with a point source hypothesis is consistent with the catalog position while the extended source hypothesis yields the location consistent with our localization with the E > 2.15 GeV cut. From this comparison, we conclude that the localization difference is mostly due to the extension of Source N. We use the average of the two localizations, (R.A., decl.) = (18^h00^m47^s, −24°05'45''), for the following analyses, with a systematic uncertainty of half of the positional difference, ∼6'.

We find no strong gamma-ray emission in the LAT energy band at the locations of HESS J1800−240A and C. We evaluated the likelihood of point sources at these locations and found likelihood ratios −2ln(L₀/L_ps) ≈ 10 (corresponding to a ∼3σ) for HESS J1800−240A and −2ln(L₀/L_ps) ≈ 9 (corresponding to a ∼3σ) for HESS J1800−240C, which confirms no significant emission at these locations.

We used the maximum likelihood fit tool, gtlike, for the spectral analysis of LAT sources in this region. Flux points were obtained by performing the maximum likelihood analysis in each energy bin. We used eight logarithmically spaced energy bins from 215 MeV to 100 GeV for Source N. We did not divide the energy range 0.2 GeV–100 GeV into bins for Source S due to poor statistics. Figure 3 shows the resulting spectral energy distributions (SEDs) for (a) Source N and (b) Source S. The 68% confidence region is illustrated for Source S up to ∼20 GeV. We do not have a sufficient number of events above ∼20 GeV to constrain the spectrum. Upper limits at 90% confidence level are calculated assuming a photon index of 2 if the detection is not significant in an energy bin, i.e., the likelihood ratio with respect to no source is less than 9. Note that the value of the spectral index has negligible effect on the upper limits. We take into account systematic errors due to uncertainties of the extension, the Galactic diffuse model, the LAT effective area, and the effect of the nearby gamma-ray pulsar PSR J1809−2332. Systematic errors associated with the extension are estimated by varying the source size by ±1σ. We also evaluated the effect of the shape by comparing the disk shape used in this paper and a circular two-dimensional Gaussian. Systematic errors due to the Galactic diffuse model are estimated by using the residual gamma-ray data with respect to the best-fit model in the region where no LAT source is present, specifically, l = 23–47 and b = −125–075. The observed residual is energy dependent and can be modeled as ∼ (116 (E/1 GeV)$^{6.45 \times 10^{-2}}-100$)% of the total Galactic diffuse flux. The normalization of the Galactic diffuse model is adjusted according to the above equation to estimate the systematic error in the source flux. We evaluated the position dependence of this residual and found a dispersion of ∼5%, which is not a large effect. Uncertainties of the Galactic diffuse model make a dominant contribution to the systematic errors in all the energy bins. Uncertainties of the LAT effective area are 10% at 100 MeV, decreasing to 5% at 500 MeV, and increasing to 20% at 10 GeV and above (Rando et al. 2009). The effect of the nearby gamma-ray pulsar was estimated by varying the flux by ±1σ.

Zoom In Zoom Out Reset image size

We evaluated a possibility of a spectral break in the LAT energy band by comparing the likelihood of the spectral fit for the LAT data between a simple power law and a broken power law as a spectral model of Source N and Source S. The fits yield the likelihood ratio −2ln(L_PL/L_BPL) ≈ 62 for Source N, where L_PL and L_BPL are the likelihoods for the simple power-law model and the broken power-law model, respectively. The likelihood ratio slightly decreases to at least ∼48 (corresponding to 6.6σ with two degrees of freedom) in the worst case, accounting for 1σ systematic uncertainties. Thus we conclude that Source N has a spectral break at 1.0 ± 0.2 GeV, where the error is dominated by statistics. Photon indices are 2.09 ± 0.08 (stat) ± 0.28 (sys) below the break and 2.74 ± 0.06 (stat) ± 0.09 (sys) above the break. Note that the photon index above the break is consistent with the HESS measurement of 2.66 ± 0.27 (Aharonian et al. 2008). On the other hand, we do not find any evidence of a spectral break for Source S. The photon index is found to be 2.19 ± 0.14 (stat) ± 0.41 (sys).

We placed upper limits on the gamma-ray flux in the LAT band at the positions of HESS J1800−240A and HESS J1800−240C. The upper limits are compared with the HESS spectra in Figures 4(a) and (b). The gamma-ray upper limits in the GeV band at the location of HESS J1800−240A appear to be inconsistent if the source spectrum is a simple power law. We evaluated the possibility of a break between the LAT and HESS energy band by comparing the likelihood of the spectral fit for the LAT data between a broken power law and a simple power law as a spectral model of this source. We fixed the spectral index and the flux values using the HESS measurements for the simple power-law model, while we fixed the spectral index above the break and the flux at the same value as the simple power law and varied the break energy and the spectral index below the break for the broken power-law model. We used the parameter values at 1σ away from the best-fit values from HESS in the direction where we expect less flux in the LAT energy band, i.e., harder spectral index and lower flux than the best-fit values. The fits yield the likelihood ratio −2ln(L_PL/L_BPL) ≈ 35 (corresponding to a 5.5σ significance with two degrees of freedom) for a break energy between 10 GeV and 100 GeV. The likelihood ratio decreases to ∼8, corresponding to a ∼2.4σ significance, if we use parameter values 2σ away from the best-fit values from HESS. We conclude that HESS J1800−240A is likely to have a spectral break between 10 GeV and 100 GeV.

Zoom In Zoom Out Reset image size

1FGL J1801.3−2322c (Source N) is found to be extended and positionally coincident with HESS J1801−233. Given the clear spatial match between the TeV source and the molecular cloud interacting SNR W28 and a smoothly connected spectrum in the GeV–TeV band, we assume that the bulk of GeV and TeV gamma-ray emission comes from the molecular cloud illuminated by particles accelerated in supernova shocks. Below, we adopt the simplest assumption that GeV and TeV gamma rays are emitted by a population of accelerated protons and electrons distributed in the same region characterized by constant density and magnetic field strength. It should be noted that the imaging resolution of the LAT is not good enough to differentiate the GeV emission region from the TeV region. We assume that the injected electrons have the same momentum distribution as the protons. This assumption implies a break in the particle momentum spectrum because the spectral index of the radio data, corresponding to lower particle momenta, is much harder than that for the gamma-rays, which correspond to higher particle momenta. Therefore, we use a broken power law to describe the particle momentum distribution in the region. Electrons suffer energy losses due to ionization (or Coulomb scattering), bremsstrahlung, synchrotron processes, and IC scattering. The modification of the electron spectral distribution due to such losses was calculated according to Atoyan et al. (1995), where electrons are assumed to be injected at t = 0 from an impulsive source. Since diffusive shock acceleration theory generally predicts particle accelerations in the Sedov phase with a typical duration of 10³–10⁴ yr, the assumption of an impulsive source would be a good approximation for SNRs with the age of 10⁴–10⁵ yr. We adopt 4 × 10⁴ yr as the age of W28 in this modeling. Note that here we consider the minimum momenta of protons and electrons to be 100 MeV c⁻¹ since the details of the proton/electron injection process are poorly known. The gamma-ray spectrum from π⁰ decay produced by the interaction of protons with ambient hydrogen is scaled by a factor of 1.84 to account for helium and heavy nuclei in the target material and cosmic-ray composition (Mori 2009).

First, we consider a π⁰-decay model to account for the broadband gamma-ray spectrum. Indeed, the hadronic scenario gives the most satisfactory explanation for the GeV gamma rays observed in other SNRs interacting with molecular gas such as W51C (Abdo et al. 2009) and W44 (Abdo et al. 2010a). The number index of protons in the high-energy regime is constrained to be α_H ≈ 2.7 from the gamma-ray spectral slope. The observed gamma-ray luminosity requires the gas density to be much larger than $\bar{n}_{\rm H}\approx 1\ {\rm cm}^{-3}$ averaged over the entire SNR shell in order not to exceed the typical kinetic energy of a supernova explosion (∼10⁵¹ erg). The resulting total proton energy, $W_{p}\sim 1.3\times 10^{49}\times (10^2\ {\rm cm}^{-3}/\bar{n}_{\rm H})\times (d/2\ {\rm kpc})^{2}$ erg, is less than 10% of the typical kinetic energy of supernova explosions and quite reasonable. Note that W_p is not the total energy of accelerated protons but that of the in situ protons in the molecular clouds. Using the parameters summarized in Table 1, we calculated radiation model curves as shown in Figure 5(a). As described above, the spectral index of the proton momentum below the break is determined to be α_L ≈ 1.7 by modeling the radio spectrum as synchrotron radiation by relativistic electrons because we assume that protons and electrons have identical injection spectra. It is difficult to derive the break point of the proton momentum spectrum from the break of the gamma-ray spectrum since it lies in the region where we expect a gamma-ray spectral curvature due to the kinematics of π⁰ production and decays. Because of this, the gamma-ray spectrum gives an upper bound for the momentum break at ∼5 GeV c⁻¹. The momentum break cannot be lower than ∼2 GeV c⁻¹ to avoid conflict with the radio data. Here, we adopt 2 GeV c⁻¹. The magnetic field strength is constrained to be B ∼ 160 μG for an electron-to-proton ratio of K_ep = 0.01, which is the ratio found in the local cosmic-ray abundance. Here, the ratio is defined at a particle momentum of 1 GeV c⁻¹. This choice of K_ep results in minor contribution from the electron bremsstrahlung emission in the gamma-ray band.

Zoom In Zoom Out Reset image size

Table 1. Model Parameters for the Fermi LAT Sources Around W28

Model	Kepa	αLb	pbc	αHd	B	$\bar{n}_{\rm H}$e	Wpf	Wef
			(GeV c−1)		(μG)	(cm−3)	(1049 erg)	(1049 erg)
Source N
(a) Pion	0.01	1.7	2	2.7	160	100	1.3	1.9 × 10−2
(b) Bremsstrahlung	1	1.7	1	2.7	4	5	1.9	4.9
(c) Inverse Comptong	1	1.7	5	3.6	0.6	0.02	54	90
Source S
(a) Pion	0.01	1.7	2	2.4	1.2 × 103	103	1.5 × 10−2	6.5 × 10−5
(b) Bremsstrahlung	1	1.7	1	2.4	4	103	1.1 × 10−3	2.1 × 10−3
(c) Inverse Comptong	1	1.7	15	3.2	2	0.05	1.7	2.3

Notes. ^aThe ratio of electron and proton distribution functions at 1 GeV c⁻¹. ^bA momentum distribution of particles is assumed to be a broken power law, where the indices and the break momentum are identical for both accelerated protons and electrons. The α_L is the photon index below the momentum break. ^cThe p_b is the momentum break for particle distribution. ^dThe photon index for broken power-law functions above the momentum break. ^eAverage hydrogen number density of ambient medium. ^fThe distance is assumed to be 2 kpc. The total energy is calculated for particles to be > 100 MeV c⁻¹. ^gSeed photons for IC scattering of electrons include the CMB, two infrared (T_IR = 29, 4.9 × 10² K, U_IR = 0.29, 5.3 × 10⁻² eV cm⁻³, respectively), and two optical components (T_opt = 3.6 × 10³, 1.0 × 10⁴ K, U_opt = 0.37, 0.13 eV cm⁻³, respectively) in the vicinity of W28, assuming a distance of 2 kpc.

Download table as:  ASCIITypeset image

On the other hand, it is difficult to model the GeV–TeV spectrum by the electron bremsstrahlung component alone because the break in the electron spectrum corresponding to the gamma-ray spectrum will appear in the radio data as shown in Figure 5(b) although it might be contaminated by a thermal component (Dubner et al. 2000). Moreover, the magnetic field strength is constrained to be less than B ∼ 4 μG for the age of 4 × 10⁴ yr since no apparent synchrotron cooling is observed at TeV energies. The low magnetic field requires the gas density to be smaller than $\bar{n}_{\rm H}\sim 5\ {\rm cm}^{-3}$ from the flux ratio of the synchrotron component in the radio band and the bremsstrahlung in the gamma-ray band. These low values for the magnetic field and the gas density are in disagreement with our assumption that the gamma-ray emission comes from the dense molecular cloud.

The gamma-ray spectrum is formally reproduced by IC emission from accelerated electrons in Figure 5(c). The interstellar radiation field for the IC scattering (see Table 1) is comprised of infrared, optical, and the cosmic microwave background (CMB). The infrared and optical components are taken from the interstellar radiation field at the location of W28 in the GALPROP code (Porter et al. 2008). Since the spectral shape of the non-CMB component is very complex, it is approximated by two infrared and two optical blackbody components. The flux ratio between the IC and the synchrotron components constrains the magnetic field to be less than 0.6 μG. The total energy in the electrons is calculated to be W_e ≈ 9.0 × 10⁵⁰ × (d/2 kpc)² erg for an energy density of ∼1.8 eV cm ⁻³ for the interstellar radiation field. This model requires a rather low gas density of $\bar{n}_{\rm H}\sim 2 \times 10^{-2}\, {\rm cm}^{-3}$ to suppress the electron bremsstrahlung, which contradicts our assumption.

Observations with Fermi have demonstrated that bright, extended gamma-ray sources coincident with middle-aged SNRs interacting with molecular clouds, such as W44 (Abdo et al. 2010a), W51C (Abdo et al. 2009), and IC 443 (Abdo et al. 2010b), exhibit spectral breaks from 9 GeV c⁻¹ to a few tens of GeV c⁻¹ in the proton momentum spectrum. Our observations of 1FGL J1801.3−2322c in the vicinity of W28, combined with the radio data, constrain the proton momentum break to be in the range of 2–5 GeV c⁻¹. The observed energy distribution of relativistic particles (most likely protons) could be influenced greatly by the diffusive transport of particles. If so, the relationship between the observed particle spectrum and the acceleration spectrum should be rather complicated. The steep particle spectrum, α ∼ 2.7, deduced for 1FGL J1801.3−2322c does not necessarily represent the acceleration index. Gabici & Aharonian (2007) discussed the time evolution of non-thermal emission from molecular clouds illuminated by cosmic rays from a nearby SNR and predicted a steep gamma-ray spectrum for an old SNR due to the energy-dependent diffusion of cosmic rays. More detailed studies of the properties of interactions between SNRs and molecular clouds are required for a comprehensive description of the differences in the break momenta and photon indices above the break among the above SNRs.

1FGL J1800.5−2359c (Source S) was found to have no significant extension and spatially coincides with the TeV source HESS J1800−240B, which is associated with molecular clouds that contain the ultra-compact H ii region W28A2 (G5.89−0.39). The observed GeV–TeV spectrum for Source S can be formally described by π⁰-dominated, bremsstrahlung-dominated, and IC-dominated models in a similar manner to Source N, as shown in Figure 6. Free–free emission (Gómez et al. 1991) is responsible for the radio spectrum of W28A2, which provides only upper limits in the SED. Thus, values of B and $\bar{n}_{\rm H}$ are less constrained by pion and bremsstrahlung models. We adopt $\bar{n}_{\rm H}= 10^3$ cm⁻³ (Aharonian et al. 2008) for both models. From the positional coincidence, W28A2 can be invoked as a possible source of high-energy particles responsible for the gamma-ray emission. W28A2 exhibits very energetic bipolar molecular outflows (Harvey & Forceille 1988; Acord et al. 1997; Sollins et al. 2004), which would arise from the accretion of matter by a stellar progenitor. Klaassen et al. (2006) estimated the total kinetic energy of the outflow of W28A2 to be 3.5 × 10⁴⁶ erg, which requires very high matter density to account for the observed gamma-ray flux in both pion and bremsstrahlung models. If the matter density is as high as ∼10⁷ cm⁻³ as suggested by Klaassen et al. (2006), the resulting total energy would be less than a few percent of the total kinetic energy of the outflow of W28A2.

Zoom In Zoom Out Reset image size

Cosmic rays that escaped from the SNR W28 in earlier epochs would be another possible explanation given the energetics. Note that Source S would lie at a projected distance of ∼10 × (d/2 kpc) pc from the southern circular boundary of W28. If we assume that cosmic rays are uniformly radiated by the SNR, the cosmic-ray flux bombarding the molecular cloud from the SNR can be scaled by the fraction of the solid angle of the molecular cloud, ∼8 × 10⁻³. Using this scaling factor, the total energies in the entirety of SNR W28 could be calculated from the total energy observed in this source to be ∼2 × 10⁴⁹ erg for protons in the π⁰ model. This is a reasonable value for an SNR. The energy density of such cosmic rays is enhanced to be at least ∼2 eV cm⁻³ at the molecular clouds under the assumption that all of the associated clouds (M ∼ 4 × 10⁴ M_☉ at 2 kpc) are interacting with energetic particles. The molecular clouds can accumulate those cosmic rays in the past since the diffusion coefficient is expected to be low in dense environments (Aharonian & Atoyan 1996).

While we detected LAT sources that are spatially coincident with HESS J1801−233 and HESS J1800−230B, no strong LAT counterpart can be found at the positions of HESS J1800−230A and HESS J1800−230C. We found evidence of a spectral break between 10 GeV and 100 GeV for HESS J1800−230A. The velocity with respect to the local standard of rest of the molecular emission positionally coincident with this source peaks at ∼16 km s⁻¹, corresponding to a distance of ∼4 kpc (Aharonian et al. 2008). However, the velocity difference between this cloud and the one associated with HESS J1800−230B is only ∼10 km s⁻¹ and can be attributed to the proper motion of the cloud complex that contains both clouds. If the distance to the cloud associated with HESS J1800−230A is 4 kpc, then it is unlikely that this source is associated with W28 since the distance to W28 is ∼2 kpc. The H ii regions G6.1−0.6 and/or G6.225−0.569 could be a source of accelerated particles. If the cloud associated with HESS J1800−230A is part of the cloud complex that includes the cloud associated with HESS J1800−230B, then the distance would be 2 kpc, which is very similar to the clouds associated with W28. In this case, particles that escaped from W28 may be responsible for the gamma-ray emissions at both HESS J1800−230A and B, and further studies of the origin of the difference in the spectral shapes of these sources may provide good constraints on the particle diffusion process.