An Investigation of the Interstellar Environment of Supernova Remnant CTB 87

We present a new millimeter CO-line observation towards supernova remnant (SNR) CTB 87, which was regarded purely as a pulsar wind nebula (PWN), and an optical investigation of a coincident surrounding superbubble. The CO observation shows that the SNR delineated by the radio emission is projectively covered by a molecular cloud (MC) complex at V$_{\rm {LSR}}$ = $-60$ to $-54$ km s$^{-1}$. Both the symmetric axis of the radio emission and the trailing X-ray PWN appear projectively to be along a gap between two molecular gas patches at $-58$ to $-57$ km s$^{-1}$. Asymmetric broad profiles of $^{12}$CO lines peaked at $-58$ km s$^{-1}$ are found at the eastern and southwestern edges of the radio emission. This represents a kinematic signature consistent with an SNR-MC interaction. We also find that a superbubble, $\sim 37'$ in radius, appears to surround the SNR from HI 21cm (V$_{\rm {LSR}} \sim -61$ to $-68$ km s$^{-1}$), WISE mid-IR, and optical extinction data. We build a multi-band photometric stellar sample of stars within the superbubble region and find 82 OB star candidates. The likely peak distance in the stars' distribution seems consistent with the distance previously suggested for CTB 87. We suggest the arc-like radio emission is mainly the relic of the part of blastwave that propagates into the MC complex and is now in a radiative stage while the other part of blastwave has been expanding into the low-density region in the superbubble. This scenario naturally explains the lack of the X-ray emission related to the ejecta and blastwave. The SNR-MC interaction also favours a hadronic contribution to the {\gamma}-ray emission from the CTB 87 region.


Introduction
The progenitors of core-collapse (CC) supernovae (SNe) are massive stars, which are expected to be born mostly in OB associations, and their environments are altered by energy feedback processes dominated by energetic stellar winds and SN explosions as well as strong ionizing radiation. Such processes can produce superbubbles filled with low-density material and significantly affect the evolution of supernova remnants (SNRs) therein. On the other hand, remnants of CC SNe, are often close (a few pc) to molecular clouds (MCs), the birthplace of the progenitor stars that end their evolution in short lifetimes. So far about 70 SNRs are confirmed or suggested to be associated with MCs (Jiang et al. 2010;Chen et al. 2014b) among the known ∼ 300 SNRs in the Milky Way (Ferrand & Safi-Harb 2012;Green 2014). Six types of observational evidence for SNR-MC interaction are summarized, including the 1720 MHz OH maser, molecular line broadening, morphological agreement of molecular features with SNR features, etc. (Jiang et al. 2010;Chen et al. 2014b). In this paper, we investigate a perplexing SNR, CTB 87, that may interact with an adjacent MC but may be located near the inner edge of a superbubble.
CTB 87 has been classified as a filled-center type SNR, with a radio size of about 8 ′ × 6 ′ , centered at R.A.= 20 h 16 m 02 s , decl.= 37 • 12 ′ 00 ′′ . It was first cataloged in a radio survey at 960 MHz of the Galactic plane conducted by the Owens Valley Radio Observatory (Wilson & Bolton 1960). The distance to it was suggested to be 12 kpc according to HI absorption measurements (Green & Gull 1989;Wallace et al. 1997). Based on the extinction-distance relation (Foster & Routledge 2002), a new distance of 6.1 ± 0.9 kpc was established, which implies that CTB 87 is located in the Perseus spiral arm (Kothes et al. 2003). In X-rays, CTB 87 is centrally brightened with a size of about 5 ′ as observed by the Einstein satellite observatory (Wilson 1980). A detailed Chandra ACIS analysis shows that the pulsar wind nebula (PWN) harbours a point source, which may putatively be a pulsar, 100 ′′ southeast away from the radio emission peak (Matheson et al. 2013). In γ-rays, the TeV point source VER J2016+371, spatially coincident with CTB 87, is resolved by the VERI-TAS telescope system (Aliu et al. 2014). The GeV Fermi -LAT source 3FGL J2015.6+3709 is spatially close to VER J2016+371, with its origin still under debate (Kara et al. 2012;Acero et al. 2015Acero et al. , 2016Saha 2016). Recently, the multiwavelength spectrum obtained from the CTB 87 region has been interpreted to be comprised of contributions from a Maxwellian distribution of electrons and a broken power-law distribution of electrons (Saha 2016).
CTB 87 has been proposed to be associated with molecular gas at a systemic local standard of rest (LSR) velocity V LSR of about −58 km s −1 (Huang & Thaddeus 1986;Cho et al. 1994;Kothes et al. 2003), while it has also been argued to lie at the inner boundary of an expanding HI superbubble at a V LSR of about −70 km s −1 (Wallace et al. 1997). The molecular gas was resolved into several subclumps, with CTB 87 seemingly located between two of them (Cho et al. 1994). However, no kinematic evidence has yet been found for the SNR-MC interaction. Also, no 1720 MHz OH maser, a reliable signpost of interaction between SNR and MC, has been detected for CTB 87 (Frail et al. 1996). About 30% of the SNRs that are confirmed to be in physical contact with MCs do not show OH masers (Jiang et al. 2010), most probably because the shocked molecular gas does not satisfy the appropriate physical conditions, namely density of order ∼ 10 5 cm −3 and temperature in the range of 50-125 K (Lockett et al. 1999).
In previous work, Cho et al. (1994) mapped the CTB 87 region in CO lines with a half-power beamwidth of 2.7 ′ and a grid spacing of 3 ′ , and Kothes et al. (2003) presented a mapping in CO (J= 1-0) with a beamwidth of 2 ′ and a grid spacing of 1 ′ , as well as a partial mapping in CO (J= 3-2). In this paper, we present a new CO-line observation toward the CTB 87 region and an optical investigation of the coincident superbubble, aiming to explore the environment of CTB 87 and the interaction of the SNR with the adjacent MC. We first describe the CO observations and data reduction process in §2; we analyze the CO-line data in detail in §3 and discuss the main result in §4. Finally, we present our summary in §5.

