<figure><img src="image_1.png" /></figure>

# On the Origin of the Gamma-Ray Emission toward SNR CTB 37A with Fermi-LAT

Soheila Abdollahi¹ Ⓓ, Jean Ballet² Ⓓ, Yasushi Fukazawa¹ Ⓓ, Hideaki Katagiri³ Ⓓ, and Benjamin Condon⁴ Ⓓ

¹ Department of Physical Sciences, Hiroshima University, Higashi-Hiroshima, Hiroshima 739-8526, Japan; abdollahi.soheila@gmail.com

² AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, F-91191 Gif-sur-Yvette, France

³ College of Science, Ibaraki University, 2-1-1, Bunkyo, Mito 310-8512, Japan

⁴ Centre d'Études Nucléaires de Bordeaux Gradignan, IN2P3/CNRS, Université Bordeaux1, BP120, F-33175 Gradignan Cedex, France

Received 2020 March 17; revised 2020 April 26; accepted 2020 May 8; published 2020 June 16

## Abstract

The middle-aged supernova remnant (SNR) CTB 37A is known to interact with several dense molecular clouds through the detection of shocked H₂ and OH 1720 MHz maser emission. In the present work, we use eight years of Fermi-LAT Pass 8 data, with an improved point-spread function and an increased acceptance, to perform detailed morphological and spectral studies of the γ-ray emission toward CTB 37A from 200 MeV to 200 GeV. The best fit of the source extension is obtained for a very compact Gaussian model with a significance of 5.75σ and a 68% containment radius of 0.116 ± 0.014stat ± 0.017sys above 1 GeV, which is larger than the TeV emission size. The energy spectrum is modeled as a LogParabola, resulting in a spectral index $\alpha = 1.92 \pm 0.19$ at 1 GeV and a curvature $\beta = 0.18 \pm 0.05$, which becomes softer than the TeV spectrum above 10 GeV. The SNR properties, including a dynamical age of 6000 yr, are derived assuming the Sedov phase. From the multiwavelength modeling of emission toward the remnant, we conclude that the nonthermal radio and GeV emission is mostly due to the reacceleration of preexisting cosmic rays (CRs) by radiative shocks in the adjacent clouds. Furthermore, the observational data allow us to constrain the total kinetic energy transferred to the trapped CRs in the clouds. Based on these facts, we infer a composite nature for CTB 37A to explain the broadband spectrum and to elucidate the nature of the observed γ-ray emission.

*Unified Astronomy Thesaurus concepts: Supernova remnants (1667); Galactic cosmic rays (567); Non-thermal radiation sources (1119)*

## 1. Introduction

Supernova remnants (SNRs) have been regarded as the most promising candidates for the bulk of Galactic cosmic rays (CRs). Energetic particles are produced at the shock waves associated with SNRs through the diffusive shock acceleration process (DSA; Axford et al. 1977; Krymskii 1977; Bell 1978a, 1978b; Blandford & Ostriker 1978). The main phenomenological argument in favor of this hypothesis is that SNRs are able to provide the total energy budget necessary to maintain the Galactic population of CRs, if approximately 10% of the kinetic energy released by supernova (SN) explosions can be transferred to CRs at SNR shocks (Ginzburg & Syrovatskii 1964). Moreover, the DSA mechanism can explain the hard power-law energy spectrum of CRs at their source with differential spectral index close to 2 (e.g., Fermi 1949; Bell 1978a, 1987).

The recent observations of several (young and middle-aged) SNRs in very high-energy (VHE) γ-rays (≥100 GeV) as well as nonthermal X-ray emission imply effective production of relativistic particles in support of the DSA paradigm, although the radiation mechanism responsible for the GeV/TeV emission is still under debate. This is because γ-rays can be produced by energetic hadrons (protons and nuclei) through inelastic collisions with interstellar gas, and/or by energetic electrons through nonthermal bremsstrahlung and inverse Compton (IC) scattering of ambient radiation fields. Therefore, to identify the nature of the γ-ray emission and particle species, detailed morphological and spectral studies correlated with other multiwavelength data are crucial.

Middle-aged SNRs interacting with nearby dense molecular clouds (MCs) show strong evidence of relativistic protons. The characteristic neutral pion-decay signature in the γ-ray

spectrum below ~200 MeV (often called “pion bump”) in three SNRs, IC 443, W44, and W51, associated with MCs is explicitly linked to the hadronic acceleration (Ackermann et al. 2013; Jogler & Funk 2016). Additional support for a hadronic origin of the γ-ray emission from SNRs interacting with MCs (often termed SNR/MC interactions) comes from observations of W28 (Abdo et al. 2010), W41, MSH 17–39, G337.7–0.1 (Castro et al. 2013), and G5.7–0.1 (Joubert et al. 2016). The SNR/MC interactions have been confirmed through the detection of 1720 MHz OH maser emission toward shocked regions of these remnants (e.g., Frail et al. 1996). Given the ambient gas density in the middle-aged SNRs, the total energy budget in accelerated particles is estimated to be $\approx 10^{50}$ erg.

SNR CTB 37A (also known as G348.5+0.1) was initially identified as a discrete source in radio surveys by Milne & Hill (1969). It is a middle-aged remnant, ~10⁴ yr old (Sezer et al. 2011; Pannuti et al. 2014), located in the CTB 37 complex region near two other remnants, CTB 37B (G348.7+0.3, associated with HESS J1713–381) and G348.5–0.0. The distance to the remnant has been estimated to be in the range between 6.3 and 9.5 kpc, based on H I absorption measurement along the line of sight in high-resolution Southern Galactic Plane Survey data by Tian & Leahy (2012). The distances obtained for CTB 37B and G348.5–0.0 (~13.2 kpc and ≤6.3 kpc, respectively) by those authors imply that the three remnants of the CTB 37 complex merely lie along similar lines of sight while they differ in distance. Association of CTB 37A with several nearby dense MCs has been firmly established based on observations of several OH (1720 MHz) maser spots detected toward the remnant (Frail et al. 1996) and also shocked clumps of clouds with high column
densities (Reynoso & Mangum 2000; Maxted et al. 2013; Braiding et al. 2018).

CTB 37A has been observed extensively across a wide range of energies, from radio to VHE γ-rays. Radio observations of the SNR (Milne & Dickel 1975; Dawnes 1984; Kassim et al. 1991) have revealed a shell structure with bright emission in the northern and eastern rims, and also a faint extension suggestive of a “breakout” into an inhomogeneous medium with a large-scale density gradient toward the southwest (Kassim et al. 1991). The angular size of the remnant (including the breakout) is estimated to be 19 × 16 arcmin² at $\nu = 843$ MHz by Whiteoak & Green (1996). From the same radio observations, a flux density of $S_\nu = 71$ Jy is obtained by those authors. Kassim et al. (1991) derived a typical spectral index of $\alpha_r = -0.5 \pm 0.1$ (where $S_\nu \propto \nu^{\alpha_r}$) at frequencies above 330 MHz. In the X-ray band, in addition to the extended soft thermal component dominated by emission lines of highly ionized species of Mg, Si, S, Ar, and Ca, a compact nonthermal hard X-ray emission has also been detected from the northwest of the SNR, which might be associated with the emission of a pulsar wind nebula (PWN) as suggested by Aharonian et al. (2008), Sezer et al. (2011), and Yamauchi et al. (2014). Searches for γ-ray pulsations in this region by Saz Parkinson et al. (2018) led to the discovery of pulsar PSR J1714−3830 coincident with SNR CTB 37A. X-ray and radio observations of this source indicate that CTB 37A is a mixed-morphology SNR (Aharonian et al. 2008) as characterized by a center-filled thermal X-ray emission surrounded by a shell-like radio structure. At GeV energies, it was detected for the first time as an extended source of Gaussian width 0.13 with a significance of ~4.5σ by Brandt & Fermi-LAT Collaboration (2013). In the previous study, by Castro & Slane (2010), a detailed characterization of the spectral properties was not possible due to limited photon statistics. The TeV γ-ray source HESS J1714−385 with extension of ~0.07 is positionally coincident with the SNR (Aharonian et al. 2008), though the hadronic or leptonic nature of the γ-ray emission toward this source still remains elusive.

