ALMA CO Observations of Gamma-Ray Supernova Remnant N132D in the Large Magellanic Cloud: Possible Evidence for Shocked Molecular Clouds Illuminated by Cosmic-Ray Protons

Observations of ¹²CO(J = 3–2) line emission at λ = 0.87 mm wavelength were conducted in 2014 September 1–3 (PI: H. Sano, proposal No. AC141006) using the ASTE 10 m radio telescope (Ezawa et al. 2004). The telescope is installed at an altitude of 5000 m in the Atacama Desert in Chile, operated by the Chile Observatory of the National Astronomical Observatory of Japan (NAOJ). We used the on-the-fly mapping mode with Nyquist sampling, and the effective observation area was $3.6^{\prime} \times 3.6^{\prime}$ centered at (${\alpha }_{{\rm{J}}2000}$, ${\delta }_{{\rm{J}}2000}$) ∼ (${05}^{{\rm{h}}}{25}^{{\rm{m}}}2$8, $-69^\circ 38^{\prime} 35^{\prime\prime}$). The front end was a sideband-separating superconductor–insulator–superconductor (SIS) mixer receiver "CATS 345" (Inoue et al. 2008). We utilized an XF-type digital spectrometer "MAC" (Sorai et al. 2000) as the back end. The bandwidth of MAC is 128 MHz with 1024 channels, corresponding to a spectral resolution of 0.125 MHz. The velocity coverage and resolution are thus ∼111 km s⁻¹ and ∼0.11 km s⁻¹, respectively. The typical system temperature was ∼300 K, including the atmosphere in the single-side band. To derive the main beam efficiency, we observed N 159W (${\alpha }_{{\rm{J}}2000},{\delta }_{{\rm{J}}2000}={05}^{{\rm{h}}}{40}^{{\rm{m}}}3$7, $-68^\circ 47^{\prime} 00^{\prime\prime} ;$ Minamidani et al. 2011), obtaining a main beam efficiency of 0.67 ± 0.08. We also observed the M-type AGB star R Dor every hour to satisfy pointing offset accuracy within $2^{\prime\prime}$. After convolution with a two-dimensional Gaussian kernel, we obtained the cube data with a beam size of ∼$23^{\prime\prime}$ (∼5.6 pc at the LMC distance). The typical noise fluctuations are ∼0.046 K at the velocity resolution of 0.4 km s⁻¹.

Observations of ¹²CO(J = 1–0) line emission at λ = 2.6 mm wavelength were carried out using the ALMA Band 3 (86–116 GHz) in Cycle 2 as an early science project (PI: H. Sano, proposal No. 2013.1.01042.S). We utilized 40 antennas of the 12 m array, nine antennas of the 7 m array, and three antennas of the total power (TP) array. The effective observation area was a $150^{\prime\prime} \times 150^{\prime\prime}$ rectangular region centered at (${\alpha }_{{\rm{J}}2000}$, ${\delta }_{{\rm{J}}2000}$) ∼ (${05}^{{\rm{h}}}{25}^{{\rm{m}}}2$79, $-69^\circ 38^{\prime} 34\buildrel{\prime\prime}\over{.} 2$). The combined baseline length of the 12 m and 7 m arrays ranges from 7.2 to 215.7 m, corresponding to u–v distances from 2.8 to 82.9 $k\lambda$. The two quasars J0334−4008 and J0635−7516 were used for complex gain calibrators. Another two quasars J0601−7036 and J0526−6749 were observed as phase calibrators. We also observed Callisto, Uranus, and a quasar J0519-454 as flux calibrators. The data reduction was performed using the Common Astronomy Software Application (CASA; McMullin et al. 2007) package version 5.4.0. We used the multiscale CLEAN task implemented in the CASA package (Cornwell 2008). The scale parameters are $0^{\prime\prime}$, 318, and 954 for the 12 m array and $0^{\prime\prime}$, 1464, and 4392 for the 7 m array. To improve the imaging quality, we also applied a uvtaper during the clean procedure for the 12 m array data. The u–v tapering applies a multiplicative Gaussian taper to the spatial frequency space to downweight high spatial frequencies. This can suppress artifacts—for example, strong side lobes—arising from poorly sampled areas near and beyond the maximum spatial frequency. We finally combined the cleaned data of the 12 and 7 m array data sets and calibrated the TP array data by using the feather task. The final beam size of the feathered data is $5\buildrel{\prime\prime}\over{.} 26\times 4\buildrel{\prime\prime}\over{.} 99$, with a position angle of $59\buildrel{\circ}\over{.} 60$, corresponding to a spatial resolution of ∼1.2 pc at the LMC distance. The typical noise fluctuations of the feathered data are ∼0.22 K at a velocity resolution of 0.4 km s⁻¹.

To investigate the CO gas distribution at larger spatial scales, we used the Magellanic Mopra Assessment Data Release 1 (MAGMA DR1, Wong et al. 2011). MAGMA is a ¹²CO(J = 1–0) mapping survey of the LMC using the Mopra 22 m radio telescope of the Australia Telescope National Facility (ATNF). The angular resolution is 45'', corresponding to the spatial resolution of ∼11 pc at the LMC distance. The typical noise fluctuations of a region surrounding N132D are ∼0.26 K at the velocity resolution of 0.53 km s⁻¹. We applied additional spatial smoothing with a two-dimensional Gaussian kernel. The angular resolution of the smoothed data is ∼60'' (∼15 pc at the LMC distance), which is the same resolution as for the H i survey data of the LMC (see Section 3.2).