Observations and data reduction
Our observations of millimeter emission toward SNR CTB 87 were made simultaneously in 12 CO (J= 1-0), 13 CO (J= 1-0), and C 18 O(J =1-0) lines during 2013 May 11-18 using the 13.7m millimeter-wavelength telescope of the Purple Mountain Observatory at Delingha. We used the new 3 × 3 pixel Superconducting Spectro-scopic Array Receiver as the front end, which is constructed with Superconductor-Insulator-Superconductor mixers using the sideband separating scheme (Zuo et al. 2011;Shan et al. 2012). We did on-the-fly mapping toward a 34 ′ ×60 ′ field that covers the full angular extent of SNR CTB 87 with a grid spacing of 30 ′′ . The half-power beamwidth of the telescope is about 52 ′′ and the pointing accuracy of the telescope is better than 4 ′′ . The typical system temperature is about 220K for 115.2GHz and 130K for 110.2GHz. The bandwidth of the spectrometer is 1 GHz, with 16,384 channels and a spectral resolution of 61 KHz. Thus the velocity resolution is 0.16 km s −1 for 12 CO and 0.17 km s −1 for 13 CO and C 18 O. All the CO-line data were reduced using the GILDAS/CLASS package developed by the IRAM observatory 1 .
For a multiwavelength investigation of the environment of this SNR, we also used Chandra X-ray (ObsID: 11092, PI: Safi-Harb; which we reduced with CIAO ver. 4.7), and WISE 22 µm mid-infrared (IR) (WISE Science Data Center, IPAC, Caltech) data. The HI line emission data were obtained from the Canadian Galactic Plane Survey (CGPS; Taylor et al. 2003), and the 1.4 GHz radio continuum emission data were from the NRAO VLA Sky Survey (NVSS; Condon et al. 1998). Optical and near-IR data from the Panoramic Survey Telescope and Rapid Response System (Pan-Starrs; Kaiser et al. 2002), the Isaac Newton Telescope (INT) Photometric Hα Survey of the Northern Galactic Plane (IPHAS; Drew et al. 2005), and the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006) were used to investigate the OB star candidates within the coincident superbubble. Figure 1 shows the averaged CO spectra from a region covering SNR CTB 87. There are two prominent 12 CO (J= 1-0) emission peaks, at around −58 km s −1 and −38 km s −1 . The 13 CO (J= 1-0) emission is only prominent at ∼ −58 km s −1 , and the C 18 O (J= 1-0) emission is not significant across the velocity range.