In this paper, benefiting from a significant improvement in the LAT sensitivity implemented in Pass 8 event selection/reconstruction as well as increased photon statistics, we report a detailed analysis of 8 yr of Fermi-LAT γ-ray data around the SNR CTB 37A, and discuss the morphological and spectral characteristics of the γ-ray emission toward the remnant, which are crucial for distinguishing between hadronic and leptonic scenarios. In Section 2, observations and data reduction are briefly described. The analysis procedures and results are given in Section 3, where the spatial extension and spectrum of the remnant are explained. In Section 4, we argue that the SNR is a composite system. In addition, the crushed clouds scenario is examined to explain the multiwavelength emission toward the system through the reacceleration of trapped CRs by the radiative shocks and the subsequent adiabatic compression. A comparison to other SNRs is presented in Section 5. Finally, the conclusions are summarized in Section 6.

## 2. Observations and Data Reduction

The Fermi Large Area Telescope (LAT) is a pair-conversion detector, designed to survey the high-energy γ-ray sky in the energy range from ~20 MeV to above 300 GeV. The LAT is equipped with a tracker/converter for direction reconstruction of incident γ-rays, a CsI(Tl) crystal calorimeter for

measurement of the energy deposition, and a surrounding anticoincidence detector for rejection of the charged particle background. Full details of the LAT instrument and data processing can be found in Atwood et al. (2009), and information regarding the on-orbit calibration is given in Abdo et al. (2009a).

Recently, new event reconstruction and classification known as Pass 8 (Atwood et al. 2013), released by the Fermi-LAT collaboration, has provided substantial improvements in the instrument response functions (IRFs). Pass 8 presents significantly increased effective area (~8000 cm² on-axis above 1 GeV), an improved angular resolution as given by the point-spread function (PSF, with a 68% containment radius of ~0.8 at 1 GeV), a reduced energy dispersion (<10% between 1 GeV and 100 GeV), and a wider field of view (~2.4 sr at 1 GeV). Taking advantage of the Pass 8 data combined with the increased photon statistics collected by the LAT, more information on source extension and spectral properties will be revealed.

The LAT Pass 8 data used for the following analysis were collected in sky-survey mode during the first eight years of scientific operations, which began on 2008 August 4. The γ-rays in the 0.2–200 GeV energy range within a region of interest (ROI) with a radius of 15° centered on the position of CTB 37A are selected for the binned maximum likelihood analysis. The 200 MeV lower limit was chosen to reduce the contamination from underpredicted γ-ray emission at low energies in the Galactic plane. The event selection is based on the low background Pass 8 source event class and the corresponding IRF is P8R2_SOURCE_V6. A zenith angle cut of 90° and 105° for events below and above 1 GeV, respectively, is applied to minimize the contamination from the Earth limb (Abdo et al. 2009b).

## 3. Analysis and Results

The γ-ray data were analyzed using the LAT Science Tools software package (v11r05p02), publicly available from the Fermi Science Support Center (FSSC).<fnref n="5" /> Two different tools (pointlike and gtlike) were used to determine the morphological and spectral characteristics of the remnant, respectively, under the maximum likelihood fitting technique (Mattox et al. 1996). pointlike is an alternative package for fast binned likelihood analysis (Kerr 2010) validated by Lande et al. (2012) and is optimized to evaluate the best-fit position and extension of the source before performing a more accurate fit of the spectrum using gtlike. These tools fit a γ-ray emission model to the LAT data along with models for the instrumental, extragalactic, and Galactic components of the background. The Galactic diffuse emission is modeled using the LAT standard diffuse emission model gll_iem_v06. fits, and the residual instrumental background and extragalactic γ-ray radiation are combined into a single isotropic component with a spectral shape described by a tabulated model iso_P8R2_SOURCE_V6_v06.txt. The models and their detailed descriptions are available from the FSSC. In addition to the two aforementioned diffuse models and all background sources within 20° around CTB 37A listed in the third Fermi-LAT catalog (3FGL, Acero et al. 2015), two

<footnote n="5">⁵ The Science Tools package and supporting documents are distributed by the Fermi Science Support Center and can be accessed at http://fermi.gsfc.nasa.gov/ssc/data/analysis/software/.</footnote>
statistically significant point sources at a distance <0.5 from CTB 37A, not already detected by the 3FGL, are included in the region model. One of them, the CTB 37B source (3FHL J1714.0–3811) with a power-law spectrum, detected by the LAT at energies above 10 GeV, is taken from the third catalog of hard Fermi-LAT sources (3FHL, Ajello et al. 2017). The other one, the FL8Y J1714.8–3850 source with a LogParabola spectral model, located in the southeast of the remnant, is given by the preliminary Fermi-LAT eight-year list of sources (FL8Y).<fnref n="6" /> The detected pulsar PSR J1714-3830 is not included in the model because the ephemerides are not long enough to cover the entire LAT data set. Therefore, it is assumed that the entire γ-ray emission is originated in the SNR, which establishes an upper limit on the GeV flux of the remnant. Detailed analyses of off-pulse emission will be discussed in future work. Background sources over an area 5° larger than the ROI are contained in the model to account for the contamination from their emission.

As a first step, to assess the morphology of CTB 37A, the position and extension of the remnant, together with the position of the two nearby point sources, CTB 37B and FL8Y J1714.8–3850, were fitted in an iterative process. Moreover, in this approach, a preliminary estimate of the spectral values for sources within 2° of CTB 37A, which may affect its flux, was provided by allowing their spectral parameters (normalization, index, and cutoff energy) to vary. Only for the source 3FGL J1718.1–3825, located 0.7° away from the remnant and associated with the pulsar PSR J1718–3825, our checks showed that the pulsar is not bright enough to justify freeing the exponential index (b) of an exponentially cutoff power-law model, and so, b was fixed to 1. All remaining spatial and spectral parameters of sources in the initial ROI model were kept fixed at the 3FGL values. Once the spatial characteristics of the SNR had been approximated, as a second step, a global fit by `gtlike` was performed using the best-fit spatial model of the first step to determine a slightly more accurate estimate of the spectral parameters previously fitted by `pointlike`. Then, as a final step, to measure the energy spectrum of CTB 37A, the γ-ray data were fitted to the model in narrow energy bins. All parameters except the normalization of the bright sources within 1°3 from the ROI center (including CTB 37A, CTB 37B, RX J1713.7–3946, FL8Y J1714.8–3850, 3FGL J1718.1–3825, and 3FGL J1718.0–3726) were fixed to those previously found in the global fit. Doing so helps with avoiding numerical instabilities resulting from the fine binning in energy. During all the analyses, the normalizations of both the Galactic and isotropic diffuse emissions were left free to account for uncertainties in the diffuse emission.

For the following analyses, three primary sources of systematic errors have been taken into account: (1) uncertainties due to imperfect modeling of the Galactic diffuse emission, $\sigma_{IEM}$; (2) uncertainties on the source spatial model, $\sigma_{model}$; and (3) uncertainties on the PSF, $\sigma_{PSF}$. Systematic errors in the Galactic diffuse emission model were evaluated by going over the whole process using an alternative diffuse emission model (the Sample Model used in Ackermann et al. 2017). The systematic errors associated with spatial modeling of the source were obtained by comparing the SNR's estimated properties for the disk and Gaussian models provided by `pointlike`. Finally, for the P8R2_SOURCE_V6 event class

in our analysis, systematics on the PSF 68% containment radius were estimated to be <5% between 100 MeV and 10 GeV, increasing to 25% at 1 TeV as explained in the caveat page of the FSSC,<fnref n="7" /> which results in a systematic uncertainty of 0.005 on the 68% containment radius following Ackermann et al. (2018). Further details on the analyses and results are given in the subsequent subsections.

## 3.1. Morphological Analysis

The spatial analysis of the source was performed using all events above 1 GeV, taking advantage of a narrower instrumental PSF in conjunction with the significantly reduced contamination by the Galactic diffuse emission. The position and possible spatial extension of CTB 37A were examined using `pointlike` assuming a LogParabola spectral shape, $dN/dE = N_0(E/E_0)^{-(\alpha + \beta \ln(E/E_0))}$, as described in the 3FGL. Starting with the point-source hypothesis as our baseline model, we individually fitted the position and spectrum of CTB 37A and the two closest sources, CTB 37B and FL8Y J1714.8–3850. The iterative fitting process was progressively continued until changes in the estimated parameters from two consecutive fits became sufficiently small, and the derivative of the log-likelihood was close to zero (here, delta log-likelihood <0.1). To examine the possible angular extension of the SNR, we then replaced the point source on CTB 37A with an extended source (using a uniform disk and a 2D symmetric Gaussian model) and repeated the procedure until the desired fit with stable estimates was achieved. The significance of the source extension is quantified as the test statistics $TS_{ext} = 2 \log(\mathcal{L}_{ext}/\mathcal{L}_{ps})$, which compares the likelihood of an extended source hypothesis with a point-source one. For one additional degree of freedom, the extension significance is $\sqrt{TS}$ in units of $\sigma$.