To better understand the distribution of neutral atomic hydrogen toward N132D, we used archival survey data of the H i line at λ = 21 cm wavelength published by Kim et al. (2003). The survey data were obtained using the Australia Telescope Compact Array (ATCA) and Parkes 64 m telescopes operated by ATNF. The angular resolution of the survey data is $60^{\prime\prime}$, corresponding to the spatial resolution of ∼15 pc at the LMC distance. The typical noise fluctuations of brightness temperature are ∼2.4 K at the velocity resolution of 1.689 km s⁻¹.

We used archival X-ray data obtained by Chandra, for which the observation IDs are 5532, 7259, and 7266 (PI: K. J. Borkowski, proposal No. 06500305), which have been published in previous papers (e.g., Borkowski et al. 2007; Xiao & Chen 2008; Schenck et al. 2016; Sharda et al. 2020). The data sets were taken with the Advanced CCD Imaging Spectrometer S-array (ACIS-S3). Table 1 lists the details of the observations. We utilized Chandra Interactive Analysis of Observations (CIAO, Fruscione et al. 2006) software version 4.12 with CALDB 4.9.1 (Graessle et al. 2007) for data reduction. All of the data sets were reprocessed using the chandra_repro task. We then created exposure-corrected, energy-filtered images using the fluximage task in the energy bands 0.35–0.85 keV, 0.5–1.2 keV (soft band), 0.85–1.6 keV, 1.2–2.0 keV (medium band), 1.6–6.0 keV, 2.0–7.0 keV (hard band), and 0.5–7.0 keV (broad band). The total effective exposure is ∼89.3 ks. For the spectral analysis, we used HEASOFT (version 6.25), including spectral fitting with XSPEC (version 12.10.1f, Arnaud 1996). We fit the spectra in the energy band 0.3–10.0 keV, and the errors of the fitted parameters are quoted at the 1σ confidence level unless specified otherwise. We fit the unbinned spectra to preserve the maximum spectral information. Following Sharda et al. (2020), we explicitly model the background as opposed to subtracting it. We use the C statistic (Cash 1979) as the minimization statistic to avoid the well-known bias introduced by the χ² statistic in the case of a low number of counts per spectral bin (Kaastra 2017), and we report the Pearson χ² (weighting by the model) to evaluate goodness of fit. We used ATOMDB version 3.0.9 (Foster et al. 2013) and nonequilibrium ionization (NEI) version 3.0.4 for the NEI models (Borkowski et al. 2001). We used the cosmic abundance given by Wilms et al. (2000) as the baseline abundance for all our analysis and the cross sections given by Verner et al. (1996).

Table 1.  Chandra ACIS-S Observation Log of SNR N132D