Spatial distribution of the ambient clouds
We made 12 CO emission channel maps around the two peaks with velocity interval 0.5 km s −1 (see Figure 2 and Figure 3) to examine the spatial distribution of the two molecular components. The ∼ −38 km s −1 12 CO component seems to overlap the SNR at its southeastern "apex" by projection in interval −39 to −37.5 km s −1 , but there is no systematic morphological feature correspondence between the molecular gas and the SNR (Figure 2).
The ∼ −58 km s −1 component appears to spread over a large area covering the SNR, and a few LSR-velocity dependent structures are noteworthy. (1) In velocity interval −55.5 to −54.5 km s −1 , there is a linear structure of 12 CO (J= 1-0) emission along the eastern edge of the SNR (Figure 3), somewhat similar to that described in Kothes et al. (2003) (in passing, note that a Galactic coordinate system was used there while an equatorial coordinate system is used here, and the shell-like structure shown there is not complete due to the limited field of view). The structure can also be discerned in the 13 CO (J= 1-0) channel map at −55 km s −1 (Figure 4). (2) In velocity interval −56.5 to −55.5 km s −1 , a bar-like molecular structure appears to pass through the SNR along the northeast-southwest orientation and, in particular, through the brightness peak of the radio emission. (3) In interval −58 to −57 km s −1 , both 12 CO and 13 CO emissions show two patches of molecular material in the eastern and southwestern edges of the SNR, similar to the sub-clumps described in Cho et al. (1994). The integrated 12 CO (J= 1-0) intensity map in the main velocity range of the ∼ −58 km s −1 component, −60 to −54 km s −1 , is presented in Figure 5 (in green). It is overlaid with Chandra 0.5-7 keV X-ray image of CTB 87 (in blue) and the radio continuum image (in red). The X-ray emitting part of the PWN, with a southeast-northwest oriented elongation, appears to be located in a gap of molecular gas, where CO emission is weak, between the two patches at −58 to −57 km s −1 .