After accounting for the systematic uncertainties associated with the localization due to the Galactic diffuse model and the LAT PSF (as $\sqrt{\sigma_{IEM}^2 + \sigma_{PSF}^2}$), we then implemented two empirical corrections by multiplying the total uncertainty by a correction factor of 1.05 (as in 3FGL) and adding the absolute 95% error of 0.0075 (as in 3FHL) in quadrature. The spatial properties of CTB 37A for the three tested models, along with their statistical and systematic errors, are summarized in Table 1, in which the systematic effects dominate the total uncertainties. Table 2 lists the best-fit position and preliminary spectral properties of CTB 37A and the two closest sources, CTB 37B and FL8Y J1714.8–3850. Our analysis confirms a 2D Gaussian morphology with a 68% containment radius of 0.116 and an extension significance of $5.75\sigma$. The obtained source size is smaller than the previously reported value by Li et al. (2017) by ~40%. This discrepancy can be explained by the contamination by the newly detected nearby source FL8Y J1714.8–3850 in our sky model. Our preliminary analysis based on the same data set but omitting FL8Y J1714.8–3850 also found a larger extension (Abdollahi et al. 2017). Figure 1 depicts the far and near vicinity of CTB 37A. The best-fit Gaussian extension (of 0.116 in radius) is coincident with the bright part of the radio shell toward the east, and it is larger than the TeV emission of Gaussian width 0.067 by ~75%. Additionally, both X-ray and TeV sources are offset to the west from the geometric center of the remnant (bottom left panel of Figure 1). The spatial properties of the remnant cannot

<footnote n="6">6 https://fermi.gsfc.nasa.gov/ssc/data/access/lat/fl8y/</footnote>

<footnote n="7">7 https://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_caveats.html</footnote>
<table><caption>**Table 1**
Best-fit Spatial Properties of CTB 37A for Different Morphological Models Above 1 GeV</caption><thead><tr><th rowspan="2">Spatial Model</th><th>R.A.</th><th>Decl.</th><th>r<sub>68</sub></th><th>TS</th><th>TS<sub>ext</sub></th><th>N<sub>dof</sub></th></tr><tr><th>(deg) ± stat ± sys</th><th>(deg) ± stat ± sys</th><th>(deg) ± stat ± sys</th><th></th><th></th><th></th></tr></thead><tbody><tr><td>Point Source</td><td>258.619 ± 0.006 ± 0.008</td><td>-38.515 ± 0.005 ± 0.008</td><td>...</td><td>1116 (1343)</td><td>...</td><td>5</td></tr><tr><td>Disk</td><td>258.642 ± 0.008 ± 0.019</td><td>-38.529 ± 0.008 ± 0.018</td><td>0.142 ± 0.025 ± 0.022</td><td>1511 (1902)</td><td>19.61 (19.91)</td><td>6</td></tr><tr><td>Gaussian</td><td>258.625 ± 0.007 ± 0.017</td><td>-38.513 ± 0.008 ± 0.011</td><td>0.116 ± 0.014 ± 0.017</td><td>1298 (1640)</td><td>33.10 (42.48)</td><td>6</td></tr></tbody></table>

**Note.** The best-fit positions for three tested models are given in J2000 epoch. The 68% containment radius for the disk and Gaussian models is defined as $r_{68} = 0.82r$ and $r_{68} = 1.51r$, respectively, where *r* is the disk radius or the Gaussian σ. The first and second errors on the spatial parameters correspond to the statistical and systematic, respectively. The test statistic is evaluated from the likelihood ratio between two models with and without the source of interest, TS = 2 log($L_{on}/L_{off}$). A comparison of source extension for different hypotheses is provided by TS<sub>ext</sub>. TS and TS<sub>ext</sub> values for the analysis using the alternative diffuse model are given in parentheses. N<sub>dof</sub> corresponds to the number of degrees of freedom for each model.

preclude an association with a PWN. Further discussions in this regard will be given in Section 4.

## 3.2. Spectral Analysis

After a global fit over the full energy range (0.2–200 GeV) using the best-fit model derived from the morphological analysis, the ROI model was reasonably optimized for generating the γ-ray spectrum of the SNR. A binned maximum likelihood analysis was performed in the full energy range, divided into nine logarithmically spaced energy bins by combining the four P8R2_SOURCE_V6 PSF event types in a joint likelihood approach. These four event types are based on the quality of the reconstructed direction, from the worst, PSF0, to the best, PSF3. Among different tested data sets, the data sample with a 68% PSF containment radius better than 1°7 (hereafter optimized data set) strikes a fair balance between minimizing the Galactic diffuse emission contribution and maintaining sufficient photon statistics. The optimized data set corresponds to PSF0 events with $E \ge 1$ GeV, PSF1 events with $E \ge 500$ MeV, PSF2 events with $E \ge 316$ MeV, and PSF3 events with $E \ge 200$ MeV. The top right panel of Figure 1 depicts a zoomed view of the ROI surrounding the SNR for this optimized data set. A power law with the spectral index fixed at 2.0 is used to measure the spectral points, which makes the results independent of the spectral model in each energy bin. As described above, only the normalizations of CTB 37A and the nearby bright sources with $>9\sigma$ significance (CTB 37B, RX J1713.7–3946, FL8Y J1714.8–3850, 3FGL J1718.1–3825, and 3FGL J1718.0–3726), as well as those of the Galactic and isotropic diffuse components, are adjusted. Fixing the nearby sources results in underestimating the uncertainties due to the diffuse emission below 500 MeV. The correction for energy dispersion is enabled for all model components except the Galactic and isotropic diffuse models, whose spectra are data-based.

Figure 2 shows the spectral energy distribution (SED) of CTB 37A. The statistical upper limit is calculated at 95% confidence level using a Bayesian method (see, e.g., Helene 1983) when the detection is not significant (TS < 4). Total systematic errors on the SED are included by adding in quadrature the uncertainties due to the Galactic diffuse model and the source spatial model as $\sqrt{\sigma_{IEM}^2 + \sigma_{model}^2}$. The dominant uncertainties on the SED are the systematic errors due to the Galactic diffuse emission below 1 GeV and the statistical errors due to smaller statistics above 20 GeV, which together can explain the mismatch between two best fits by pointlike and gtlike (see Figure 2). The total γ-ray energy flux of CTB 37A in the full energy range is calculated to be

($1.08 \pm 0.15_{\text{stat}} \pm 0.19_{\text{sys}}) \times 10^{-10} \text{ erg cm}^{-2} \text{ s}^{-1}$. The derived spectrum from the global fit with a spectral index of $1.92 \pm 0.19$ and a curvature of $0.18 \pm 0.05$ becomes softer above 10 GeV than the HESS spectrum with a spectral index of $2.3 \pm 0.13$ (Aharonian et al. 2008), which suggests two or more particle populations at the origin of the GeV and TeV emission. The quoted errors on the spectral index are statistical only.

## 4. Discussion

In the following discussion, based on the multiwavelength morphological and spectral characteristics of the remnant, we argue that CTB 37A should be classified as a composite SNR. We then review our model, which is constructed to examine the reacceleration and compression of preexisting CRs in the adjacent molecular clouds as the primary source of the observed emission toward the SNR, following an analytical approach given by Uchiyama et al. (2010). The naima package (Zabalza 2015) is used to perform the multiwavelength spectral fitting. We note that in all estimates presented below, a mean distance of 7.9 kpc to the SNR (Tian & Leahy 2012) is adopted.

### 4.1. CTB 37A: A New Composite SNR

The LAT observation of the GeV γ-rays toward CTB 37A reveals a spatial extension approximately two times larger than that reported in the TeV range by Aharonian et al. (2008). Furthermore, the LAT spectrum steepens above 10 GeV, so it is difficult to extend it smoothly into the harder TeV spectrum of HESS J1714–385. Taking these two facts together, it is very likely that CTB 37A is a composite system confining a nonthermal nebula inside its shell. Also, the VHE γ-ray source has an integrated γ-ray luminosity of $L_{\gamma(1-10 \text{ TeV})} = 1.73 \times 10^{34} \text{ erg s}^{-1}$ in the energy range from 1 to 10 TeV for the mean distance of 7.9 kpc. It is in good agreement with the representative value $L_{\gamma(1-10 \text{ TeV})} \approx 7 \times 10^{34} \text{ erg s}^{-1}$ for the observed middle-aged PWNe (between 7 and 23 kyr) as reported by HESS Collaboration et al. (2018).