ObsID	Observation Date	Exposure	${\alpha }_{{\rm{J}}2000}$	${\delta }_{{\rm{J}}2000}$
		(ks)	(h m s)	(° ' '')
05532	2006 Jan 09	44.59	05 25 02.28	−69 38 37.32
07259	2006 Jan 10	24.85	05 25 02.28	−69 38 37.32
07266	2006 Jan 15	19.90	05 25 02.28	−69 38 37.32

Download table as:  ASCIITypeset image

To investigate the origin of the hard X-ray emission in N132D, we also used a map of the hard-band X-ray (E: 10–15 keV) obtained with NuSTAR (Bamba et al. 2018). The NuSTAR observations were executed on 2015 December 10–11. The total effective exposure is 62.3 ks. The angular resolution is $18^{\prime\prime}$ (half-power beam width, HPBW) or $58^{\prime\prime}$ (half-power diameter, HPD). To improve the signal-to-noise ratios of the map, we smoothed the data with a two-dimensional Gaussian kernel of 20''.

To compare the spatial distribution with the interstellar medium (ISM) environment of N132D, we also used an excess count map of teraelectronvolt gamma rays obtained by the High Energy Stereoscopic System (H.E.S.S. Collaboration et al. 2015). The angular resolution is ∼3' for the point-spread function (PSF, 68% containing radius) or ∼7' for the FWHM, corresponding to a spatial resolution of ∼44 pc for the PSF or ∼100 pc for the FWHM.

The optical data of the Hα and [O iii] emission lines are used to derive the spatial distributions of the ionized gas and shocked ejecta. We utilized the Hubble Space Telescope (HST) Wide Field Planetary Camera 2 (WFPC2) and the Advanced Camera for Survey (ACS) images from the Hubble Legacy Archive.¹⁷ The observations were carried out using the ACS Wide Field Channel F658N for Hα and the WFPC2 F502N for [O iii], which have been published by Morse et al. (1996) and Borkowski et al. (2007). For further details about the data reductions and pipeline processes, we refer the reader to the HST Data Handbook.¹⁸

Figure 1 shows a map of ATCA and Parkes H i intensity overlaid with the Mopra ¹²CO(J = 1–0) intensity (black dashed contours), Chandra X-ray boundary of N132D (black solid contours), and the H.E.S.S. teraelectronvolt gamma rays (white solid contours). An H i cloud appears projected onto the SNR that is elongated to the southwest direction with a hollow structure along the X-ray shell boundary. Three to four GMCs are located near the local intensity peaks of the H i cloud. One of them is possibly associated with the southern shell boundary of the SNR, which is consistent with previous CO studies (Banas et al. 1997; Desai et al. 2010; Sano et al. 2015b). Note that there are no dense molecular and atomic clouds toward the northeast outside the SNR. Further, note that teraelectronvolt gamma rays are emitted from the SNR itself rather than from the surrounding GMCs and H i cloud, even after taking into consideration the large PSF of the gamma-ray data.

Zoom In Zoom Out Reset image size

Figure 2(a) shows a red/green/blue (RGB) image of N132D obtained with Chandra. The X-ray shell shows an incomplete elliptical morphology, slightly elongated in the southwestern direction, with a breakout structure in the northeast. Many filamentary structures of X-rays appear not only in the shell boundary but also inside the SNR. The hard-band X-rays (E: 2.0–7.0 keV) are brighter in the southeastern shell.

Zoom In Zoom Out Reset image size

Figures 2(b) and (c) show the integrated intensity maps of ASTE ¹²CO(J = 3–2) and ATCA and Parkes H i. Using high-sensitivity and full-spatial-sampling observations of CO line emission, we found a molecular cloud toward the center of the SNR (hereafter "N132D MC-center"). Note that N132D MC-center is significantly detected because the CO integrated intensity of 1.12 K km s⁻¹ represents the ∼10σ level. The spatially resolved MC-center cloud is more extended than the beam size, with 10σ or higher significance in integrated intensity. We also confirm the presence of the previously identified GMC (hereafter "N132D GMC-south") in contact with the southeastern edge of the SNR. The peak velocities of the clouds are V_LSR ∼ 264 km s⁻¹ for N132D MC-center and V_LSR ∼ 266 km s⁻¹ for N132D GMC-south, and the latter is roughly consistent with the previous CO observations using SEST (Banas et al. 1997). On the other hand, the overall distribution of H i tends to encircle the X-ray shell except for northeast at the same velocity range of CO (V_LSR = 256.8–271.2 km s⁻¹). We also find that diffuse H i gas with an intensity of ∼300 K km s⁻¹ fills the interior of the X-ray shell.

Figure 3 shows a position–velocity diagram of CO and H i. We find an intensity dip at the velocity of ∼266 km s⁻¹, which is roughly centered at the position of the SNR in decl. On the other hand, the CO clouds appear projected onto the edge of the H i dip at the intensity level of 0.3 K degree (dashed contour centered at ∼266 km s⁻¹). Figure 4 shows averaged line profiles of CO and H i. The velocity range of the H i cloud at V_LSR ∼ 250–280 km s⁻¹ contains that of CO clouds at V_LSR ∼ 260–270 km s⁻¹. A strong absorption line of H i is detected at the velocity of V_LSR ∼ 266 km s⁻¹ toward only the SNR direction (blue, inside the SNR), suggesting that the absorption line was caused by the proximity of strong radio continuum radiation from the SNR (Yamane et al. 2018; Sano et al. 2018, 2019a). We therefore focus on both the CO and H i clouds around V_LSR ∼ 266 km s⁻¹ that are likely related to the SNR.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5 shows an RGB image of N132D composed of a combination of HST Hα (red), ALMA ¹²CO(J = 1–0) integrated intensity (green), and Chandra broadband X-rays (blue). We spatially resolved nine molecular clouds, named A to I, within the X-ray shell of N132D. Cloud A is located in the breakout region with very faint X-rays. Clouds H and I lie on the edge of the southwestern shell. The other clouds, B to G, corresponding to N132D MC-center, are concentrated in the center of the SNR. In other words, N132D MC-center is split into clouds B to G owing to high-resolution observations using ALMA. Note that clouds B, C, E, and F are located in the vicinity of Hα blobs or filaments, as shown in red.

Zoom In Zoom Out Reset image size

To derive the masses of these molecular clouds, we utilize the following equations:

$$\begin{eqnarray}&&M={m}_{{\rm{p}}}\mu {\rm{\Omega }}{D}^{2}\displaystyle \sum _{i}{N}_{i}({{\rm{H}}}_{2}),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&N({{\rm{H}}}_{2})=A\cdot W{[}^{12}\mathrm{CO}(J=1-0)],\end{eqnarray} \tag{ 2 }$$

where m_p is the mass of atomic hydrogen, μ = 2.72 is the mean molecular weight, Ω is the solid angle of each data pixel, D is the distance to the LMC (=50 kpc), ${N}_{i}({{\rm{H}}}_{2})$ is the column density of molecular hydrogen for each data pixel i, A is the CO-to-H₂ conversion factor, and $W{[}^{12}\mathrm{CO}(J=1-0)]$ is the integrated intensity of ¹²CO(J = 1–0) line emission. Here, we use the CO-to-H₂ conversion factor $A=7.0\times {10}^{20}$ cm⁻² (K km s⁻¹)⁻¹ (Fukui et al. 2008). The size of each molecular cloud is defined as an effective diameter, determined by the contour of the half-level of maximum integrated intensity. The detailed definitions and physical properties of molecular clouds are summarized in Table 2. The typical cloud masses and sizes are ∼50–100 M_⊙ and ∼1.5–2.0 pc, respectively. Note that the mass, n(H₂), and size of each cloud have ∼30% relative errors, due to uncertainties in the CO-to-H₂ conversion factor and distance to the LMC. There are no broad-line features with a velocity width more than 10 km s⁻¹, whereas the line widths of cloud A (ΔV = 5.7 km s⁻¹) and cloud B (ΔV = 4.4 km s⁻¹) are significantly larger than that of the other clouds (ΔV ∼ 1–2 km s⁻¹).

Table 2.  Physical Properties of Molecular Clouds Associated with N132D

Cloud Name	${\alpha }_{{\rm{J}}2000}$	${\delta }_{{\rm{J}}2000}$	${T}_{\mathrm{mb}}$	${V}_{\mathrm{LSR}}$	ΔV	Size	Mass	$n({{\rm{H}}}_{2})$
	(h m s)	(° ' '')	(K)	(km s−1)	(km s−1)	(pc)	(M⊙)	(cm−3)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
A .........	05 25 05.94	−69 37 52.9	0.79 ± 0.08	263.5 ± 0.3	5.7 ± 0.7	1.6	90	910
B .........	05 25 04.42	−69 38 17.8	0.52 ± 0.09	263.5 ± 0.4	4.4 ± 0.9	1.8	70	450
C .........	05 25 01.17	−69 38 14.1	0.74 ± 0.18	263.6 ± 0.1	1.1 ± 0.3	2.2	40	130
D .........	05 24 59.95	−69 38 22.1	1.34 ± 0.15	262.6 ± 0.1	1.5 ± 0.2	1.7	50	410
E .........	05 25 04.62	−69 38 31.6	1.86 ± 0.12	264.1 ± 0.1	1.7 ± 0.1	1.4	60	790
F .........	05 25 02.69	−69 38 35.8	2.79 ± 0.16	263.4 ± 0.1	1.7 ± 0.1	1.8	130	930
G .........	05 25 05.23	−69 38 51.2	1.39 ± 0.12	264.7 ± 0.1	2.1 ± 0.2	1.2	40	820
H .........	05 25 08.58	−69 38 57.0	1.66 ± 0.17	265.8 ± 0.1	0.9 ± 0.1	1.4	30	370
I .........	05 25 01.88	−69 39 16.1	2.95 ± 0.12	266.7 ± 0.1	2.0 ± 0.1	2.1	240	990

Notes. Column (1): Cloud name. Columns (2–9): Observed physical properties of the clouds obtained by single or double Gaussian fitting with ¹²CO(J = 1–0) emission line. Columns (2)–(3): Position of the clouds. Column (4): Maximum radiation temperature. Column (5): Central velocity of CO spectra. Column (6): FWHM line width of CO spectra ΔV. Column (7): Diameter of clouds defined as ${(S/\pi )}^{0.5}\times 2$, where S is the surface area of clouds surrounded by contours of the half-level of maximum integrated intensity. Column (8): Mass of clouds derived from equation $N({{\rm{H}}}_{2})/W(\mathrm{CO})=7.0\times {10}^{20}$ (K km s⁻¹)⁻¹ cm⁻², where $N({{\rm{H}}}_{2})$ is the column density of molecular hydrogen and W(CO) is the integrated intensity of ¹²CO(J = 1–0) (Fukui et al. 2008). Column (9): Number density of molecular hydrogen $n({{\rm{H}}}_{2})$.

Download table as:  ASCIITypeset image

Figure 6 shows an intensity ratio map of ¹²CO(J = 3–2)/¹²CO(J = 1–0) (hereafter R_3–2/1–0) using ASTE and ALMA, overlaid with Chandra X-ray contours. The intensity ratio reflects the CO rotational excitation states of the molecular clouds, and hence a high intensity ratio R_3–2/1–0 indicates a high temperature of the cloud. We find a high intensity ratio R_3–2/1–0 of ∼1.5–2.0 within the X-ray shell boundary. On the other hand, the intensity ratio R_3–2/1–0 of N132D GMC-south, south of the SNR, is ∼0.4, corresponding to the typical values of quiescent molecular clouds without any embedded OB association or shocks (e.g., Celis Peña et al. 2019).

Zoom In Zoom Out Reset image size

Figure 7(a) shows an enlarged view of the central region of N132D containing molecular clouds (green contours) and O-rich ejecta as seen by optical [O iii] emission (red) and O vii plus O viii images of X-rays (blue, see the Appendix and Figure 11). The optical [O iii] emission is especially bright toward clouds B and C, also known as Lasker's Bowl (Morse et al. 1996). The intercloud region between clouds D and F is also bright in both the [O iii] and O vii plus O viii emission, whereas no bright O-rich ejecta is detected toward the center of cloud F. We find no apparent trend between the spatial distributions of the molecular clouds and O-rich ejecta. Figure 7(b) shows radial velocity distributions of the O-rich knots presented by Law et al. (2020), overlaid with the ALMA CO contours. Major O-rich knots—B1, B2, B3, B4, R1, R2, and RK (runaway knot)—defined by Morse et al. (1995) are also indicated. We find that clouds B to F are projected onto the O-rich knots; cloud E is in contact with the blueshifted O-rich knots, whereas clouds B to D lie in the redshifted O-rich knots. Cloud F is possibly associated with both the blue and redshifted O-rich knots. It is noteworthy that cloud F shows a good spatial coincidence with the kinematic center of O-rich knots (marked with star symbol). Note that there are no CO counterparts of B1 and RK.

Zoom In Zoom Out Reset image size

To investigate the relationship between the molecular clouds and the X-ray emission of N132D, we derive the absorbing column densities toward two regions: one at the position of cloud F, and the other at a reference position south of cloud F (see Figure 7(a)). The absorbing column density is useful for constraining the origin of the X-ray emission. This might provide evidence for a possible positional relationship between the cloud and the shock front (Sano et al. 2015a, 2019a; Y. Yamane et al. submitted to ApJ).

We extracted ACIS-S3 spectra from the regions labeled as the "center region" and "southwest region" in Figure 7(a). We extracted the background spectrum from two rectangular regions with a total area of 2.8 arcmin² to provide sufficient statistics. One region was located to the southwest of the remnant and the other to the northeast; both were positioned to include the contribution from the transfer streak of the CCD. We use the background model of Sharda et al. (2020) to fit both the source and the background spectra for each region.

We used a two-component absorption model composed of the Milky Way absorption ${N}_{{\rm{H}},\mathrm{MW}}$ (TBabs) and the LMC absorption ${N}_{{\rm{H}},\mathrm{LMC}}$ (TBvarabs) by the ISM within the LMC along the line of sight. We fixed the hydrogen column density of the Milky Way at $1.47\times {10}^{21}$ cm⁻² (HI4PI Collaboration et al. 2016) with solar abundance (Wilms et al. 2000). Further, we set the elemental abundance for the ISM in the LMC with He = 0.9 Z_⊙ and the other elements at 0.5 Z_⊙ on the Wilms et al. (2000) scale.

Following the latest X-ray study of N132D (Sharda et al. 2020), we fitted each source spectrum with a plane-parallel shock model (vpshock, see Borkowski et al. 2001) plus two NEI components (vnei + vnei). The thermal emission from the forward shock along the outer rim has been modeled by Sharda et al. (2020) with a vpshock model with a temperature of ${{kT}}_{{\rm{e}}}=0.86\,\,\mathrm{keV}$ and an ionization timescale of $1.94\times {10}^{11}\,{\mathrm{cm}}^{-3}\,{\rm{s}}$. We fix these parameters at these values in our fits and only allow the normalization to vary. The vpshock component is intended to represent any emission from the forward shock that may contribute to our spectra along the line of sight. The two vnei components are intended to represent emission from a shock/cloud interface or shock-heated ejecta.

Figure 8 and Table 3 show the spectral fit results and the best-fit parameters, respectively. For the center region (cloud F), the C statistic is 1350 with 1302 degrees of freedom (DOF), and the Pearson reduced χ² is 0.99. The normalization of the vpshock component went to 0.0, and the LMC absorption went to $1.03\pm 0.10\times {10}^{21}$ cm⁻². One vnei component goes to a moderate temperature (∼0.82 keV) and has the abundances of O, Ne, Mg, Si, S, and Fe free to vary. There is marginal evidence for enhanced O, Ne, Mg, Si, and S abundances, but the Fe abundance is consistent with LMC values. The other vnei goes to a high temperature (3.36 keV) and has O, Ne, Mg, and Fe free. The Ne, Mg, and Fe abundances are significantly enhanced compared to mean local LMC values. Note that the ionization timescale for the 0.82 keV component goes to a value consistent with collisional ionization equilibrium (CIE; $3.5\times {10}^{13}$ cm⁻³ s), while the ionization timescale for the 3.36 keV component goes to a low value of $6.6\,\times {10}^{9}$ cm⁻³ s, indicating that the plasma producing this emission has been shocked relatively recently. For the southwest region (the reference region), the C statistic is 1427 with 1302 DOF, and the Pearson reduced χ² is 1.12. The vpshock component is now a significant contributor (see magenta line in Figure 8(b)). The LMC absorption is now 0.0, with an upper limit of $1.3\times {10}^{20}$ cm⁻². The fitted parameters for the two vnei components for the center and southwest regions are similar to each other; most values are within $1.0\sigma$ of each other. The major difference is the absence of the vpshock component in the center spectrum and the additional absorption for the center spectrum. This can be seen as the difference between the two spectra in the 0.35–1.0 keV energy range. The southwest spectrum has more emission at these lowest energies than the center spectrum, and this emission is modeled by the vpshock component. This indicates that there is additional absorption along the line of sight to the center region than toward the southwest region.

Zoom In Zoom Out Reset image size

Table 3.  Best-fit X-Ray Spectral Parameters

	Parameter	Center	Southwest
${N}_{{\rm{H}}}$	${N}_{{\rm{H}},\mathrm{LMC}}\,({10}^{21}\,{\mathrm{cm}}^{-2})$	${1.04}_{-0.11}^{+0.18}$	$\leqslant 0.13$
	${N}_{{\rm{H}},\mathrm{MW}}\,({10}^{21}\,{\mathrm{cm}}^{-2})$	1.47 (fixed)	1.47 (fixed)
vnei 1	${{kT}}_{e}\,(\mathrm{keV})$	${0.82}_{-0.02}^{+0.03}$	${0.79}_{-0.02}^{+0.04}$
	${Z}_{{\rm{O}}}\,(\mathrm{solar})$	${1.18}_{-0.48}^{+0.33}$	${1.92}_{-0.50}^{+0.29}$
	${Z}_{\mathrm{Ne}}\,(\mathrm{solar})$	${2.36}_{-0.27}^{+0.25}$	${1.88}_{-0.25}^{+0.16}$
	${Z}_{\mathrm{Mg}}\,(\mathrm{solar})$	${0.94}_{-0.13}^{+0.11}$	${0.70}_{-0.10}^{+0.12}$
	${Z}_{\mathrm{Si}}\,(\mathrm{solar})$	${0.86}_{-0.10}^{+0.08}$	${0.80}_{-0.09}^{+0.07}$
	${Z}_{{\rm{S}}}\,(\mathrm{solar})$	${0.70}_{-0.12}^{+0.11}$	${0.65}_{-0.06}^{+0.06}$
	${Z}_{\mathrm{Fe}}\,(\mathrm{solar})$	${0.30}_{-0.05}^{+0.04}$	${0.28}_{-0.07}^{+0.05}$
	${n}_{e}t\,({10}^{13}\,{\mathrm{cm}}^{-3}\,{\rm{s}})$	$\geqslant 3.50$	$\geqslant 3.50$
	norm (${10}^{-2}\,{\mathrm{cm}}^{-5}$)	${8.27}_{-0.76}^{+0.68}$	${13.17}_{-1.82}^{+1.12}$
vnei 2	${{kT}}_{e}\,(\mathrm{keV})$	${3.36}_{-0.46}^{+0.47}$	${2.44}_{-0.33}^{+0.49}$
	${Z}_{{\rm{O}}}\,(\mathrm{solar})$	${0.87}_{-0.18}^{+0.33}$	${0.64}_{-0.14}^{+0.22}$
	${Z}_{\mathrm{Ne}}\,(\mathrm{solar})$	${1.38}_{-0.30}^{+0.25}$	${1.70}_{-0.61}^{+0.36}$
	${Z}_{\mathrm{Mg}}\,(\mathrm{solar})$	${2.05}_{-0.20}^{+0.51}$	${2.58}_{-0.62}^{+0.79}$
	${Z}_{\mathrm{Fe}}\,(\mathrm{solar})$	${2.36}_{-0.47}^{+1.21}$	${6.38}_{-1.94}^{+4.43}$
	${n}_{e}t\,({10}^{10}\,{\mathrm{cm}}^{-3}\,{\rm{s}})$	${0.66}_{-0.06}^{+0.06}$	${0.47}_{-0.03}^{+0.01}$
	norm (${10}^{-3}\,{\mathrm{cm}}^{-5}$)	${4.65}_{-2.58}^{+1.16}$	${6.64}_{-2.47}^{+2.40}$
vpshock	${{kT}}_{e}\,(\mathrm{keV})$	0.86 (fixed)	0.86 (fixed)
	${n}_{e}t\,({10}^{11}\,{\mathrm{cm}}^{-3}\,{\rm{s}})$	1.94 (fixed)	1.94 (fixed)
	norm (${10}^{-2}\,{\mathrm{cm}}^{-5}$)	0.0	${3.03}_{-0.02}^{+0.02}$
	cstat (d.o.f.)	1350 (1302)	1427 (1302)
	Pearson-χ2 (reduced)	1287 (0.99)	1460 (1.12)

Download table as:  ASCIITypeset image

Figure 9 shows an overlay map of the 10–15 keV band image obtained with NuSTAR (colored image, Bamba et al. 2018) and the ALMA CO distribution in contours. The 10–15 keV band image represents possible synchrotron X-ray emission, although Bamba et al. (2018) could not exclude the possibility of a very high temperature plasma emission. We note that the hard X-ray emission is concentrated inside the SNR, where the molecular clouds B to G are located. Although the local intensity peaks of hard X-rays appear to be offset from the center of the molecular clouds, it is not certain whether the trend is significant because of the modest angular resolution of NuSTAR $\sim 18^{\prime\prime}$ in HPBW (∼4.4 pc at the LMC distance). It is certain that the edge of N132D is not bright with the hard X-ray emission, which is not typical for young SNRs with synchrotron X-rays (Bamba et al. 2005).

Zoom In Zoom Out Reset image size

In addition to the previously known GMC, which we refer to as N132D GMC-south, we identified eight new molecular clouds toward the center and southern edge of N132D. To better understand the relationship among the clouds, high-energy radiation, and O-rich ejecta in N132D, it is essential to know which clouds are physically associated with the SNR. Here, we argue that the eight new molecular clouds resolved by ALMA are likely interacting with shock waves and lie inside a wind-blown bubble.

We first claim that the high intensity ratio of R_3-2/1-0 > 1.5–2.0 as shown in Figure 6 provides strong evidence for shock–cloud interaction. The ratio of R_3-2/1-0 is useful for measuring the degree of rotational excitation of CO molecules, because the upper state of J = 3 lies at 33.2  K from the ground state of J = 0, corresponding to ∼28 K above the state of J = 1 at 5.5 K. The higher ratio of R_3-2/1-0 can trace warm molecular clouds heated by shock interactions not only for the Galactic SNRs (e.g., W28, Arikawa et al. 1999; Kesteven 79, Kuriki et al. 2018), but also for the Magellanic SNRs (e.g., LMC SNR N49, Yamane et al. 2018; SMC SNR RX J0046.5−7308, Sano et al. 2019b, 2019b). It is noteworthy that preshocked gas in N132D GMC-south shows a significantly lower intensity of R_3-2/1-0 ∼ 0.4, which is a typical ratio of a quiescent cloud in the LMC without external heating (e.g., Celis Peña et al. 2019).

We argue that cloud F has been completely engulfed by shocks and is located on the near side of the remnant. Figure 10 shows an enlarged view of the X-ray three-color image superposed on boundaries of molecular clouds. We find an X-ray filament toward cloud F. The color changes from red/yellow to green as one moves from east to west onto the cloud. This indicates that low-energy X-rays are suppressed toward cloud F by absorption. In fact, the LMC absorption ${N}_{{\rm{H}},\mathrm{LMC}}$ of cloud F ($1.04\times {10}^{21}$ cm⁻²) is significantly higher than that of the reference region without dense clouds ($\leqslant 0.13\times {10}^{21}$ cm⁻²). In this case, the forward shock likely propagated from behind cloud F to in front of it. Then, the X-ray filament was formed behind cloud F via shock interaction. This interpretation is also consistent with the absence of the vpshock component. If the shock wave is in the process of wrapping around the cloud, the thermal emission from the forward shock would be suppressed on the near side of the cloud as the shock re-forms on that side of the cloud. In addition, any thermal emission from the forward shock on the far side of the remnant is absorbed by the cloud. Both effects lead to a reduction in the thermal emission located along the line of sight to the center of the cloud.

Zoom In Zoom Out Reset image size

The fitted LMC absorption ${N}_{{\rm{H}},\mathrm{LMC}}$ of $\sim 1\times {10}^{21}$ cm⁻² toward cloud F also suggests that the molecular cloud is engulfed by shock waves. By using Equation (2), we can derive an average proton column density of cloud F to be $\sim 5\times {10}^{21}$ cm⁻², which is five times higher than the X-ray-derived value (see Section 4.6 and Table 3). This suggests that some of the X-ray emission originates in front of the cloud and not just behind it. The evaporating cloud scenario described in Cowie & McKee (1977) and White & Long (1991) and further explored by Zhang & Chevalier (2019) can produce such a morphology. A similar discussion is also applicable for cloud I in the southern edge of the SNR. According to Sharda et al. (2020), the LMC absorption ${N}_{{\rm{H}},\mathrm{LMC}}$ toward cloud I is ∼2–3 × 10²¹ cm⁻², two times lower than the average proton column density of cloud I. Although cloud I is located on the shell boundary of the SNR, it is likely that the shock is interacting with this cloud.

Additionally, the large ${n}_{{\rm{e}}}t$ value of the CIE plasma in cloud F is possibly consistent with a long elapsed time since the cloud was heated. Our ALMA observations revealed the molecular hydrogen density of cloud F to be 930 cm⁻³. Considering the postshocked gas density equals one-quarter of the preshock gas density in the limit of large Mach number (see Rankin-Hugoniot shock jump conditions), the electron density toward the region can be derived to ∼560 cm⁻³ assuming the electron-to-proton density ratio of 1.2. We then obtain the ionization time of ≳2000 yr, which is roughly consistent with the latest estimation of the SNR age of 2450 ± 195 yr (Law et al. 2020). We therefore propose a possible scenario that cloud F is completely engulfed by shocks soon after the supernova explosion.

We also argue that the molecular clouds we observe through ALMA were left behind inside a wind-blown bubble. According to Inoue et al. (2012), the surrounding ISM of a high-mass progenitor shows a highly inhomogeneous density distribution. The less dense gas such as H i clouds can be completely disrupted by the strong stellar winds, while the more dense gas such as molecular clouds can survive. As a result, a wind-blown bubble with a density of ∼0.01 cm⁻³ coexists with dense molecular clouds with density more than ∼10³ cm⁻³. For N132D, the wind cavity and dense clouds are seen in Figures 2(c) and 3. The H i cloud shows a cavity-like distribution in both the spatial and velocity planes. The expansion velocity of ∼6 km s⁻¹ is consistent with the typical gas motion seen in other core-collapse SNRs (e.g., Fukui et al. 2012; Kuriki et al. 2018; Sano et al. 2019b). We note that the wind-bubble explosion scenario is also proposed by previous optical and X-ray studies of N132D (Hughes 1987; Sutherland & Dopita 1995; Blair et al. 2000; Sharda et al. 2020) and modeled by Chen et al. (2003). Additionally, the estimated progenitor mass for N132D from the model of Chen et al. (2013) for an SNR evolving in a cavity is in good agreement with that from nucleosynthesis modeling of ejecta-rich regions (Sharda et al. 2020).

We emphasize that although the shock–cloud interaction has occurred, most of the dense molecular clouds can survive the shock erosion. When shock waves hit the dense clouds, the penetrating velocity can be described as ${V}_{\mathrm{sh}}/\sqrt{(n/{n}_{0})}$, where V_sh is the shock velocity before collision, n₀ is the ambient density inside the wind bubble (${n}_{0}\sim 0.01$ cm⁻³), and n is the number density of the molecular cloud ($n\sim 930$ cm⁻³ for the case of cloud F). Therefore, the shock waves in cloud F will be much decelerated to 1/300 V_sh, and hence the shock cannot penetrate into the cloud within a few thousand years. The numerical results also support this idea (e.g., Celli et al. 2019). Furthermore, evaporation of the shocked cloud is also negligible, due to the small thermal capacity of the SNR's shocks. Tatematsu et al. (1990) and Sasaki et al. (2006) discovered shocked molecular clouds in the middle-aged Galactic SNR G109.1−1.0 with the age of $\sim 1.4\times {10}^{4}\,\mathrm{yr}$ (Sasaki et al. 2013), which survived the encounter. The surviving clouds, with a total mass of 63 M_⊙, are associated with a thermal X-ray lobe. The authors conclude that the X-ray lobe was likely formed by the evaporation of a small outer portion of the clouds; the mass ratio of the thermal plasma relative to the molecular clouds is less than 10% (Sasaki et al. 2006). This interpretation is also supported by numerical simulations (Bolte et al. 2015). For N132D, clouds B to F have been partially evaporated because the bright Hα emission and thermal X-rays are seen in their vicinity (see Figures 2 and 5). In other words, almost all parts of the clouds in N132D likely survived shock erosion or evaporation, considering the large total cloud mass of ∼750 M_⊙ and young age of N132D (∼2500 yr). Therefore, N132D may be considered to represent an early stage of the mixed-morphology SNR (e.g., Rho & Petre 1998).

The hard X-ray enhancements around the shocked clouds also provide alternative evidence for an inhomogeneous density distribution of the ISM. When a supernova shock wave propagates into an inhomogeneous ISM with a density fluctuation of roughly 10⁵, the shock–cloud interaction may generate turbulence that enhances the magnetic field up to ∼1 mG at the surface of the shocked clouds (e.g., Inoue et al. 2009, 2012). This can be observed as synchrotron X-rays or a radio continuum enhancement around the shocked gas clouds (e.g., Sano et al. 2010, 2013, 2015a, 2017b, 2017a; Yamane et al. 2018; Okuno et al. 2018). On the other hand, the shock–cloud interaction will develop multiple reflected shock structures that can heat the gas up to high temperatures (e.g., Sano et al. 2019a). For N132D, the observational trend in Figure 9—a spatial correspondence between the 10–15 keV X-rays and molecular clouds B to F—shows possible evidence for the shock–cloud interaction with the magnetic field amplification or shock ionization. To confirm this idea, we need conclusive evidence of a synchrotron X-ray or high-temperature plasma enhancement around the molecular clouds with sufficient angular resolution. A deep exposure with Chandra offers the possibility of extracting the hard X-ray spectral component spatially coincident with the molecular clouds.

In conclusion, the eight new molecular clouds presented here are possibly located inside the wind-blown bubble formed by stellar winds from the progenitor of N132D, and these clouds are likely engulfed by supernova shock waves. On the other hand, the relationship between these clouds and the O-rich ejecta is still unknown from the current data set (see Section 4.5 and Figure 7). It is possible that the O-rich emission of optical and X-rays was efficiently produced by reverse shocks because of the shock interaction with the dense molecular clouds (e.g., Milisavljevic & Fesen 2015). Future ALMA observations with ∼0.1 pc resolution will allow us to compare spatial and kinematic distributions of the ISM/circumstellar medium (CSM) and ejecta.

N132D is thought to be a promising candidate for a hadronic gamma-ray emitter because of its bright teraelectronvolt gamma rays and very weak or absent synchrotron X-ray emission (Bamba et al. 2018). Although a detailed spatial comparison between the CO data and the gamma-ray emission could not be carried out, the presence of shocked molecular clouds provides support for the hadronic origin of gamma rays in N132D. Assuming that this hypothesis is correct, we derive the total energy of accelerated cosmic-ray protons, W_p, in N132D while taking into account the target gas density. Previous studies measured the values of W_p ∼ 10⁵⁰–10⁵¹ erg using the X-ray-based or model-dependent gas density (H.E.S.S. Collaboration et al. 2015; Bamba et al. 2018). Here, we reconsider the total energy of cosmic-ray protons in N132D using the neutral gas density, which is derived from radio observations.

It should be emphasized that teraelectronvolt gamma rays are emitted from the direction of N132D itself, rather than from the surrounding GMCs or H i cloud, even after taking into consideration the PSF of the gamma-ray image (Figure 1). This implies that the surrounding three to four GMCs and the southern H i cloud do not significantly contribute to the gamma-ray emission via the hadronic process.¹⁹ In the present study, we therefore focus on the target gas density within the shell of the SNR.

We estimate the target proton density within a wind-blown bubble. As discussed in Section 5.1, the molecular clouds are located inside the wind bubble. Since the intercloud density in the bubble is thought to be significantly low (∼0.01 cm⁻³, e.g., Weaver et al. 1977), the only plausible mechanism to produce gamma-ray emission concentrated in the center of the SNR is if the cosmic-ray protons interact with the molecular clouds. According to H.E.S.S. Collaboration et al. (2015), the total energy of cosmic-ray protons W_p can be described in the case of the hadron-dominant model as

$$\begin{eqnarray}&&{W}_{{\rm{p}}}\sim {10}^{52}{({n}_{{\rm{H}}}/1{\mathrm{cm}}^{-3})}^{-1}\,\,(\mathrm{erg}),\end{eqnarray} \tag{ 3 }$$

where n_H is the number density of interstellar protons. Adopting proton densities of molecular clouds of 260–1980 cm⁻³ (see Table 2), we then obtain W_p ∼ 0.5–3.8 × 10⁴⁹ erg. This is comparable to the values obtained for the Galactic gamma-ray SNRs (∼10⁴⁸–10⁴⁹ erg; e.g., Fukui et al. 2012, 2017; Yoshiike et al. 2013; Fukuda et al. 2014; Kuriki et al. 2018; Sano et al. 2019c). Note that the derived value gives a conservative lower limit on the total energy of cosmic-ray protons, because the hadronic gamma-ray emission can be observed only toward the gas cloud even if cosmic-ray protons have an azimuthally isotropic distribution. In other words, there are cosmic-ray protons that do not interact with the molecular clouds and do not produce gamma rays, and the value of W_p should be slightly increased. In any case, N132D can be classified as a common accelerator of cosmic-ray protons in the Local Group of galaxies.

We also discuss an alternative case that the shock front of N132D has reached the cavity wall of the wind-blown bubble. In this case, atomic hydrogen gas within the wind shell acts as the target for cosmic-ray protons. The column density of atomic hydrogen N_p(H i) is calculated using the following equation (Dickey & Lockman 1990):

$$\begin{eqnarray}&&{N}_{{\rm{p}}}({\rm{H}}\,{\rm\small{I}})=1.823\times {10}^{18}\cdot W({\rm{H}}\,{\rm\small{I}})\,({\mathrm{cm}}^{-2}),\end{eqnarray} \tag{ 4 }$$

where W(H i) is the integrated intensity of H i in units of K km s⁻¹. Since W(H i) toward N132D has a large uncertainty because of the radio continuum absorption (see Figures 3 and 4), we derive it by referring to the H i intensity surrounding the shell. The typical value of W(H i) near the shell is ∼500 K km s⁻¹ (see Figure 2(c)), and hence the average column density of atomic hydrogen is derived as $\sim 0.9\times {10}^{21}$ cm⁻². Considering the wind-bubble expansion, the atomic hydrogen gas was swept up within the thick wind shell. We here assume that the diameter and thickness of the H i wind shell are ∼25 pc and ∼5 pc, respectively. The former corresponds to an effective diameter of N132D, and the latter represents the typical thickness of a wind shell surrounding a high-mass star or core-collapse SNR (e.g., Yamamoto et al. 2006; Fukui et al. 2012, 2017; Sano et al. 2019c). We finally obtain the atomic hydrogen density within the wind shell to be ∼30 cm⁻³, corresponding to ${W}_{{\rm{p}}}\sim 3\times {10}^{50}$ erg. This energy is significantly higher than the values that are seen in Galactic gamma-ray SNRs, and hence N132D might be an energetic accelerator of cosmic-ray protons, as mentioned before (H.E.S.S. Collaboration et al. 2015; Bamba et al. 2018). To confirm the shock interaction with the wind shell, further H i observations are needed. The Australian Square Kilometre Array Pathfinder (ASKAP), MeerKAT, and the Square Kilometre Array (SKA) will be able to spatially resolve the wind-blown bubble of H i with fine angular resolution and high sensitivity.