Kinematic evidence for SNR-MC interaction ?
We have inspected the 12 CO (J= 1-0) line profiles of the two molecular components toward CTB 87. For the ∼ −38 km s −1 component, we do not find any asymmetric broad profiles of the 12 CO line in the grid of CO spectra. (This component, as mentioned in §3.1, has neither morphological feature corresponding to the SNR). For the other component, around −58 km s −1 , we show a grid of 12 CO (J= 1-0) and 13 CO (J= 1-0) spectra in the velocity range of −61 to −53 km s −1 ( Figure 6) and find asymmetric broad profiles of the 12 CO line in the eastern and southwestern edges. Figure 7 shows the average line profiles of some of the pixels at the edges of the SNR (regions "E" and "SW" marked in Figure 6). There is a secondary peak at −56 to −55 km s −1 in region "E" (also seen in most of the pixels therein) besides the main peak at −58 km s −1 , with little 13 CO counterpart ( 3σ). Although the broad red (right) wing may include a contribution from real line broadening due to shock disturbation, we cannot exclude contamination by line-of-sight emission in a wide region. We note that there is similar secondary peak on the east of the radio boundary and there is non-negligible intensity at −56 to −55 km s −1 on the northeast of the boundary. However, in the three bottom pixels of region "E", there are unique plateaus in the blue (left) wings (∼ −61 to −59 km s −1 ), which are also clearly (> 4σ) reflected in the average profile for region "E" (top panel, Figure 7). This plateau-like spectral feature is somewhat similar to that detected in the 12 CO spectra in the western edge of SNR 3C397 (Jiang et al. 2010). In the spectra of the pixels in the southwestern edge, as typified by region "SW", the line profiles are characterized by strong peaks at −58 km s −1 for both 12 CO and 13 CO and a broad red wing of 12 CO extending to ∼ −54 km s −1 without a significant 13 CO counterpart (bottom panel, Figure 7). The surrounding pixels essentially display different shapes of line profiles, and the emission outside the southern edge becomes weak. Since 13 CO emission, which is usually optically thin (like that with τ ( 13 CO)≪ 1 given in Tables 1 and 2), is yielded in quiescent, intrinsically high-column-density molecular gas, a broad 12 CO-line wing without a 13 CO counterpart is very likely to represent a perturbed gas deviating from the systemic LSR velocity. Therefore, the asymmetric 12 CO wings from the edges of CTB 87 (especially the blue wings in region "E" and the red wings in region "SW") very likely result from Doppler broadening of the −58 km s −1 line and thus might provide a kinematic evidence for interaction between the ∼ −58 km s −1 MC complex and the SNR.
We fit the CO emissions with Gaussian lines for the ∼ −58 km s −1 molecular gas in regions "R" and "L" (as defined in Figure 5), and the derived parameters, molecular column density N(H 2 ), excitation temperature T ex , and optical depth of 13 CO (J= 1-0) τ ( 13 CO) are summarized in Table 1 and Table 2. The distance to the MC/SNR is taken to be 6.1 kpc, as suggested by Kothes et al. (2003). The column density of H 2 and the mass of the molecular gas are estimated using two methods. In the first method, the conversion relation for the molecular column densities, N(H 2 ) ≈ 7 × 10 5 N ( 13 CO) (Frerking et al. 1982), is used under the assumption of local thermodynamic equilibrium for the molecular gas and optically thick condition for the 12 CO (J= 1-0) line. In the second method, a value of the CO-to-H 2 mass conversion factor, N(H 2 )/W ( 12 CO) (known as the "X-factor"), 1.8 × 10 20 cm −2 K −1 km −1 s (Dame et al. 2001) is adopted.
3.3. The superbubble at V LSR ∼ −64 km s −1 By projection, SNR CTB 87 is located within a superbubble, centered at R.A.=20 h 14 m 03 s , decl.=37 • 13 ′ 27 ′′ , with a circular boundary (brightened in the south and north) of angular radius ∼ 37 ′ , and nears its eastern edge (shown in Figure 8). We investigate the WISE mid-IR observation toward SNR CTB 87 and find the superbubble is bright at both 12 and 22µm. The HI emission in this sky region also shows a cavity at the LSR velocity around −64 km s −1 (−61 to −68 km s −1 ), which happens to be spatially coincident with the mid-IR superbubble and to have the same angular size. Furthermore, such a bubble-like structure has a counterpart in the optical extinction map (also see Figure 8).
As superbubbles are usually the products of energy and material feedback from massive OB stars into the interstellar medium, we searched for possible OB stars in the direction toward this superbubble.
We built a multi-band photometric stellar sample of over 0.18 million stars within the projected region of the superbubble (i.e., the circle shown in Figure 8) by cross-matching the photometric cataloges of the Pan-Starrs, IPHAS, and 2MASS. We used a spectral energy distribution (SED) fitting algorithm, similar to those employed by Chen et al. (2014a) and Mohr-Smith et al. (2015), to obtain the effective temperature T eff , optical extinction A V , and distance modulus µ of the individual stars. We compared the observed SED of each star to the stellar models from the Padova isochrone data base (CMD v3.0;Bressan et al. 2012;Marigo et al. 2017). We considered only the main-sequence (surface gravity log g > 3.8 dex), and Galactic thin disk metallicity ([Fe/H] = −0.2 dex) models. The observed magnitudes of a star can be modeled by m obs where m obs i and m mod i are the observed and Padova magnitudes in the ith band (i=1, 2, ..., 11 corresponding to Pan-Starrs g P , r P , i P , z P , y P , IPHAS r ′ , i ′ , Hα, and 2MASS J, H, and K S bands, respectively), A i is the extinction in the ith band, and µ is the distance modulus. We adopt the R V = 3.1 extinction law from Cardelli et al. (1989) to convert the optical extinction A V to extinction in each individual band A i . A simple Bayesian scheme based on the Markov Chain Monte Carlo sampling is adopted to obtain the parameters log T eff , A V , and µ for the individual stars. We adopt the likelihood where σ i is the uncertainty of the ith band observed magnitude, and the priors We adopt the results for the stars with high photometric precisions (σ i ≤ 0.1 mag for each band) and high posterior probabilities (log P ≥ 5, where 5 is around the median value of the logarithmic posterior probabilities). We note that, due to the lack of information in blue optical wavelengths (i.e., u band), there could be degeneracies between the effective temperatures and the amount of extinction for the individual stars.
As a result, we found 82 OB star candidates. Figure 9 shows the distribution of distance and optical extinction for all these objects. There seems to be a peak of the distribution of distance of these OB candidates at d=7.3 kpc, with a dispersion of 1.9 kpc, although the selected stars may be an incomplete and contaminated sample of the candidates owing to the lack of the u band data. Notably, this distance range covers that of the Perseus spiral arm in the direction of CTB 87 (Kothes et al. 2003, Figure 8 therein). As some massive stars should be responsible for the superbubble, it is most likely that they are located in the dispersion range of the peak. The distance range of the massive stars or the superbubble also seems to be consistent with the distance to SNR CTB 87, either 6.1 ± 0.9 kpc, that is derived from the intervening neutral hydrogen column density (Kothes et al. 2003) or ∼ 6-8 kpc estimated from the column density of the line-of-sight X-ray absorbing hydrogen atoms and molecules (see § 4.1).