Using XMM-Newton and Chandra observations of CXOU J171419.8–383023 of the northwestern rim of the SNR, Aharonian et al. (2008) derived a spectral index of $1.32_{-0.35}^{+0.39}$ for an extraction radius of 50". Yamauchi et al. (2014) analyzed the Suzaku data and obtained a softer spectral index of $1.94_{-0.14}^{+0.15}$ for an extraction region of radius ~2', quite typical index value for X-ray PWNe (see, e.g., Kargaltsev & Pavlov 2008). The apparent discrepancy in the derived X-ray spectral index could be the result of the radiative losses of very energetic electrons
<table><caption>**Table 2**
Spatial and Spectral Properties of CTB 37A and Two Nearby Sources Using the Gaussian Fit</caption><thead><tr><th>Source Name</th><th>R.A.<br>(deg) ± stat</th><th>Decl.<br>(deg) ± stat</th><th>TS</th><th>Spectral Index (α)</th><th>Curvature (β)</th></tr></thead><tbody><tr><td>CTB 37A</td><td>258.625 ± 0.007</td><td>-38.513 ± 0.008</td><td>1298</td><td>2.086 ± 0.135</td><td>0.103 ± 0.043</td></tr><tr><td>FL8Y J1714.8-3850</td><td>258.777 ± 0.026</td><td>-38.761 ± 0.018</td><td>128</td><td>2.661 ± 0.214</td><td>0.351 ± 0.245</td></tr><tr><td>CTB 37B</td><td>258.533 ± 0.017</td><td>-38.185 ± 0.014</td><td>82</td><td>1.611 ± 0.127</td><td>...</td></tr></tbody></table>

**Note.** The last two columns give the spectral index $\alpha$ at the reference energy $E_0 = 1$ GeV and the curvature $\beta$ for a LogParabola spectrum. The errors are only statistical.

with distance from the core of the emission. Not only fresh high-energy electrons but also accumulated aged cooled electrons can contribute to the Suzaku spectrum. This spectral softening away from the core detected in CXOU J171419.8-383023 is additional evidence to support a composite classification for the SNR and suggests a PWN origin for the TeV emission (see, e.g., Slane et al. 2000; Funk et al. 2007). Converting the X-ray flux to luminosity in the 0.5–10 keV band from Suzaku data yields $L_{X(0.5–10 keV)} = 8.22 \times 10^{34} \text{ erg s}^{-1}$, which is compatible with the suggested $L_{X(0.5–10 keV)} \sim 10^{35} \text{ ergs}^{-1}$ for the observed PWNe population (Gaensler & Slane 2006).

In the PWN scenario, the larger size of the VHE emission compared with the X-ray one can be explained by the different lifetime of the emitting electrons. In this picture, high-energy electrons producing the synchrotron X-ray emission suffer severe radiative losses (so-called fast cooling). In contrast, low-energy electrons can survive longer and produce the TeV emission via IC scattering of ambient photons (slow cooling). The observed offset between the peak of TeV and X-ray emission ($\delta\theta \sim 3.67'$) could arise from an asymmetric reverse shock resulting from the expansion of the SNR into an initially inhomogeneous medium so that the reverse shock typically returns faster from the direction of higher ambient density (e.g., Blondin et al. 2001).

The offset of the nonthermal X-ray emission from the center of the radio SNR ($\delta\theta \sim 2.0'$) can be explained by the proper motion of the pulsar with a transverse velocity of $\approx 750 d_{7.9} \tau_6^{-1} \text{ km s}^{-1}$, as supported by the tail-shape feature of the nonthermal X-ray emission (Aharonian et al. 2008). Here $d_{7.9} = d_{\text{SNR}}/(7.9 \text{ kpc})$ and $\tau_6 = \tau_{\text{SNR}}/(6 \text{ kyr})$ are the normalized distance and age of the remnant, respectively. The assumed age of 6 kyr for the remnant is explained in Section 4.2. The measured pulsar's proper motion is faster than the sound speed in the shocked medium ($C_s = 270 \text{ km s}^{-1}$), assuming an ideal gas with a mean temperature of $kT_X = 0.64 \text{ keV}$ for the whole SNR from Pannuti et al. (2014).

All evidence (size, offset, index, and luminosity) is reasonably consistent with a composite system, suggesting that the TeV and hard X-rays arise from the putative associated PWN, whereas the observed radio, soft X-rays, and GeV $\gamma$-ray emission originates in the SNR shock. The consistency between the center of the SNR shell and the GeV peak provides additional justification for the proposed scenario.

## 4.2. Modeling the Multiwavelength Emission from SNR
CTB 37A

Multiwavelength modeling of the remnant spectrum is conducted to probe the nature of the observed $\gamma$-ray emission toward the SNR. This approach includes the radio data at frequencies above 330 MHz taken from several surveys

(Kassim et al. 1991) and the LAT GeV spectrum from this work.

To explain the observed emission toward the system, we develop a numerical model in which acceleration of CRs by a strong shock propagating into the interstellar medium (ISM), as well as reacceleration of preexisting CRs in the nearby clouds by relatively weaker shocks, are considered. The model consists of three main parts: (i) production of energetic particles (including both protons and electrons), (ii) temporal evolution of the particle momentum distribution, and (iii) radiative mechanisms through the use of *naima*. In the following two subsections, we provide the required physical parameters and discuss the details of our model.

### 4.2.1. Physical Parameters of the Model

To specify the required physical parameters for the model, we follow the approach described in Devin et al. (2018) and the references therein. Here, we assume that the remnant is in the Sedov (adiabatic) phase of evolution. The post-shock temperature is related to the temperature of the X-ray thermal gas as $kT_s = kT_X/1.27 = 0.50 \text{ keV}$ (Borkowski et al. 2001). Moreover, the velocity of the blast wave assuming an adiabatic index of $\gamma = 5/3$ is given by:

$$v_s = \left( \frac{16 kT_s}{3 \mu} \right)^{1/2} \approx 650 \text{ km s}^{-1}, \quad (1)$$

where $\mu = 0.609 m_H$ is the mean atomic weight for a fully ionized gas. For the blast wave with a radius $R_s = 10 \text{ pc}$, the age of the remnant can then be estimated using $t_s = 2R_s/5v_s \approx 6000 \text{ yr}$. Consequently, the kinetic energy released through the SN explosion in terms of the shock dynamics is estimated by:

$$\frac{E_{\text{SN}}}{10^{51} \text{ erg}} = \left( \frac{R_s}{8.248 \text{ pc}} \right)^3 \left( \frac{n_0}{1 \text{ cm}^{-3}} \right) \left( \frac{kT_s}{1 \text{ keV}} \right) = 0.57, \quad (2)$$

where $n_0 = 0.64 \text{ cm}^{-3}$ is the ambient density. All the input parameters (radius, density, and gas temperature) are taken from Pannuti et al. (2014), assuming a distance of 7.9 kpc to the remnant as mentioned before. Equation (2) is derived from the standard Sedov solution $R_s^5 = \xi E_{\text{SN}} t_s^2 / \rho_0$, where $\xi = 2.026$ for the adiabatic index $\gamma = 5/3$, and $\rho_0 = 1.4 m_H n_0$ is the ambient mass density. The uncertainty on the distance to the SNR translates directly into the uncertainties in the remnant age and explosion energy.

The infrared observation of the SNR by Andersen et al. (2011) has revealed interactions with the surrounding clouds, in which radiative shocks occur, through the detection of [O I] emission at 63 $\mu$m, molecular and ionic lines. The ionic lines
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (Top left): Fermi-LAT counts maps in 0.2–200 GeV around the SNR CTB 37A in the entire ROI and using all data with a pixel size of 0.05 smoothed by a Gaussian kernel of 0.05. (Top right): Same as the top left panel, but for the optimized data set (68% PSF containment <1.7, as defined in the text) and a zoomed view of the inner 5° × 5° region centered on the SNR position. (Bottom left): Test Statistic (TS) map in the vicinity of CTB 37A in the 1–200 GeV energy range. The 843 MHz radio contours from MOST are overlaid in green. The position of CTB 37A, HESS J1714–385, and their statistical errors, are marked with blue and red vertical crosses, respectively. The 68% containment radii obtained by the Gaussian spatial model for CTB 37A and the HESS source (Aharonian et al. 2008) are shown with thick dashed blue and red circles, respectively. The inner and outer radii in each case represent the statistical errors on the fitted extension. The black diamond indicates the position of the X-ray source CXOU J171419.8–383023. Yellow cross is the position of FL8Y J1714.8–3850. (Bottom right): *Mopra* $^{12}$CO (J = 1–0) emission contours in green within the velocity range of −70 to −50 km s$^{-1}$ are overlaid on the TS map (same as the bottom left panel). The black crosses correspond to the position of the detected OH (1720 MHz) masers (Frail et al. 1996).</figcaption></figure>

suggest moderately fast J-type shocks with velocities of ~75–100 km s⁻¹ associated with compressed 10³.⁵–10⁴ cm⁻³ media in the northern and southern shell of the remnant. Therefore, it is reasonable to consider the contribution of both the main shock expanding into the interstellar medium and the radiative shock driven into the clouds for the γ-ray production in our model.

A dynamic pressure equilibrium between the main shock and the radiative shock is expressed as:

$$n_{0,cl} = k^2(v_s/v_{s,cl})^2 n_0. \qquad (3)$$

Thus for the main shock moving at velocity $v_s = 650$ km s⁻¹ into an ambient medium of density $n_0 = 0.64$ cm⁻³ and a radiative shock with velocity $v_{s,cl} = 100$ km s⁻¹ (Andersen et al. 2011), Equation (3) yields an upstream cloud density of $n_{0,cl} = 46$ cm⁻³, assuming a numerical factor $k = 1.3$ as in Uchiyama et al. (2010).

The downstream magnetic field strength in the cooled regions can be determined from

$$\frac{B_m^2}{8\pi} = k^2 \rho_0 v_s^2, \qquad (4)$$

balancing the magnetic pressure in the cooled regions with the ram pressure of the swept-up material in the shock (Hollenbach & McKee 1989), in which

$$B_m = \sqrt{2/3} (n_m / n_{0,cl}) B_{0,cl}. \qquad (5)$$

Thus, Equation (4) yields an amplified magnetic field strength of 520 µG in cooled regions of density $n_m = 10^{3.5}$ cm⁻³ (Andersen et al. 2011). Such a strong magnetic field is consistent with the measurements of the OH (1720 MHz) maser features spatially coincident with the associated clouds in CTB 37A (Brogan et al. 2000). The $\sqrt{2/3}$ coefficient in Equation (5) refers to the tangential component of the magnetic field, which is amplified by compression and prevents the collapse of the radiative shock by providing pressure support.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Spectral energy distribution of CTB 37A is measured by Fermi-LAT. Statistical and systematic uncertainties are given in black and red, respectively. An upper limit is calculated at 95% confidence level using a Bayesian method in which the TS is less than 4. The systematic effects on the upper limit are taken from the larger of the errors caused by the Galactic diffuse emission and source spatial model. The blue curve shows the best fit with pointlike above 1 GeV, extrapolated down to 0.2 GeV for better comparison. The purple curve indicates the global fit using gtlike over the full energy range (0.2–200 GeV). The shaded region in green represents a power-law fit of the HESS J1714–385 measurements, taking into account the statistical errors only.</figcaption></figure>

Using Equation (5), the upstream magnetic field strength in a cloud of density $n_{0,cl}$ is estimated to be $B_{0,cl} = 9.17$ µG.

The upstream magnetic field strength and density in the clouds are related by

$$B_{0,cl} = b(n_{0,cl}/\text{cm}^{-3})^{1/2} \mu\text{G}, \qquad (6)$$

where $b = V_A/(1.84 \text{ km s}^{-1})$ is a dimensionless parameter related to the Alfvén velocity $V_A$. As discussed in Hollenbach & McKee (1989), assuming that the velocity dispersion in molecular clouds is equal to the mean upstream Alfvén velocity, $b$ ranges between about 0.3 and 3. Here, the value of $b$ is constrained to ~1.4 by Equation (6), which is consistent with the expectations. Further, the downstream magnetic field (before radiative compression) can be obtained by $B_{d,cl} = \sigma_B B_{0,cl}$, where $\sigma_B = \sqrt{(2\sigma^2 + 1)/3}$ is the magnetic compression ratio in terms of the shock compression ratio $\sigma$ (Berezhko et al. 2002). It is assumed here that the magnetic field is fully isotropized. Then, the downstream magnetic field for a weakly modified shock ($\sigma = 4$) is found to be $B_{d,cl} = 30$ µG. Similarly, in the interstellar medium, the upstream magnetic field $B_{ISM} = 1.1$ µG and the downstream magnetic field $B_{d,ISM} = 3.6$ µG are derived for the ambient density of 0.64 cm$^{-3}$.

Once the physical characteristics of the main shock, the radiative shock, and the media they are moving into are set, the spectral properties of emission from both components can be inferred following the same procedure as that adopted in Parizot et al. (2006). In this approach, the acceleration timescale for $\sigma = 4$ is written as:

$$\tau_{\text{acc}} = 30.6 \times k_0(E) \times E_{\text{TeV}} B_2^{-1} v_{s,3}^{-2} \text{ yr}, \qquad (7)$$

which strongly depends on the shock velocity, and the synchrotron cooling timescale is estimated by

$$\tau_{\text{syn}} = (1.25 \times 10^3) \times E_{\text{TeV}}^{-1} B_2^{-2} \text{ yr}, \qquad (8)$$

where $B_2$ and $v_{s,3}$ are the magnetic field and the shock velocity in units of 100 µG and 1000 km s$^{-1}$, respectively. In Equation (7), $k_0$ is the deviation of the diffusion coefficient from the Bohm value$^8$ ($D/D_{Bohm}$), and it is found to be in the range from $\ge 1$ to 10. As the modeling of two young SNRs, Cas A (Zirakashvili et al. 2014) and RX J1713.7–3946 (Zirakashvili & Aharonian 2010), shows, the particle acceleration in young SNRs proceeds close to the Bohm diffusion limit for which $k_0$ is about 1. However, in the case of evolved SNRs, the particle acceleration is not efficient, which means the $k_0$ parameter can be larger. In our model, $k_0 = 10$ is used for the highest-energy electrons confined in the acceleration region. The maximum achievable energy of particles in the main shock is calculated by matching the acceleration time either to the synchrotron cooling timescale or to the shock’s age as $\tau_{\text{acc}} = \min\{\tau_{\text{syn}}, t_s\}$. In the case of the radiative shock, as in Uchiyama et al. (2010), the shock’s age is the time elapsed since the clouds were shocked. We approximate it by $t_c = t_s/2$, which leads to $\tau_{\text{acc}} = \min\{\tau_{\text{syn}}, t_c\}$. For the main shock moving into the interstellar medium of downstream magnetic field $B_{d,ISM} = 3.6$ µG, we obtain a maximum energy $E_{\text{max},e,p} \approx 300$ GeV. The slower radiative shock in the crushed clouds with downstream magnetic field of $B_{d,cl} = 30$ µG results in a lower value of $E_{\text{max},e,p} \approx 30$ GeV. Neither is limited by synchrotron cooling, so it is the same for electrons and protons.

Following Uchiyama et al. (2010), we account for a break in the spectrum of accelerated particles in the radiative shocks, where the gas is partially neutral. This feature is associated with Alfvénic turbulence and is due to particles accelerated along the upstream magnetic field in the clouds, which escape and result in a spectral steepening (Malkov et al. 2011).

The additional adiabatic compression as $s = (n_m/n_{0,cl})/\sigma$ in the radiative shocks boosts the momentum of accelerated particles by a factor $s^{1/3}$ close to 2.6. Table 3 summarizes the physical parameters in the model for both the main shock and the radiative shock in the clouds. $E_{\text{max}}$ and $\mathcal{E}_{\text{b,Alf}}$ are reported just downstream (before adiabatic compression).

### 4.2.2. Origin of the Gamma-Ray Emission: Accelerated or Reaccelerated CRs?

In the radiative shocks, the temporal evolution of the particle energy distribution is obtained by solving the transport equation, given by:

$$\frac{\partial N_i(E, t)}{\partial t} = -\frac{\partial}{\partial E} [\dot{E} N_i(E, t)] + Q_i, \qquad (9)$$

where $N_i(E, t)$ is the density of particles of species i (protons and electrons) per unit of energy, $\dot{E} \equiv dE/dt$ accounts for the total energy losses, and $Q_i$ denotes the injection rate of particles. It is assumed in the model that particles are confined to the expanding shell. Secondary $e^{\pm}$ generated through nuclear spallation are included. Contrary to the primaries, the injection rate of secondaries is time-dependent. To solve the transport equation numerically, two quantities on the right-hand

<footnote>$^8$ The Bohm value is calculated when the particle mean free path is equal to its Larmor radius $r_L = E/ZeB$, i.e., $D_{\text{Bohm}} = r_L c/3$.</footnote>
<table><caption>**Table 3**
Parameters for the Model Shown in Figure 3</caption><thead><tr><th>Emission Region</th><th>n (cm<sup>-3</sup>)</th><th>B (μG)</th><th>E<sub>max</sub> (GeV)</th><th>E<sub>b,Alf</sub> (GeV)</th><th>f</th><th>W<sub>p</sub> (×10<sup>48</sup> erg)</th><th>W<sub>e</sub> (×10<sup>47</sup> erg)</th></tr></thead><tbody><tr><td>Main Shock</td><td>0.64</td><td>3.6</td><td>300</td><td>no value</td><td>0.3</td><td>0.6 (15.8)</td><td>0.3 (7.9)</td></tr><tr><td>Radiative Shock</td><td>46</td><td>30</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Cooled Regions</td><td>3160</td><td>520</td><td>30</td><td>11<sup>†</sup></td><td>0.7<sup>*</sup></td><td>1.9</td><td>1.4</td></tr></tbody></table>

**Note.** Model parameters for both the main shock ($v_s = 650$ km s$^{-1}$) and radiative shocks ($v_{s,cl} = 100$ km s$^{-1}$) in SNR CTB 37A. The ambient medium properties (density and magnetic field) are taken from the observational data. The cutoff energy $E_{max}$ and the break energy $E_{b,Alf}$ are derived as in Parizot et al. (2006) and Malkov et al. (2011), respectively. $W_p$ and $W_e$ are the energy input to reaccelerated protons and (primary) electrons, respectively. Values in parentheses are the energy budget of CRs through the acceleration of fresh particles in the blast wave. A parameter marked with an asterisk is fit to the data. The spectral break marked with a dagger is caused by the Alfvén wave evanescence.

side of Equation (9), i.e., the energy losses and particle injection spectra, are specified as follows.

The energy loss mechanisms for electrons in the largely neutral clouds include synchrotron, bremsstrahlung, and ionization processes (Prantzos et al. 2011, and references therein). In general, ionization is the dominant loss process at low energies, and the contribution of bremsstrahlung and synchrotron to energy loss becomes important in the intermediate and high energy ranges, respectively. Protons in the radiative shocks suffer mainly from inelastic p-p interactions (Aharonian & Atoyan 1996) and ionization loss.

We used the spectra of the Galactic CR protons $n_{GCR,p}(p)$ and electrons+positrons $n_{GCR,e}(p)$ from Uchiyama et al. (2010) as the initial seed particles for considering the reacceleration process by the shocks in the model. The spectra of particles have been extended downward in energy to 20 MeV. It is much below the threshold energy of pion production, but it is necessary for computing their contribution to the total energy budget of the system. We assume that the CR properties in the clouds are the same as in the ISM.

The accelerated/reaccelerated CR spectrum $n_{acc}(p)$ (see Equation (3) in Uchiyama et al. 2010) through the DSA process (Blandford & Eichler 1987) undergoes several modifications. An exponential cutoff in momentum as a factor of $\exp[-(p/p_{max})]$ is introduced to take into account the maximum attainable energy derived using Equations (7) and (8). The DSA proton spectrum steepens by one power above the break momentum $p_{b,Alf}$ due to the evanescence of Alfvén waves in the radiative shocks. Moreover, the adiabatic compression inside cooling regions behind the radiative shock results in an enhanced energy density as $n_{ad}(p) = s^{2/3} n_{acc}(s^{-1/3} p)$.

Using the particle density in the radiative shocks, the injection rate of particles into the transport process is determined by:

$$Q_{p,e}^{r}(p) = \frac{f V_{\text{SNR}} n_{0,cl}}{n_{m} t_{c}} n_{\text{ad}}^{r}(p), \quad (10)$$

where the quantity $f$ is the filling factor of the clouds before they were crushed relative to the SNR volume $V_{\text{SNR}}$. $t_c$ represents the time spent by energetic particles in the shocked clouds. As in Uchiyama et al. (2010), it is assumed that the shock velocity remains constant inside the clouds.

In addition to the primary particles, secondary $e^{\pm}$ are also produced in the clouds of density $n_m$ through hadronic interactions and, consequently, production and decay of charged pions. Solving Equation (9) numerically, the proton spectrum in the crushed clouds $N_p^r(p, t)$ serves as the source term for calculating the spectrum of secondary $e^{\pm}$. The

injection rate of secondaries is computed as

$$Q_{e^{\pm}}^{r}(E, t) = n_{m} c \epsilon_{M} \int dT_{p} N_{p}^{r}(T_{p}, t) \frac{d\sigma(E, T_{p})}{dE}, \quad (11)$$

where $c$ is the speed of light, and $\epsilon_M$ is the nuclear enhancement factor to take into account the effects of heavier nuclei ($A > 1$) in both CRs and the target material (Mori 2009). The last term $d\sigma(E, T_p)/dE$ is the energy-dependent differential cross section of $e^{\pm}$ generation from incident protons of kinetic energy $T_p$ as parametrized by Kamae et al. (2006). Contrary to the primaries, the injection rate of secondary $e^{\pm}$ is time-dependent. A low-energy cutoff at a kinetic energy of 1 MeV is assumed for secondaries.

Using all the required ingredients, the numerical solution of the transport equation provides the final spectrum of accelerated CR species $N(p, t)$. Equation (9) is solved at $t_c$ for both preexisting particles in the clouds and the secondaries. We then apply *naima* for computing nonthermal radiation processes and fitting the multiwavelength data toward the remnant. The radiation from neutral pion-decay, bremsstrahlung, synchrotron, and IC are included. The interstellar radiation field (ISRF) at the location of the remnant is taken from GALPROP (Porter et al. 2017) as seed photons for IC. Energy density and temperature of the ISRF in the near-infrared band, far-infrared band, and the cosmic microwave background (CMB) photons are ($U_{\text{NIR}} = 1.84$ eV cm$^{-3}$, $T_{\text{NIR}} = 3484.16$ K), ($U_{\text{FIR}} = 0.81$ eV cm$^{-3}$, $T_{\text{FIR}} = 28.98$ K), ($U_{\text{CMB}} = 0.26$ eV cm$^{-3}$, $T_{\text{CMB}} = 2.72$ K), respectively. The $\gamma$-rays from p-p interactions are rescaled by the nuclear enhancement factor $\epsilon_M$.

Additionally, we considered the contribution of acceleration of fresh CRs at the blast wave. At this site, the energy losses of particles are negligible, and there are no secondaries, so we use the analytic approximation of the solution of the transport equation. In our model, the energy distribution of electrons is approximated by a power law with a break associated with the aging of the particles and an exponential cutoff as:

$$N_e(E) = A_e \exp(-E/E_{\text{max}}) \times \begin{cases} (E/E_0)^{-\Gamma_{e,1}} & E \le \mathcal{E}_{\text{b,syn}} \\ (\mathcal{E}_{\text{b,syn}}/E_0)^{\Gamma_{e,2}-\Gamma_{e,1}} (E/E_0)^{-\Gamma_{e,2}} & E > \mathcal{E}_{\text{b,syn}} \end{cases} \quad (12)$$

in which the electron spectral index changes by unity ($\Gamma_{e,2} - \Gamma_{e,1} = 1$) after the break energy $\mathcal{E}_{b,syn}$. The cutoff energy $E_{max}$ corresponds to the maximum achievable particle energy within the acceleration process, and the reference energy $E_0$ is set to 1 GeV. For the proton spectrum, a power law with spectral index $\Gamma_p$ and an exponential cutoff at $E_{max}$ in
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Acceleration/reacceleration scenario for the multiwavelength modeling of the emission toward CTB 37A, in which dashed lines correspond to the emissions from primaries (p) and dotted lines correspond to those of secondaries (s) in the radiative shocks. Dash-dotted lines represent the emissions from fresh CRs accelerated at the blast wave. The total model from the summation of all emissions is shown in black. The IC emission from the reaccelerated electrons in the clouds is insignificant, with a peak at ~2 × 10⁻¹⁵ erg cm⁻² s⁻¹. The radio data points are a combination of several observations at frequencies above 330 MHz, taken from Table 2 of Kassim et al. (1991). The Suzaku X-ray and HESS TeV spectra are from Yamauchi et al. (2014) and Aharonian et al. (2008), respectively. The error bars on the LAT measurement are calculated by adding statistical and systematic errors in quadrature.</figcaption></figure>

the form

$$N_p(E) = A_p (E/E_0)^{-\Gamma_p} \exp(-E/E_{\text{max}}) \quad (13)$$

is assumed. The proton and electron spectral indices $\Gamma_p = \Gamma_{e,1} = 2$ are set as predicted by the DSA mechanism. Equations (12) and (13) are in terms of energy, as defined in naima. The $\gamma$-rays from p-p interactions are enhanced by the factor $\epsilon_M$. $\mathcal{E}_{b,syn}$ is set by equating the synchrotron cooling timescale with the shock's age through $\tau_{syn} = t_s$ in the main shock. Due to the weak downstream magnetic field at the main shock, the synchrotron cooling is insignificant, and the resulting spectral break is beyond the cutoff energy ($\mathcal{E}_{b,syn} \propto B_{d,ISM}^{-2}$). So, the electron distribution is a single power law as Equation (13). The total kinetic energy of the accelerated protons is set to $10\%(1-f)E_{SN}$, i.e., $W_p = 1.58 \times 10^{49}$ erg. The energy budget of accelerated electrons is set to 5% that of the protons.

Because all physical parameters (except the energy budget of CRs through the acceleration process) are constrained by observations, the level of gamma-ray emission at the radiative shocks is defined only by the filling factor $f$ of the clouds and is fit to the data. The resulting factor of 0.72 is large but not impossible. The fit requires 0.36% of $E_{SN}$ to be transferred to the radiative shocks. In this model, the maximum energy $E_{max}$ with the assumed $k_0 = 10$ provides a reasonable fit to the data. The resulting spectrum using these values implies that the radiative shells driven into the dense clouds are the dominant contributor to the observed radio and GeV emission (see dashed-line curves in Figure 3), and the emission from the blast wave is negligible, similar to SNR W44 as tested by our model. Secondary particles are found to contribute insignificantly to the total radio and GeV emission. Moreover, the model predicts that 13% of the energy injected into the protons was lost since the clouds were shocked. Because the physical parameters are entirely constrained in our model, it is impossible to fit the ratio of the radio to gamma-ray data exactly. The model ratio on Figure 3 is about twice too low.

<caption>**Table 4**</caption>

Physical Properties of Three Middle-aged SNRs Explained with the Crushed Clouds Scenario

<table><thead><tr><th></th><th>CTB 37A</th><th>MSH<br>15-56</th><th>W44</th></tr></thead><tbody><tr><td colspan="4"><strong>Physical Parameters</strong></td></tr><tr><td colspan="4"><strong>SNR Dynamics</strong></td></tr><tr><td>τ (kyr)</td><td>6</td><td>16.5</td><td>10</td></tr><tr><td>d (kpc)</td><td>7.9</td><td>4.1</td><td>2.9</td></tr><tr><td>R (pc)</td><td>10</td><td>21</td><td>12.5</td></tr><tr><td>E<sub>SN</sub> (×10<sup>51</sup> erg)</td><td>0.57</td><td>0.5</td><td>5</td></tr><tr><td colspan="4"><strong>Gas Properties in Clouds</strong></td></tr><tr><td>v<sub>s,cl</sub> (km s<sup>-1</sup>)</td><td>100</td><td>150</td><td>100</td></tr><tr><td>n<sub>0,cl</sub> (cm<sup>-3</sup>)</td><td>46 (183)</td><td>2 (8)</td><td>200 (800)</td></tr><tr><td>B<sub>0,cl</sub> (μG)</td><td>9 (30)</td><td>4 (14)</td><td>25 (83)</td></tr><tr><td>n<sub>m</sub> (cm<sup>-3</sup>)</td><td>3160</td><td>90</td><td>10620</td></tr><tr><td>B<sub>m</sub> (μG)</td><td>520</td><td>160</td><td>1080</td></tr><tr><td colspan="4"><strong>Gas Properties in Blast Wave</strong></td></tr><tr><td>v<sub>s</sub> (km s<sup>-1</sup>)</td><td>650</td><td>500</td><td>490</td></tr><tr><td>n<sub>0</sub> (cm<sup>-3</sup>)</td><td>0.64<br>(2.56)</td><td>0.1 (0.4)</td><td>5 (20)</td></tr><tr><td>B<sub>ISM</sub> (μG)</td><td>1.1 (3.6)</td><td>3.0 (10.0)</td><td>4.0 (13.1)</td></tr><tr><td colspan="4"><strong>Spectral Parameters in Radiative<br>Shocks</strong></td></tr><tr><td>p<sub>c</sub> (GeV/c)</td><td>30</td><td>82.7</td><td>122</td></tr><tr><td>p<sub>b</sub> (GeV/c)</td><td>11.3</td><td>15.2</td><td>7</td></tr><tr><td>α<sub>r</sub></td><td>0.50</td><td>0.34</td><td>0.37</td></tr><tr><td><sup>a</sup>References</td><td>(1)</td><td>(2)</td><td>(3)</td></tr></tbody></table>

**Note.** $\tau$ is the age, and $d$ is the distance to SNR. Upstream density and magnetic field in the clouds are given by $n_0$, $c_l$ and $B_0$, $c_l$, respectively. These parameters in the ISM are presented by $n_0$ and $B_{ISM}$, respectively. The numbers in the parentheses correspond to those in the downstream regions. The two parameters $n_m$ and $B_m$ correspond to those in the cooled regions. $\alpha_r$ refers to the radio spectral index. The spectral break $p_b$ is due to Alfvén wave damping in CTB 37A and W44, and to synchrotron cooling in MSH 15–56. All other parameters are the same as those explained in the text.

<sup>a</sup> The numbers refer to the following references: (1) This work, (2) Devin et al. (2018), (3) Uchiyama et al. (2010).

Having a larger magnetic field in the radiative shocks (i.e., a larger pressure in the SNR) would improve that. We note that for the results presented here, the distance to the remnant is the main source of uncertainty in our calculations.

## 5. Comparison with Other Middle-aged SNRs

As reported in Table 5 of Acero et al. (2016), among all surveyed Galactic SNRs by Fermi-LAT, only 11 of them have associated OH masers, which are signposts of SNR/MC interaction. The measured magnetic field strength of ~mG toward those masers is consistent with those for shock-compressed clouds (see, e.g., Brogan et al. 2000). Table 4 lists three interacting SNRs along with their physical parameters, for which the reacceleration of particles, followed by compression in radiative shocks ($v_{s,cl} \ge 100$ km s⁻¹), is believed to be at the origin of hadronic $\gamma$-ray emission. The compressed gas and amplified magnetic field by the fast shocks in the interaction region ($v_{s,cl} \gtrsim 100$ km s⁻¹) play a crucial role in the hadronic nature of their $\gamma$-ray emission.

As shown in Table 4, the maximum achievable energy in CTB 37A is lower than the two other SNRs. It is mainly due to
a shorter acceleration time compared with the two other SNRs, a lower downstream magnetic field in its associated clouds compared with that in W44, and a slower shock in the clouds relative to MSH 15–56 ($E_{max} \propto t_c v_{s,cl}^2 B_{d,cl}$).

In CTB 37A, the emission from secondaries does not exceed the primaries' spectra, contrary to the results in Uchiyama et al. (2010) for W44. Interestingly, as in Lee et al. (2015), our model results in a subdominant contribution of secondaries in SNR W44 using the same physical parameters as listed in Table 4. Lee et al. (2015) have attributed the observed difference to the fact that a time-dependent shock velocity is assumed in their model. However, we have followed the same approach as Uchiyama et al. (2010), in which the dynamics of the remnant are not taken into account. The difference in the flux of secondary electrons may stem from how the nonthermal cooling in the radiative shock is computed in Uchiyama et al. (2010).

## 6. Conclusions

Using 8 yr of Pass 8 Fermi-LAT data, we have studied the nature of the γ-ray emission in the direction of the CTB 37A system. The morphological analysis of the source using all data above 1 GeV revealed an extended emission that is best modeled by a Gaussian distribution with $r_{68} = 0.116 \pm 0.014_{\text{stat}} \pm 0.017_{\text{sys}}$ at a significance of 5.75σ. The measured angular extension is comparable with the radio size, while it is larger than the extension of the TeV emission from HESS J1714–385 by ~75%. The GeV emission is centered on the radio SNR and offset from both the nonthermal X-ray and the TeV emission. The spectral analysis of the remnant using an optimized data set with $r_{68}$ less than 1.7 over the full energy range (0.2–200 GeV) showed that the GeV spectrum steepens above 10 GeV compared with the HESS spectrum, which strengthens two different origins for the GeV and TeV emission toward this system. The detected pulsar inside the system confirms a composite class of the SNR, as we proposed.

Assuming the SNR is in the Sedov stage, we examined a scenario in which both the acceleration of Galactic CRs and the reacceleration of preexisting CRs are considered. The maximum energy of particles through the DSA process reaches 300 GeV in the blast wave and 30 GeV in the radiative shocks. The low-density ambient medium leads to faint emission from the blast wave. On the contrary, the dense clouds interacting with the SNR can explain well both the radio and GeV spectra through the reacceleration of CRs followed by radiative compression. For all physical parameters fixed at their values from the observational data, a reasonable fit to the radio and GeV spectra can be obtained if the clouds occupied most of the SNR volume prior to explosion. The energy left in the CR protons implies that 13% of the injected energy was lost since the clouds were shocked. Moreover, the contribution of secondaries is subdominant compared with the primaries.

The accumulated Fermi-LAT data will provide more candidates in which the reacceleration of CRs is the dominant process. The Cherenkov Telescope Array will resolve the CTB 37A system and clarify the relation between the SNR and the PWN.

The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the

National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K.A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden.

Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515.

Facility: Fermi-LAT.

Software: Fermi Science Tools<fnref n="9" /> (version v11r05p02), naima package (Zabalza 2015).

## ORCID iDs

Soheila Abdollahi https://orcid.org/0000-0002-6803-3605

Jean Ballet https://orcid.org/0000-0002-8784-2977

Yasushi Fukazawa https://orcid.org/0000-0002-0921-8837

Hideaki Katagiri https://orcid.org/0000-0003-2347-8819

Benjamin Condon https://orcid.org/0000-0002-9850-5108

## References

* Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009a, ApH, 32, 193

* Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009b, ApJS, 183, 46

* Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 718, 348

* Abdollahi, S., Mizuno, T., Fukazawa, Y., et al. 2017, ICRC (Busan), 301, 743

* Acero, F., Ackermann, M., Ajello, M., et al. 2015, ApJS, 218, 23

* Acero, F., Ackermann, M., Ajello, M., et al. 2016, ApJS, 224, 8

* Ackermann, M., Ajello, M., Albert, A., et al. 2017, ApJ, 840, 43

* Ackermann, M., Ajello, M., Allafort, A., et al. 2013, Sci, 339, 807

* Ackermann, M., Ajello, M., Baldini, L., et al. 2018, ApJ, 237, 32

* Aharonian, F. A., & Atoyan, A. M. 1996, A&A, 309, 917

* Aharonian, F. A., Akhperjanian, A. G., Barres de Almeida, U., et al. 2008, A&A, 490, 685

* Ajello, M., Atwood, W. B., Baldini, L., et al. 2017, ApJS, 232, 18

* Andersen, M., Rho, J., Reach, W. T., Hewitt, J. W., & Bernard, J. P. 2011, ApJ, 742, 7

* Atwood, W., Albert, A., Baldini, L., et al. 2013, arXiv:1303.3514

* Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071

* Axford, W. I., Leer, E., & Skadron, G. 1977, ICRC (Budapest), 11, 132

* Bell, A. R. 1978a, MNRAS, 182, 147

* Bell, A. R. 1978b, MNRAS, 182, 443

* Bell, A. R. 1987, MNRAS, 225, 615

* Berezhko, E. G., Ksenofontov, L. T., & Völk, H. J. 2002, A&A, 395, 943

* Blandford, R. D., & Eichler, D. 1987, PhR, 154, 1

* Blandford, R. D., & Ostriker, J. P. 1978, ApJL, 221, L29

* Blondin, J. M., Chevalier, R. A., & Frierson, D. M. 2001, ApJ, 563, 806

* Borkowski, K. J., Lyerly, W. J., & Reynolds, S. P. 2001, ApJ, 548, 820

* Braiding, C., Wong, G. F., Maxted, N. I., et al. 2018, PASA, 35, 29

* Brandt, T. J. & Fermi-LAT Collaboration 2013, AdSpR, 51, 247

* Brogan, C. L., Frail, D. A., Goss, W. M., & Troland, T. H. 2000, ApJ, 537, 875

* Castro, D., & Slane, P. 2010, ApJ, 717, 372

* Castro, D., Slane, P., Carlton, A., & Figueroa-Feliciano, E. 2013, ApJ, 774, 36

* Dawnes, A. 1984, MNRAS, 210, 845

* Devin, J., Acero, F., Ballet, J., & Schmid, J. 2018, A&A, 617, A5

* Fermi, E. 1949, PhRv, 75, 1169

* Frail, D. A., Goss, W. M., Reynoso, E. M., et al. 1996, AJ, 111, 1651

* Funk, S., Hinton, J. A., Moriguchi, Y., et al. 2007, A&A, 470, 249

* Gaensler, B. M., & Slane, P. O. 2006, ARA&A, 44, 17

<footnote n="9">⁹ https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/</footnote>
Ginzburg, V. L., & Syrovatskii, S. I. 1964, The Origin of Cosmic Rays
(Oxford: Elsevier)

Helene, O. 1983, NIMPR, 212, 319

HESS Collaboration, Abdalla, H., Abramowski, A., et al. 2018, A&A, 612, A2

Hollenbach, D. J., & McKee, C. F. 1989, ApJ, 342, 306

Jogler, T., & Funk, S. 2016, ApJ, 816, 100

Joubert, T., Castro, D., Slane, P., & Gelfand, J. 2016, ApJ, 816, 63

Kamae, T., Karlsson, N., Mizuno, T., et al. 2006, ApJ, 647, 692

Kargaltsev, O., & Pavlov, G. G. 2008, in AIP Conf. Proc. 983, 40 Years of Pulsars: Millisecond Pulsars, Magnetars, and More, ed. C. Bassa et al. (Melville, NY: AIP), 171

Kassim, N. E., Baum, S. A., & Weiler, K. W. 1991, ApJ, 374, 212

Kerr, M. 2010, PhD thesis, Univ. Washington, Seattle, arXiv:1101.6072v1

Krymskii, G. F. 1977, DoSSR, 234, 1306

Lande, J., Ackermann, M., Allafort, A., et al. 2012, ApJ, 756, 5

Lee, Sh., Patnaude, D. J., Raymond, J. C., et al. 2015, ApJ, 806, 71

Li, J., Rea, N., Torres, D. F., & de Oña-Wilhelmi, E. 2017, ApJ, 835, 30

Malkov, M. A., Diamond, P. H., & Sagdeev, R. Z. 2011, NatCo, 2, 194

Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396

Maxted, N. I., Rowell, G. P., Dawson, B. R., et al. 2013, MNRAS, 434, 2188

Milne, D. K., & Dickel, J. R. 1975, AuJPh, 28, 209

Milne, D. K., & Hill, E. R. 1969, AuJPh, 22, 211

Mori, M. 2009, APh, 31, 341

Pannuti, T. G., Rho, J., Heinke, C. O., & Moffitt, W. P. 2014, AJ, 147, 55

Parizot, E., Marcowith, A., Ballet, J., & Gallant, Y. A. 2006, A&A, 453, 387

Porter, T. A., Jóhannesson, G., & Moskalenko, A. V. 2017, ApJ, 846, 67

Prantzos, N., Boehm, C., Bykov, A. M., et al. 2011, RvMP, 83, 1001

Reynoso, E. M., & Mangum, J. G. 2000, ApJ, 545, 874

Saz Parkinson, P. M., Limyansky, B., Clark, C. J., et al. 2018, in 8th Int. Fermi Symp., https://fermi.gsfc.nasa.gov/science/mtgs/symposia/2018/program/115/

Sezer, A., Gök, F., Hudaverdi, M., & Ercan, E. N. 2011, MNRAS, 417, 1387

Slane, P., Yang, C., Schulz, N. S., et al. 2000, ApJL, 533, L29

Tian, W. W., & Leahy, D. A. 2012, MNRAS, 421, 2593

Uchiyama, Y., Blandford, R. D., Funk, S., Tajima, H., & Takaaki, T. 2010, ApJL, 723, L122

Whiteoak, J. B. Z., & Green, A. J. 1996, A&AS, 118, 329

Yamauchi, S., Minami, S., Ota, N., & Koyama, K. 2014, PASJ, 66, 2

Zabalza, V. 2015, ICRC (The Hague), 34, 922

Zirakashvili, V. N., & Aharonian, F. A. 2010, ApJ, 708, 965

Zirakashvili, V. N., Aharonian, F. A., Yang, R., Oña-Wilhelmi, E., & Tuffs, R. J. 2014, ApJ, 785, 130