SNR Environment
Our new CO-line observation shows spatial correspondence of the SNR CTB 87 with an MC complex at LSR velocities around −58 km s −1 (−54 to −60 km s −1 ). We have also demonstrated a probable kinematic evidence of the SNR shocking against molecular patches at V LSR ∼ −58 km s −1 . Hereafter we parameterize the distance to CTB 87 and the associated MC as d = 6.1d 6.1 kpc. We can make a rough estimate of the density of the associated molecular gas. We adopt angular sizes of the two patches of molecular gas in regions "R" and "L" as ∼ 7.5 ′ × 7.5 ′ and ∼ 3.9 ′ × 3.5 ′ (i.e., ∼ 13.3d 6.1 pc × 13.3d 6.1 pc and ∼ 6.9d 6.1 pc × 6.2d 6.1 pc), respectively, and assume line-of-sight sizes 13.3 pc and 6.6 pc for them. The estimated gas density, n(H 2 ), and gas mass, M, for the two patches are given in Tables 1 and 2, respectively.
The multiwavelength spectrum from radio to γ-rays for CTB 87, which incorporates the TeV emission of VER J2016+371 (Aliu et al. 2014) and the GeV emission from 3FGL J2015.6+3709 (Kara et al. 2012;Acero et al. 2015Acero et al. , 2016, has been explained using a leptonic scenario by Saha (2016), with a combination of a broken power law and a Maxwellian distribution of electrons. Nonetheless, it is also mentioned that the observed GeV-TeV data can be explained by a hadronic scenario as long as there is a dense ambient target medium with a mean hydrogen atom density ∼ 20 cm −3 or higher (depending on detailed hadronic interaction mechanisms) (Saha 2016). Actually, it is notable that the centroid of the VER-ITAS very high energy γ-ray emission (Aliu et al. 2014) is essentially coincident with both the radio emission of the SNR and the ∼ −58 km s −1 MC complex. The possible SNR-MC interaction revealed here provides a likely hotbed for the production of hadronic γ-rays. The gas density of the two molecular patches, a few tens cm −3 , estimated here, appears to satisfy the hadronic scenario.
CTB 87 is very likely to be located within (but near the eastern boundary of) the revealed superbubble, in view of the following issues. (1) The LSR velocity of the 12 CO emission of the associated MC, −54 to −60 km s −1 , is very similar to that of the HI emission for the superbubble, −61 to −68 km s −1 .
(2) The distance to the SNR, derived both from the foreground intervening neutral hydrogen column density, 6.1 ± 0.9 kpc (Kothes et al. 2003), and from the X-ray-absorbing hydrogen column density, 6 to 8 kpc (see below), is consistent with the likely peak distance of the OB star candidates that are projected inside the superbubble, 7.3±1.9 kpc, which can be regarded as the distance to the superbubble and the OB star candidates inside. The hydrogen (including atoms and molecules) column density for the entire diffuse X-ray emission region of CTB 87 is N H = 1.21-1.57 × 10 22 cm −2 (Matheson et al. 2013). The optical extinction can then be estimated to be A V = 4.8 ± 0.7 from the empirical relation N H = (2.87 ± 0.12) × 10 21 A V cm −2 (Foight et al. 2016), which indicates a distance of ∼ 6-8 kpc to CTB 87 using the A V -distance relation toward this SNR given in Kothes et al. (2003). (3) There is little detection of thermal X-rays related to the SN ejecta and blast wave, which sets an upper limit 0.2 cm −3 to the density of the medium (Matheson et al. 2013). This can be readily explained if the SN exploded material is blown into the tenuous and hot gas within the superbubble (in CTB 87, another part of the gas may have shocked into the MC; see §4.2). In the superbubble, the SN ejecta can freely move away at a very high speed, and the blast shock will propagate with a small Mach number so that it can be expected to be weak and become thermalized, with most of the SN kinetic energy carried by this part of the material being deposited in the thermal energy of the hot medium (Tang & Wang 2005).
The MC with which the SNR probably interacts may survive the strong radiation of the OB stars in the superbubble. The photodissociation timescale of the MC can be comparable to, or even larger than, the age of the superbubble. For simplicity, the dissociating far ultraviolet (FUV) radiation of the OB stars (of number N OB ), to which the MC is exposed, is approximated as being from the bubble center, and the molecular column density along the bubble radial is approximated to the line-of-sight column density N(H 2 ). The distance from the MC to the center, R MC , is not less than the projected distance to the center (∼ 22.5 ′ or ∼ 40d 6.1 pc). Each dissociating photon is absorbed by a hydrogen molecule with a photodissociation probability, p D ∼ 0.15 (Draine 2011). Thus, the photodissociation timescale of the MC is given by τ MC N(H 2 )(4πR 2 MC )/(p D N OB S D ), where S D is the mean production rate of Lyman band dissociating photons of an OB star (∼ 10 47.5 s −1 , Diaz-Miller et al. 1998). It is estimated to be τ MC 6.4 × 10 6 (N(H 2 )/1 × 10 21 cm −2 )(N OB /20) −1 (S D /10 47.5 s −1 ) −1 d 2 6.1 yr, where the reference value 1 × 10 21 cm −2 is adopted for N(H 2 ) according to Table 1 and Table 2, and N OB ∼ 20 is adopted as a reference value considering this approximate number of stars are possibly inside and responsible for the superbubble (based on Figure 9). On the other hand, the age of the superbubble is estimated to be (Weaver et al. 1977) t SB ∼ 5.8 × 10 6 d 5/3 6.1 (N OB /20) −1/3 × (L w /10 34 erg s −1 ) −1/3 (n 0 /0.6cm −3 ) 1/3 yr, where L w is the power of the wind of an OB star, and n 0 is the atomic number density of the environment medium into which the superbubble expands. Here we use a density of 0.6 cm −3 typical for a warm interstellar HI gas (which takes the largest volume fraction next to the coronal gas in the Galactic disk) to be a reference value for n 0 (Draine 2011). Thus their ratio is τ MC /t SB 1.1 × (N(H 2 )/1 × 10 21 cm −2 )(N OB /20) −2/3 d 1/3 6.1 (S D /10 47.5 s −1 ) −1 (L w /10 34 erg s −1 ) 1/3 (n 0 /0.6cm −3 ) −1/3 . Here we have ignored the formation of new molecules and the extinction of the dissociating photons by dust grains in the molecular gas. This estimate indicates that, with some typical or proper values of parameters adopted, the photodissociation timescale appears not less than the age of the superbubble, and thus the survival of the MC in the superbubble could be reasonable.
It has been suggested that CTB 87 is inside another superbubble found in HI emission at V LSR ∼ −70 km s −1 , with a radius in the range ∼ 38 ′ -70 ′ , centered at approximately R.A.= 20 h 16 m 15 s , decl.=36 • 40 ′ (J2000) (Wallace et al. 1997). This superbubble, however, is not as favored as the ∼ −64 km s −1 one in view of the bigger offset of the LSR velocity ∼ −70 km s −1 from the CO-line velocity of the associated MC (∼ −58 km s −1 ) than that of ∼ −64 km s −1 . Incidentally, we do not find the counterpart of this HI superbubble in the WISE near-and mid-IR observations.

Radio and X-ray Configuration
The X-ray observation of CTB 87 shows a cometary-like trailing structure, ∼ 200 ′′ ×300 ′′ or ∼ 5.9d 6.1 pc × 8.9d 6.1 pc in size at 0.3-7 keV (Matheson et al. 2013). The radio emission takes a blow-out, arc-like shape, with the X-ray trail at its symmetric axis, and has a remarkably larger size (∼ 6 ′ × 8 ′ or 10.6d 6.1 pc × 14.2d 6.1 pc, even up to 16 ′ ) than the X-ray emission and a different brightness peak location from the X-ray one (Kothes et al. 2003;Matheson et al. 2013). The apparent one-sided confinement of the X-ray emission represents a typical structure for a relative oriented motion between the ambient gas and the PWN. But why is the radio emission much more extended than the X-ray nebula, with different locations for the brightness peaks in the two bands? Matheson et al. (2013) suggest that the radio emission of CTB 87 is a relic PWN that was crushed by the reverse shock propagating back from somewhere external, and the X-ray nebula is the new trailing PWN after the passage of the reverse shock and the subsequent oscillation of the nebula. The pulsar has moved for about 10 kyr southeastward from the explosion site, assumed to be at the location of the radio brightness peak. However, as the authors point out, no sign of a forward shock (especially in the nearby southeastern region) has been found in this scenario. Also, according to the simulation in Blondin et al. (2001), van der Swaluw et al. (2004), and Kolb et al. (2017), relic radio PWN are apparently swept/left behind the head of the new, small X-ray nebula. In particular, given the physical contact of the eastern and southwestern edges of the extended radio emission with the MC, this seems to leave no room for the reverse shock that shocks against and crushes the suggested radio relic PWN.
With the probable interaction of CTB 87 with the ambient MC revealed here, we suggest instead that the radio emission is mainly a remnant of the blastwave that propagates into the MC at V LSR ∼ −58 km s −1 , although it may include a contribution from the PWN. This scenario also allows for a reverse shock to have moved backward from the remnant's edge and crushed the PWN, forming the trailing morphology of X-ray emission, as in the case of other, similar PWNe (e.g. N157B, van der Swaluw et al. 2004;Chen et al. 2006). This scenario naturally addresses the question, as noted in Matheson et al. (2013), of the lack of the X-ray emission related to the ejecta and blastwave.
The progenitor exploded somewhere along the symmetric axis of the radio emission or the X-ray trail, with a part of the blastwave shocking against the MC and the other part expanding, and quickly becoming sufficiently faint in a low-density region, possibly a portion of the superbubble. The shock in the MC decelerated and entered the intense radiative stage shortly thereafter, with the shocked gas at that point at too low a temperature to emit X-rays. Actually, the cooling time scale of the shock propagating into the molecular gas is t cool = 2.8× 10 3 (E/10 51 erg) 0.24 (n a /80 cm −3 ) −0.52 yr (Falle 1981), where E is the SNR's explosion energy and n a is the preshock H atom density. For an average molecular density n(H 2 ) ∼ 40 cm −3 (see Tables 1 and 2), the cooling time is much shorter than the spin down time ∼ 1 × 10 4 yr (Matheson et al. 2013). Although the SNR may have a part blown out, for simplicity we crudely estimate the evolution of the radiative part in the MC as a complete sphere with the radius as r s ∼ 5(ǫ/0.24) 5/21 (E/10 51 erg) 5/21 (n a /80 cm −3 ) −5/21 (t/10 4 yr) 2/7 pc, where ǫ ∼ 0.24 is the energy fraction (McKee & Ostriker 1977;Blinnikov et al. 1982). If we approximate the pulsar's spindown time of ∼ 10 4 yr as the age of the remnant and adopt an explosion energy ∼ 10 51 erg, then the radius is ∼ 5 pc, very similar to the size of the radio emission. The shock velocity is v s ∼ (2/7)r s /t ∼ 140(ǫ/0.24) 5/21 (E/10 51 erg) 5/21 (n a /80 cm −3 ) −5/21 (t/10 4 yr) −5/7 km s −1 . Since MCs are usually highly clumpy, the observed shocked molecular gas showing broad COline wings should be in dense clumps, in which the shock velocity can be well below 50 km s −1 (Draine & McKee 1993) so that the molecules are not dissociated. The radio emission as a relic of the blastwave in this scenario indicates a composite nature for SNR CTB 87.

Summary
We have performed a new millimeter CO-line observation toward the region of CTB 87, which was thought to be a filled-center type SNR, and an optical investigation of the coincident superbubble. The CO observation shows that the SNR delineated by the radio emission is projectively covered by a bar-like molecular structure at V LSR = −56.5 to −55.5 km s −1 . Both the symmetric axis of the radio emission and the trailing X-ray PWN appear projectively to be at a gap between two molecular gas patches at −58 to −57 km s −1 . Asymmetric broad line profiles of the ∼ −58 km s −1 12 CO-line are obtained from the molecular gas at the eastern and southwestern boundary of the radio emission. This could well be a kinematic evidence of the physical contact between CTB 87 and the ∼ −58 km s −1 (−60 to −54 km s −1 ) MC complex. A superbubble, ∼ 37 ′ in angular radius and centered at (20 h 14 m 03 s , 37 • 13 ′ 27 ′′ , J2000), seemingly surrounding the SNR, is found in HI 21cm (V LSR ∼ −61 to −68 km s −1 ), WISE mid-IR, and optical extinction observations. We built a multi-band photometric stellar sample of over 0.18 million stars within the superbubble region and found 82 OB star candidates. The distribution of the stars' distances is likely peaked at 7.3 kpc (with a dispersion of 1.9 kpc) and seems consistent with the previously suggested distance of 6.1 ± 0.9 kpc for CTB 87. We suggest the arc-like radio emission is mainly a relic of the part of the blastwave that is driven into the MC complex, and is now in a radiative stage, while the other part of the blastwave has been expanding into the low-density region, very likely in the superbubble. This scenario naturally explains the lack of X-ray emission related to the ejecta and blastwave. The SNR-MC interaction also favors a hadronic contribution to the γ-ray emission from the CTB 87 region.
We are grateful to the staff of Qinghai Radio Observing Station at Delingha for their help during the observation. We highly appreciate Xin Zhou and Gao-Yuan Zhang for the advice on the data analysis, and Xiao Zhang for the discussion about PWN physics. This work is supported by the 973 Program grants 2015CB857100 and2017YFA0402600 andNSFC grants 11233001, 11633007, 11773014, 11503008, and11590781. This research has made use of the NVSS data and NASA's Astrophysics Data System 2 . Center ( km s −1 ) FWHM ( km s −1 ) T peak (K) W (K km s −1 ) 12 CO(J= 1-0) −57.5 2.6 3.6 9.9 13 CO(J= 1-0) −57.1 2.0 0.6 1.3
b See the text for the two estimating methods.
c The excitation temperature is calculated from the maximum 12 CO(J= 1-0) emission point of the −58 km s −1 component.   Table 1 Fig. 1.-Average CO spectra from a 6.5 ′ × 6.5 ′ region (centered at R.A.=20 h 16 m 06 s .7, decl.=37 • 12 ′ 32 ′′ ) covering the pulsar wind nebula, in the velocity range of −80 km s −1 -20 km s −1 . The black line is for 12 CO (J= 1-0), the green line for 13 CO (J= 1-0), and the red line for C 18 O(J= 1-0).  Black lines denote 12 CO spectra, red lines denote 13 CO, and dashed lines denote the 0 K main-beam temperature. The size of each pixel is 30 ′′ × 30 ′′ . The radio contours are the same as those in Figure 2. Two regions delineated by blue rectangles (labelled as "E" and "SW") are used to extract CO spectra, in which redward broadened wings are shown for the 12 CO emission at systemic velocity ∼ −58 km s −1 (see Figure 7).  Figure 6. The black line denotes the 12 CO spectra, the green line represents 13 CO, and the dotted line represents the 0 K main-beam temperature. Fig. 9.-Distribution of distance (top panel) and optical extinction (bottom panel) for the OB star candidates (log P ≥ 5 and log T eff ≥ 4.3) within the projected region of the superbubble (i.e., the circle shown in Figure 8). The red curve in the top panel is the Gaussian fitting result of the peak-like component.