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

0000-0003-2062-5692]H. Sano
National Astronomical Observatory of Japan, Mitaka, Tokyo 181-8588, Japan<EMAIL_ADDRESS>Department of Physics, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan
Institute for Advanced Research, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan

0000-0003-1415-5823]P. P. Plucinsky
Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA

0000-0003-0890-4920]A. Bamba
Department of Physics, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan
Research Center for the Early Universe, School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

0000-0003-3347-7094]P. Sharda
Research School of Astronomy and Astrophysics, Australian National University, Canberra, ACT 2611, Australia
Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia

0000-0002-4990-9288]M. D. Filipović
Western Sydney University, Locked Bag 1797, Penrith South DC, NSW 1797, Australia

0000-0003-1413-1776]C. J. Law
Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA

0000-0001-5609-7372]R. Z. E. Alsaberi
Western Sydney University, Locked Bag 1797, Penrith South DC, NSW 1797, Australia

0000-0001-8296-7482]Y. Yamane
Department of Physics, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan

0000-0002-2062-1600]K. Tokuda
National Astronomical Observatory of Japan, Mitaka, Tokyo 181-8588, Japan<EMAIL_ADDRESS>Department of Physical Science, Graduate School of Science, Osaka Prefecture University, 1-1 Gakuen-cho, Naka-ku, Sakai 599-8531, Japan

0000-0002-6606-2816]F. Acero
Service d’Astrophysique, CEA Saclay, F-91191 Gif-sur-Yvette Cedex, France

0000-0001-5302-1866]M. Sasaki
Dr. Karl Remeis-Sternwarte, Erlangen Centre for Astroparticle Physics, Friedrich-Alexander-Universit$\ddot{a}$t Erlangen-N$\ddot{u}$rnberg, Sternwartstra$\beta$e 7, D-96049 Bamberg, Germany

0000-0002-4708-4219]J. Vink
Anton Pannekoek Institute/GRAPPA, University of Amsterdam, P.O. Box 94249, 1090 GE Amsterdam, The Netherlands

T. Inoue
Department of Physics, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan

0000-0003-4366-6518]S. Inutsuka
Department of Physics, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan

0000-0003-3383-2279]J. Shimoda
Department of Physics, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan

0000-0002-2794-4840]K. Tsuge
Department of Physics, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan

K. Fujii
Department of Physics, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

F. Voisin
School of Physical Sciences, The University of Adelaide, North Terrace, Adelaide, SA 5005, Australia

0000-0003-2762-8378]N. Maxted
School of Science, University of New South Wales, Australian Defence Force Academy, Canberra, ACT 2600, Australia

0000-0002-9516-1581]G. Rowell
School of Physical Sciences, The University of Adelaide, North Terrace, Adelaide, SA 5005, Australia

0000-0001-7826-3837]T. Onishi
Department of Physical Science, Graduate School of Science, Osaka Prefecture University, 1-1 Gakuen-cho, Naka-ku, Sakai 599-8531, Japan

0000-0001-7813-0380]A. Kawamura
National Astronomical Observatory of Japan, Mitaka, Tokyo 181-8588, Japan<EMAIL_ADDRESS>
N. Mizuno
National Astronomical Observatory of Japan, Mitaka, Tokyo 181-8588, Japan<EMAIL_ADDRESS>
H. Yamamoto
Department of Physics, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan

0000-0002-1411-5410]K. Tachihara
Department of Physics, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan

Y. Fukui
Department of Physics, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan
Institute for Advanced Research, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan

N132D is the brightest gamma-ray supernova remnant (SNR) in the Large Magellanic Cloud (LMC). We carried out $^{12}$CO($J$ = 1–0, 3–2) observations toward the SNR using the Atacama Large Millimeter/submillimeter Array (ALMA) and Atacama Submillimeter Telescope Experiment. We find diffuse CO emission not only at the southern edge of the SNR as previously known, but also inside the X-ray shell. We spatially resolved nine molecular clouds using ALMA with an angular resolution of $5''$, corresponding to a spatial resolution of $\sim$1 pc at the distance of the LMC. Typical cloud sizes and masses are $\sim$2.0 pc and $\sim$100 $M_\odot$, respectively. High-intensity ratios of CO $J$ = 3–2 / 1–0 $> 1.5$ are seen toward the molecular clouds, indicating that shock-heating has occurred. Spatially resolved X-ray spectroscopy reveals that thermal X-rays in the center of N132D are produced not only behind a molecular cloud, but also in front of it. Considering the absence of a thermal component associated with the forward shock towards one molecular cloud located along the line of sight to the center of the remnant, this suggests that this particular cloud is engulfed by shock waves and is positioned on the near side of remnant. If the hadronic process is the dominant contributor to the gamma-ray emission, the shock-engulfed clouds play a role as targets for cosmic-rays. We estimate the total energy of cosmic-ray protons accelerated in N132D to be $\sim$0.5–$3.8 \times 10^{49}$ erg as a conservative lower limit, which is similar to that observed in Galactic gamma-ray SNRs.

§ INTRODUCTION

It has been a long-standing question how cosmic rays, consisting of mainly relativistic protons, are accelerated in interstellar space. Supernova remnants (SNRs) are promising candidates for acceleration sites of Galactic cosmic rays below the knee energy ($\sim$3$ \times 10^{15}$ eV), through the mechanism of diffusive shock acceleration <cit.> at their shocks. A conventional value of the total energy of Galactic cosmic-rays accelerated in a SNR is thought to be $\sim10^{49}$–$10^{50}$ erg, corresponding to $\sim$1–10% of the typical kinematic energy released by a supernova explosion <cit.>. One of the current challenges is to verify these predictions experimentally.

A young SNR (a few thousand years old) with bright TeV gamma-ray emission is a potential source for accelerating cosmic rays close to knee energy <cit.>. TeV gamma-rays from young SNRs can be generally produced by two different mechanisms: hadronic and leptonic processes <cit.>. For the hadronic process, interaction between cosmic-ray proton and interstellar proton creates a neutral pion that decays into two gamma-ray photons (refer to as “hadronic gamma-rays”). For the leptonic process, cosmic-ray electron energizes an interstellar photon to gamma-ray energy via inverse Compton scattering (refer to as “leptonic gamma-rays”). To establish the SNR origin of cosmic-ray protons, an observational detection of hadronic gamma-ray is needed. However, it is difficult to distinguish the hadronic/leptonic processes from spectral modeling alone <cit.>.

Investigating the interstellar gas associated with gamma-ray SNRs holds a key to solving this problem. If the hadronic process dominates, the gamma-ray flux is proportional to the number density of the interstellar gas assuming azimuthally-isotropic distribution of cosmic rays. This implies that the presence of gamma rays should be spatially coincident with the interstellar gas. [Fukui et al., 2012] demonstrated the spatial correspondence using CO/ datasets as interstellar molecular and atomic gas tracers and TeV gamma-ray data toward the Galactic SNR RX J1713.7$-$3946. The authors also derived the total energy of cosmic rays to be $\sim$10$^{48}$ erg, by adopting the number density of interstellar gas that interacts with the SNR. Subsequent studies toward the young TeV gamma-ray SNRs HESS J1731$-$347, Vela Jr., and RCW 86 in the Milky Way show similar values of $\sim10^{48}$–$10^{49}$ erg <cit.>. To better understand the origin of cosmic rays and their energy budget, we need to study not only Galactic SNRs, but also extragalactic sources such as gamma-ray bright SNRs in the nearby Large Magellanic Cloud (LMC) and the Small Magellanic Cloud (SMC).

§ OVERVIEW OF THE MAGELLANIC SNR N132D

N132D (LHA 120-N 132D) is the brightest X-ray and TeV gamma-ray SNR in the LMC <cit.>. The shell-like morphology of this remnant is clearly resolved in radio, infrared, optical, and X-ray wavelengths <cit.>, with a size of $114'' \times 90''$[See also, catalog papers of the LMC SNRs by [Badenes et al., 2010] and [Bozzetto et al., 2017].] or $\sim$25 pc at the distance of the LMC <cit.>. The age of the SNR is estimated to be $\sim$2500 yr [Morse et al., 1995, Hughes et al., 1998, Vogt & Dopita, 2011, Law et al., 2020]. N132D is also categorized as an Oxygen-rich (or O-rich) SNR, which most likely originated from a core-collapse supernova explosion <cit.>. Detailed optical studies revealed kinematic motions of O-rich ejecta using Doppler reconstructions with an average expansion velocity 1745 km s$^{-1}$ <cit.>.

Since the detection of TeV gamma ray emission associated with N132D, it has received much attention as a possible efficient accelerator of cosmic rays. The [H.E.S.S. Collaboration et al., 2015] first reported the significant detection of TeV gamma-rays toward three sources in the LMC, including the superbubble 30 Doradus C, and two SNRs N157B and N132D. The authors derived the 1–10 TeV gamma-ray luminosity of $(0.9 \pm 0.2) \times 10^{35}$ erg s$^{-1}$ for N132D at the assumed distance of 50 kpc, which is an order of magnitude higher than that of the young ($\sim$1600 yr) TeV gamma-ray SNR RX J1713.7$-$3946 in the Galactic plane. A subsequent GeV gamma-ray study using Fermi-LAT reported a 1–100 GeV gamma-ray luminosity of $\sim$10$^{36}$ erg s$^{-1}$, indicating that N132D is the brightest GeV gamma-ray SNR not only in the Magellanic Clouds, but also in the Local Group galaxies [Ackermann et al., 2016, Acero et al., 2016]. [Bamba et al., 2018] discovered hard X-ray emission ($E$: 10–15 keV) using NuSTAR. The authors derived an upper limit on the synchrotron X-ray flux of $2.0 \times 10^{35}$ erg s$^{-1}$ in the 2–10 keV band using Suzaku and NuSTAR, and argued that a high flux ratio of TeV gamma-ray and synchrotron X-rays is consistent with the hadronic origin of gamma-rays. However, to estimate the total energy of cosmic rays, the number density of interacting molecular and atomic clouds is needed.

N132D is also believed to be associated with a giant molecular cloud (GMC) that might be a possible target for cosmic-ray protons. [Banas et al., 1997] discovered a GMC toward the south of N132D using $^{12}$CO($J$ = 2–1) line emission with the Swedish-ESO Submillimeter Telescope (SEST). The GMC has a size of $\sim$22 pc and a virial mass of at least $\sim2 \times 10^5$ $M_\odot$. The authors suggested that a part of the GMC is interacting with the southern edge of the SNR. This interpretation was further supported by the presence of shock-heated dust components in the southeastern shell of N132D <cit.>. Subsequently, [Desai et al., 2010] and [Sano et al., 2015] presented a CO map using archival $^{12}$CO($J$ = 1–0) line emission data which was taken by the Mopra 22-m radio telescope as part of the Magellanic Mopra Assessment project <cit.>. A diffuse part of the GMC is possibly aligned with the southern shell of the SNR, while no dense clouds are found inside the shell. Owing to the modest angular resolution of the CO data of $\sim$23$''$–$45''$ (or $\sim$6–11 pc at the LMC distance) and lack of higher excitation line data (e.g., $^{12}$CO $J$ = 3–2, 4–3), there is no conclusive evidence for shock-heated molecular clouds in the existing data of this remnant.

In this study, we report new millimeter/submillimeter observations using $^{12}$CO($J$ = 1–0, 3–2) line emission with the Atacama Submillimeter Telescope Experiment (ASTE) and the Atacama Large Millimeter/submillimeter Array (ALMA). The high angular resolution of 5$''$ in the ALMA CO data will allow us to resolve molecular clouds illuminated by shock waves and cosmic-ray protons in N132D. Section <ref> gives details about the observations, data reductions, and archival data. Sections <ref> and <ref> show a large-scale view of the CO, , X-ray, and TeV gamma-ray emission; Section <ref> presents ALMA CO results and basic properties of the molecular clouds; Section <ref> discusses the excitation condition of the molecular clouds; Section <ref> gives a detailed comparison with the O-rich ejecta; Sections <ref> and <ref> present X-ray spectroscopy and a comparison with hard X-ray emission. Discussion and conclusions are given in Sections <ref> and <ref>, respectively.

§ OBSERVATIONS, DATA REDUCTIONS, AND ARCHIVAL DATA

§.§ CO

Observations of $^{12}$CO($J$ = 3–2) line emission at $\lambda$=0.87 mm wavelength were conducted in 2014 September 1–3 (PI: H. Sano, proposal# 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' \times 3.6'$ centered at ($\alpha_\mathrm{J2000}$, $\delta_\mathrm{J2000}$) $\sim$ ($05^\mathrm{h}25^\mathrm{m}2$$8$, $-69\degr38\arcmin35\arcsec$). 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 backend. 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 $\sim$111 km s$^{-1}$ and $\sim$0.11 km s$^{-1}$, respectively. The typical system temperature was $\sim$300 K, including the atmosphere in the single-side band. To derive the main beam efficiency, we observed N 159W <cit.>, obtaining a main beam efficiency of 0.67$\pm$0.08. We also observed the M-type AGB star R Dor every hour to satisfy pointing offsets accuracy within $2''$. After convolution with a two-dimensional Gaussian kernel, we obtained the cube data with the beam size of $\sim$$23''$ ($\sim$5.6 pc at the LMC distance). The typical noise fluctuations are $\sim$0.046 K at the velocity resolution of 0.4 km s$^{-1}$.

Observations of $^{12}$CO($J$ = 1–0) line emission at $\lambda$=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# 2013.1.01042.S). We utilized 40 antennas of 12-m array, 9 antennas of 7-m array, and 3 antennas of total power (TP) array. The effective observation area was a $150'' \times 150''$ rectangular region centered at ($\alpha_\mathrm{J2000}$, $\delta_\mathrm{J2000}$) $\sim$ ($05^\mathrm{h}25^\mathrm{m}2$$79$, $-69\degr38\arcmin34\farcs2$). The combined baseline length of 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 <cit.> package version 5.4.0. We used the multiscale CLEAN task implemented in the CASA package [Cornwell, 2008]. The scale parameters are $0''$, 318, and 954 for the 12-m array and $0''$, 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 — e.g., strong side-lobes — arising from poorly sampled areas near and beyond the maximum spatial frequency. We finally combined the cleaned data of 12- and 7-m array datasets and calibrated the TP array data by using the feather task. The final beam size of the feathered data is $5\farcs26 \times 4\farcs99$, with a position angle of $59\fdg60$, corresponding to a spatial resolution of $\sim$1.2 pc at the LMC distance. The typical noise fluctuations of the feathered data are $\sim$0.22 K at a velocity resolution of 0.4 km s$^{-1}$.

To investigate the CO gas distribution at larger spatial scales, we used the Magellanic Mopra Assessment Data Release 1 <cit.>. MAGMA is a $^{12}$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 $\sim$11 pc at the LMC distance. The typical noise fluctuations of a region surrounding N132D are $\sim$0.26 K at the velocity resolution of 0.53 km s$^{-1}$. We applied additional spatial smoothing with a two-dimensional Gaussian kernel. The angular resolution of smoothed data is $\sim$60$''$ ($\sim$15 pc at the LMC distance), which is the same resolution with the  survey data of the LMC (see Section <ref>).

To better understand the distribution of neutral atomic hydrogen toward N132D, we used an archival survey data of the  line at $\lambda$ = 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 Australia Telescope National Facility (ATNF). The angular resolution of the survey data is $60''$, corresponding to the spatial resolution of $\sim$15 pc at the LMC distance. The typical noise fluctuations of brightness temperature are $\sim$2.4 K at the velocity resolution of 1.689 km s$^{-1}$.

§.§ X-rays

We used archival X-ray data obtained by Chandra, for which the observation IDs are 5532, 7259, and 7266 (PI: K. J. Borkowski, proposal# 06500305), which have been published in previous papers <cit.>. The datasets were taken with the Advanced CCD Imaging Spectrometer S-array (ACIS-S3). Table <ref> lists the details of the observations. We utilized Chandra Interactive Analysis of Observations <cit.> software version 4.12 with CALDB 4.9.1 [Graessle et al., 2007] for data reduction. All the datasets were reprocessed using the chandra_repro task. We then created exposure-corrected, energy-filtered images using the fluximage task in the energy bands of 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 $\sim$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 $\sigma$ 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 $\chi^2$ statistic in the case of a low number of counts per spectral bin [Kaastra, 2017] and we report the Pearson $\chi^2$ (weighting by the model) to evaluate goodness of fit. We used ATOMDB version 3.0.9 [Foster et al., 2013] and non-equilibrium 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].

l c c c c

Chandra ACIS-S observation log of SNR N132D.
ObsID     Observation Date     Exposure     $\alpha_{\mathrm{J2000}}$     $\delta_{\mathrm{J2000}}$

(ks)     ($^{\mathrm{h}}$ $^{\mathrm{m}}$ $^{\mathrm{s}}$)     ($^{\circ}$ $\arcmin$ $\arcsec$)
05532   Jan 09, 2006   44.59     05 25 02.28     $-$69 38 37.32

07259   Jan 10, 2006   24.85     05 25 02.28     $-$69 38 37.32

07266   Jan 15, 2006   19.90     05 25 02.28     $-$69 38 37.32

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''$ (half-power beam-width; HPBW) or $58''$ (half-power diameter; HPD). To improve signal to noise ratios of the map, we smoothed the data with a two-dimensional Gaussian kernel of $20''$.

§.§ TeV Gamma-Rays

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

Integrated intensity map of the ATCA & Parkes  <cit.> overlaid with the Mopra $^{12}$CO($J$ = 1–0) integrated intensity (black dashed contours), Chandra X-ray boundary of N132D <cit.>, and the H.E.S.S. TeV gamma-ray excess counts <cit.>. The CO data have been spatially smoothed to match the FWHM of  ($\sim$$60''$). The integration velocity ranges of CO and  are from 240 to 290 km s$^{-1}$, which cover 74 % of total integrated intensity of . The lowest contour level and intervals of CO are 6.3 K km s$^{-1}$ ($\sim$9 $\sigma$) and 2.1 K km s$^{-1}$ ($\sim$3 $\sigma$), respectively. The contour level of X-ray boundary is $0.3 \times 10^{-6}$ counts pixel$^{-1}$ s$^{-1}$. The contour levels of TeV gamma-rays are 10, 14, and 18 excess counts. The dashed white rectangles indicate the observed areas of CO. The FWHM beam size of CO/ and the PSF of TeV gamma-rays are also shown in top right corner.

§.§ H$\alpha$ and [Oiii]

The optical data of H$\alpha$ and [Oiii] 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[hla.stsci.edu]. The observations were carried out using the ACS Wide Field Channel F658N for H$\alpha$ and the WFPC2 F502N for [Oiii], 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[www.stsci.edu/hst/HST_overview/documents/datahandbook].

§ RESULTS

(a) RGB X-ray image of N132D obtained with Chandra <cit.>. The red, green, and blue colors represent the energy bands of 0.5–1.2 keV, 1.2–2.0 keV, and 2.0–7.0 keV, respectively. (b) Intensity map of $^{12}$CO($J$ = 3–2) obtained with ASTE overlaid with the Chandra X-ray contours. The integration velocity range is from 256.8 km s$^{-1}$ to 271.2 km s$^{-1}$. The black contours represent the X-ray intensity. The contour levels of X-rays are 0.3, 0.6, 1.8, 3.6, 7.2, and $14.4 \times 10^{-6}$ counts pixel$^{-1}$ s$^{-1}$. (c) Intensity map of  obtained with ATCA and Parkes [Kim et al., 2003] overlaid with the Chandra X-ray contours. The integration velocity range and contour levels are the same as in Figure <ref>(b). The beam size and scale bar are also shown.

§.§ Large-Scale Distribution of CO, , X-ray, and TeV Gamma-Rays

Figure <ref> shows a map of ATCA & Parkes  intensity overlaid with the Mopra $^{12}$CO($J$ = 1–0) intensity (black dashed contours), Chandra X-ray boundary of N132D (black solid contours), and the H.E.S.S. TeV gamma-rays (white solid contours). An  cloud appears projected onto the SNR, which 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  cloud. One of them is possibly associated with the southern shell boundary of the SNR, which is consistent with previous CO studies <cit.>. Note that there are no dense molecular and atomic clouds toward Northeast outside the SNR. Further, note that TeV gamma-rays are emitted from the SNR itself rather than from the surrounding GMCs and  cloud, even after taking into consideration the large PSF of gamma-ray data.

§.§ CO and  Clouds toward the SNR

Figure <ref>(a) shows an 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.

Figures <ref>(b) and <ref>(c) show the integrated intensity maps of ASTE $^{12}$CO($J$ = 3–2) and ATCA & Parkes . Because of 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$^{-1}$ represents $\sim$10$\sigma$ level. The spatially-resolved MC-center cloud is more extended than the beam size, with 10$\sigma$ 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_\mathrm{LSR} \sim$264 km s$^{-1}$ for N132D MC-center and $V_\mathrm{LSR} \sim$266 km s$^{-1}$ 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  tends to encircle the X-ray shell except for northeast at the same velocity range of CO ($V_\mathrm{LSR} = 256.8$–271.2 km s$^{-1}$). We also find that diffuse  gas with an intensity of $\sim$300 K km s$^{-1}$ fills the interior of the X-ray shell.

Position-velocity diagram of H i. Superposed white solid contours indicate $^{12}$CO($J$ = 3–2) intensity. The integration range in Right Ascension is from $05^\mathrm{h}24^\mathrm{m}57^\mathrm{s}$ to $05^\mathrm{h}25^\mathrm{m}09^\mathrm{s}$. The lowest contour and contour intervals of CO are 0.003 K degree and 0.005 K degree, respectively. The dashed white line delineates an  cavity with the  intensity of 0.3 K degree. The beam size and velocity resolution are also shown in bottom right corner. Horizontal dashed lines represent observed boundaries of CO.

Averaged line profiles of $^{12}$CO($J$ = 3–2) (red) and  (blue and cyan) toward the SNR N132D. The blue and cyan spectra were extracted from regions of the inside and outside of the circle with a radius of $30''$ at the center position of ($\alpha_\mathrm{J2000}$, $\delta_\mathrm{J2000}$) $\sim$ ($05^\mathrm{h}25^\mathrm{m}02\fs88$, $-69\degr38\arcmin34\farcs8$). The red spectrum was produced by averaging the whole spectra shown in Figure <ref>b.

Figure <ref> shows a position-velocity diagram of CO and . We find an intensity dip at the velocity of $\sim$266 km s$^{-1}$, which is roughly centered at the position of the SNR in declination. On the other hand, the CO clouds appear projected onto the edge of the  dip at the intensity level of 0.3 K degree (dashed contour centered at $\sim$266 km s$^{-1}$). Figure <ref> shows averaged line profiles of CO and . The velocity range of  cloud at $V_\mathrm{LSR} \sim$250–280 km s$^{-1}$ contains that of CO clouds at $V_\mathrm{LSR} \sim$260–270 km s$^{-1}$. A strong absorption line of  is detected at the velocity of $V_\mathrm{LSR} \sim$266 km s$^{-1}$ toward only the SNR direction (blue, inside the SNR), suggesting that the absorption line was caused by vicinity of strong radio continuum radiation from the SNR <cit.>. We therefore focus on both the CO and  clouds around $V_\mathrm{LSR} \sim$266 km s$^{-1}$ that are likely related with the SNR.

RGB image of N132D obtained with the HST H$\alpha$ (red), ALMA $^{12}$CO($J$ = 1–0) (green), and Chandra X-rays in the energy band of 0.5–7.0 keV (blue). The integration velocity range is from 262.0 km s$^{-1}$ to 268.4 km s$^{-1}$. The contours represent the integrated intensity of CO, whose levels are 1.2 ($\sim3\sigma$), 1.5, 2.0, 3.0, 4.0, 6.0, 9.0, 13.0, 18.0, and 24.0 K km s$^{-1}$. The region enclosed by the dashed line indicates the observed area with ALMA. The CO clouds, named A–I, discussed in Section <ref> are indicated.

§.§ Detailed CO Distribution with ALMA

Figure <ref> shows an RGB image of N132D composed of a combination of HST H$\alpha$ (red), ALMA $^{12}$CO($J$ = 1–0) integrated intensity (green), and the Chandra broad-band X-rays (blue). We spatially resolved nine molecular clouds, named A–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–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–G owing to high-resolution observations using ALMA. Note that clouds B, C, E, and F are located in the vicinity of H$\alpha$ blobs or filaments as shown in red.

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

\begin{eqnarray}
M = m_{\mathrm{p}} \mu \Omega D^2 \sum_{i} N_i(\mathrm{H}_2),\\
N(\mathrm{H}_2) = A \cdot W[\mathrm{^{12}CO}(J = 1\text{--}0)],
\label{eq1}
\end{eqnarray}

where $m_\mathrm{p}$ is the mass of atomic hydrogen, $\mu = 2.72$ is the mean molecular weight, $\Omega$ is the solid angle of each data pixel, $D$ is the distance to the LMC (= 50 kpc), $N_i(\mathrm{H}_2)$ is the column density of molecular hydrogen for each data pixel $i$, $A$ is the CO-to-H$_2$ conversion factor, and $W[\mathrm{^{12}CO}(J = 1\text{--}0)]$ is the integrated intensity of $^{12}$CO($J$ = 1–0) line emission. Here, we use the CO-to-H$_2$ conversion factor $A = 7.0 \times 10^{20}$ cm$^{-2}$ (K km s$^{-1})^{-1}$ [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 <ref>. The typical cloud masses and sizes are $\sim$50–100 $M_\sun $ and $\sim$1.5–2.0 pc, respectively. Note that the mass, $n$(H$_2$), and size of each cloud have $\sim$30 % relative errors due to uncertainties in the CO-to-H$_2$ conversion factor and distance to the LMC. There are no broadline features with a velocity width more than 10 km s$^{-1}$, whereas the linewidths of cloud A ($\Delta V$ = 5.7 km s$^{-1}$) and cloud B ($\Delta V$ = 4.4 km s$^{-1}$) are significantly larger than that of the other clouds ($\Delta V$ $\sim$1–2 km s$^{-1}$).

Physical properties of molecular clouds associated with N132D

1cCloud name     $\alpha_{\mathrm{J2000}}$     $\delta_{\mathrm{J2000}}$     $T_{\mathrm{mb}} $     $V_{\mathrm{LSR}}$     $\Delta V$     Size      Mass     $n(\mathrm{H_2})$

($^{\mathrm{h}}$ $^{\mathrm{m}}$ $^{\mathrm{s}}$)     ($^{\circ}$ $\arcmin$ $\arcsec$)     (K)     0.9[1](km $\mathrm{s^{-1}}$)     0.9[1](km $\mathrm{s^{-1}}$)     (pc)      ($M_\sun $)     (cm$^{-3}$)

1c(1)     (2)     (3)     (4)     (5)     (6)     (7)     (8)     (9)
A .................     05 25 05.94     $-$69 37 52.9     $0.79 \pm 0.08$     $263.5 \pm 0.3$     $5.7 \pm 0.7$     1.6     090     910

B .................     05 25 04.42     $-$69 38 17.8     $0.52 \pm 0.09$     $263.5 \pm 0.4$     $4.4 \pm 0.9$     1.8     070     450

C .................     05 25 01.17     $-$69 38 14.1     $0.74 \pm 0.18$     $263.6 \pm 0.1$     $1.1 \pm 0.3$     2.2     040     130

D .................     05 24 59.95     $-$69 38 22.1     $1.34 \pm 0.15$     $262.6 \pm 0.1$     $1.5 \pm 0.2$     1.7     050    410

E .................     05 25 04.62     $-$69 38 31.6     $1.86 \pm 0.12$     $264.1 \pm 0.1$     $1.7 \pm 0.1$     1.4     060     790

F .................     05 25 02.69     $-$69 38 35.8     $2.79 \pm 0.16$     $263.4 \pm 0.1$     $1.7 \pm 0.1$     1.8     130     930

G .................     05 25 05.23     $-$69 38 51.2     $1.39 \pm 0.12$     $264.7 \pm 0.1$     $2.1 \pm 0.2$     1.2     040     820

H .................     05 25 08.58     $-$69 38 57.0     $1.66 \pm 0.17$     $265.8 \pm 0.1$     $0.9 \pm 0.1$     1.4     030     370

I ..................     05 25 01.88     $-$69 39 16.1     $2.95 \pm 0.12$     $266.7 \pm 0.1$     $2.0 \pm 0.1$     2.1     240     990

Col. (1): Cloud name. Cols. (2–9): Observed physical properties of the clouds obtained by single or double Gaussian fitting with $^{12}$CO($J$ = 1–0) emission line. Cols. (2)–(3): Position of the clouds. Col. (4): Maximum radiation temperature. Col. (5): Central velocity of CO spectra. Col. (6): FWHM linewidth of CO spectra $\Delta V$. Col. (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. Col. (8): Mass of clouds derived by an equation of $N(\mathrm{H_2}) / W(\mathrm{CO}) = 7.0 \times 10^{20}$ (K km s$^{-1}$)$^{-1}$ cm$^{-2}$, where $N(\mathrm{H_2})$ is the column density of molecular hydrogen and $W$(CO) is the integrated intensity of $^{12}$CO($J$ = 1–0) [Fukui et al., 2008]. Col. (9): Number density of molecular hydrogen $n(\mathrm{H_2})$

§.§ CO 3–2 / 1–0 Ratio

Intensity ratio map of $^{12}$CO($J$ = 3–2) / $^{12}$CO($J$ = 1–0) using ASTE and ALMA. The ALMA data was smoothed to match the effective beam size of the ASTE data (an effective beam size of $23''$). The beam size and scale bar are also shown in top right corner. The velocity range is same as in Figure <ref>. White dashed contours represent the X-ray intensity, whose contour levels are the same as in Figure <ref>b. The gray areas represent that the $^{12}$CO($J$ = 1–0) and/or $^{12}$($J$ = 3–2) data show the low significance of $\sim$8$\sigma$ or lower.

Figure <ref> shows an intensity ratio map of $^{12}$CO($J$ = 3–2) / $^{12}$CO($J$ = 1–0) (hereafter $R_{3\text{--}2/1\text{--}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 high intensity ratio $R_{3\text{--}2/1\text{--}0}$ indicates high temperature of the cloud. We find a high intensity ratio $R_{3\text{--}2/1\text{--}0}$ of $\sim$1.5–2.0 within the X-ray shell boundary. On the other hand, the intensity ratio $R_{3\text{--}2/1\text{--}0}$ of N132D GMC-south, south of the SNR, is $\sim$0.4, corresponding to the typical values of quiescent molecular clouds without any embedded OB association and/or shocks <cit.>.

§.§ Comparison with O-rich Ejecta

(a) RGB image of N132D obtained with the HST [O iii] (red), ALMA $^{12}$CO($J$ = 1–0) (green contours), and Chandra O viii emission produced using a novel method for source separation (blue, see Appendix). The integration velocity range and contour levels of CO are the same as in Figure <ref>. Regions enclosed by solid yellow lines are used for the X-ray spectral analysis (see Section <ref>). The dashed white box represents the area shown in Figure <ref>(b). (b) Radial velocity distributions of O-rich knots presented by [Law et al., 2020]. Superposed contours indicate the ALMA $^{12}$CO($J$ = 1–0), whose contour levels are the same as in Figure <ref>. The star and filled circle represent the center of O-rich knots and that of X-ray shell. Major knots — B1, B2, B3, B4, R1, R2, and RK (runaway knot) — from [Morse et al., 1995] are also indicated.

Figure <ref>a shows an enlarged view of the central region of N132D containing molecular clouds (green contours) and O-rich ejecta as seen by optical [Oiii] emission (red) and Ovii plus Oviii image of X-rays (blue). The optical [Oiii] 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 [Oiii] and Ovii plus Oviii 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 <ref>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–F are projected onto the O-rich knots; Cloud E is in contact with the blue-shifted O-rich knots, whereas clouds B–D lie in the red-shifted O-rich knots. Cloud F is possibly associated with both the blue and red-shifted O-rich knots. It is noteworthy that cloud F shows a good spatial coincidence with the kinematic center of O-rich knots (marked as star symbol). Note that there are no CO counterparts of B1 and RK.

§.§ X-ray Spectral Analysis

(a) ACIS-S spectra of the center (cloud F) and southwest (reference), with the best-fit models shown in the top panels. The magenta, green, red, and blue lines represent vpshock, vnei 1, vnei 2, and background components, respectively. The black lines represent the sum of the components. The bottom panels in each figure indicate (data $-$ model) / error.

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 <ref>a). The absorbing column density is useful to constrain the origin of the X-ray emission. This might provide evidence for a possible positional relationship between the cloud and the shock front <cit.>.

We extracted ACIS-S3 spectra from the regions labelled as the “center region” and “southwest region” in Figure <ref>a. We extracted the background spectrum from two rectangular regions with a total area of $2.8~\mathrm{arcmin^2}$ 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 comprised of the Milky Way absorption $N_\mathrm{H,MW}$ (TBabs) and the LMC absorption $N_\mathrm{H,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$^{-2}$ [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_\odot$ and the other elements = 0.5 $Z_\odot$ 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 <cit.> 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_\mathrm{e} = 0.86$ keV and an ionization timescale of $1.94\times10^{11}~\mathrm{cm^{-3}~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 and/or shock-heated ejecta.

Best-fit X-ray Spectral Parameters

$N_{\rm H}$     $N_{\rm H, LMC}~(10^{21}~\rm cm^{-2})$     1.04 $_{-0.11}^{+0.18}$     $\leq 0.13$

$N_{\rm H, MW}~(10^{21}~\rm cm^{-2})$     1.47 (fixed)     1.47 (fixed)

vnei 1     $kT_e~\rm (keV)$     0.82 $_{-0.02}^{+0.03}$     0.79 $_{-0.02}^{+0.04}$

$Z_{\rm O}~\rm (solar)$     1.18 $_{-0.48}^{+0.33}$     1.92 $_{-0.50}^{+0.29}$

$Z_{\rm Ne}~\rm (solar)$     2.36 $_{-0.27}^{+0.25}$     1.88 $_{-0.25}^{+0.16}$

$Z_{\rm Mg}~\rm (solar)$     0.94 $_{-0.13}^{+0.11}$     0.70 $_{-0.10}^{+0.12}$

$Z_{\rm Si}~\rm (solar)$     0.86 $_{-0.10}^{+0.08}$     0.80 $_{-0.09}^{+0.07}$

$Z_{\rm S}~\rm (solar)$     0.70 $_{-0.12}^{+0.11}$     0.65 $_{-0.06}^{+0.06}$

$Z_{\rm Fe}~\rm (solar)$     0.30 $_{-0.05}^{+0.04}$     0.28 $_{-0.07}^{+0.05}$

$n_et~\rm (10^{13}~cm^{-3}~s)$     $\geq 3.50$     $\geq 3.50$

norm ($10^{-2}~\rm cm^{-5}$)     8.27 $_{-0.76}^{+0.68}$     13.17 $_{-1.82}^{+1.12}$

vnei 2     $kT_e~\rm (keV)$     3.36 $_{-0.46}^{+0.47}$     2.44 $_{-0.33}^{+0.49}$

$Z_{\rm O}~\rm (solar)$     0.87 $_{-0.18}^{+0.33}$     0.64 $_{-0.14}^{+0.22}$

$Z_{\rm Ne}~\rm (solar)$     1.38 $_{-0.30}^{+0.25}$     1.70 $_{-0.61}^{+0.36}$

$Z_{\rm Mg}~\rm (solar)$     2.05 $_{-0.20}^{+0.51}$     2.58 $_{-0.62}^{+0.79}$

$Z_{\rm Fe}~\rm (solar)$     2.36 $_{-0.47}^{+1.21}$     6.38 $_{-1.94}^{+4.43}$

$n_et~\rm (10^{10}~cm^{-3}~s)$     0.66 $_{-0.06}^{+0.06}$     0.47 $_{-0.03}^{+0.01}$

norm ($10^{-3}~\rm cm^{-5}$)     4.65 $_{-2.58}^{+1.16}$     6.64 $_{-2.47}^{+2.40}$

vpshock     $kT_e~\rm (keV)$     0.86 (fixed)     0.86 (fixed)

$n_et~\rm (10^{11}~cm^{-3}~s)$     1.94 (fixed)     1.94 (fixed)

norm ($10^{-2}~\rm cm^{-5}$)     0.0     3.03 $_{-0.02}^{+0.02}$

cstat (d.o.f.)     1350 (1302)     1427 (1302)

Pearson-$\chi^2$ (reduced)     1287 (0.99)     1460 (1.12)

Figure <ref> and Table <ref> 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 $\chi^2$ 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$^{-2}$. One vnei component goes to a moderate temperature ($\sim$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$^{-3}$ s), while the ionization timescale for the 3.36 keV component goes to a low value of $6.6 \times 10^{9}$ cm$^{-3}$ 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 $\chi^2$ is 1.12. The vpshock component is now a significant contributor (see magenta line in Figure <ref>b). The LMC absorption is now 0.0, with an upper limit of $1.3 \times 10^{20}$ cm$^{-2}$. 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 the one toward the southwest region.

§.§ Comparison with hard X-ray emission

Map of hard X-ray emission ($E$: 10–15 keV) obtained with NuSTAR [Bamba et al., 2018]. The white and yellow contours indicate the ALMA $^{12}$CO($J$ = 1–0) and the boundary of X-ray shell. The contour levels and integration velocity range of CO are the same as in Figure <ref>. The beam size and scale bar are also shown in bottom right corner for each panel. The CO clouds A–I are also indicated.

Figure <ref> shows an overlay map of 10–15 keV band image obtained with NuSTAR <cit.> 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 are concentrated inside the SNR, where the molecular clouds B–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 $\sim18''$ in HPBW ($\sim$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].

§ DISCUSSION

§.§ Molecular Clouds associated with N132D

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 shockwaves and lie inside a wind-blown bubble.

We first claim that the high-intensity ratio of $R_{3\text{--}2/1\text{--}0} > 1.5$–2.0 as shown in Figure <ref> provides strong evidence for shock-cloud interaction. The ratio of $R_{3\text{--}2/1\text{--}0}$ is useful to measure 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 $\sim$28 K above the state of $J$ = 1 at 5.5 K. The higher ratio of $R_{3\text{--}2/1\text{--}0}$ can trace warm molecular clouds heated by shock interactions not only for the Galactic SNRs (e.g., W28, [Arikawa et al., 1999] [Arikawa et al., 1999]; Kesteven 79, [Kuriki et al., 2018] [Kuriki et al., 2018]), but also for the Magellanic SNRs (e.g., LMC SNR N49, [Yamane et al., 2018] [Yamane et al., 2018]; SMC SNR RX J0046.5$-$7308, [Sano et al., 2019] [Sano et al., 2019]). It is noteworthy that pre-shocked gas in N132D GMC-south shows significantly lower intensity of $R_{3\text{--}2/1\text{--}0} \sim$0.4, which is typical ratios of a quiescent cloud in the LMC without external heating <cit.>.

We argue that cloud F has been completely engulfed by shocks and is located on the near side of remnant. Figure <ref> 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 due to absorption. In fact, the LMC absorption $N_\mathrm{H,LMC}$ of cloud F ($1.04 \times 10^{21}$ cm$^{-2}$) is significantly higher than that of the reference region without dense clouds ($\leq0.13 \times 10^{21}$ cm$^{-2}$). 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 reforms 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.

RGB X-ray image toward the center of N132D. The red, green, and blue colors represent the energy bands of 0.35–0.85 keV, 0.85–1.6 keV, and 1.6–6.0 keV, respectively. Superposed contours indicate boundaries of ALMA CO clouds B–I as shown in Figure <ref>, whose contour levels are 1.2 K km s$^{-1}$. Note that cloud A is located outside the field of view. The dashed line indicates the X-ray filament, discussed in Section <ref>.

The fitted LMC absorption $N_\mathrm{H,LMC}$ of $\sim1 \times 10^{21}$ cm$^{-2}$ toward cloud F also suggests that the molecular cloud is engulfed by shockwaves. By using equation (<ref>), we can derive an average proton column density of cloud F to be $\sim5 \times 10^{21}$ cm$^{-2}$, which is five times higher than the X-ray derived value (see Section <ref> and Table <ref>). 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_\mathrm{H,LMC}$ toward cloud I is $\sim$2–$3 \times 10^{21}$ cm$^{-2}$, 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_\mathrm{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$^{-3}$. Considering the post-shocked gas density equals to 1/4 of the pre-shock gas density in the limit of large Mach number (c.f., Rankin-Hugnoiot shock jump conditions), the electron density toward the region can be derived to $\sim$560 cm$^{-3}$ assuming the electron-to-proton density ratio of 1.2. We then obtain the ionization time of $\ga 2000$ yr, which is roughly consistent with the latest estimation of the SNR age of $2450 \pm 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  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 $\sim$0.01 cm$^{-3}$ coexists with dense molecular clouds with density more than $\sim10^3$ cm$^{-3}$. For N132D, the wind cavity and dense clouds are seen in Figures <ref>c and <ref>. The  cloud shows a cavity-like distribution in both the spatial and velocity planes. The expansion velocity of $\sim$6 km s$^{-1}$ is consistent with the typical gas motion seen in other core-collapse SNRs <cit.>. We note that the wind-bubble explosion scenario is also proposed by previous optical and X-ray studies of N132D <cit.> 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_\mathrm{sh}$ is the shock velocity before collision, $n_0$ is the ambient density inside the wind-bubble ($n_0 \sim 0.01$ cm$^{-3}$), and $n$ is the number density of molecular cloud ($n \sim 930$ cm$^{-3}$ for the case of cloud F). Therefore, the shock waves in cloud F will be much decelerated to $1/300$ $V_\mathrm{sh}$, and hence the shock cannot penetrate into the cloud within a few thousand years. The numerical results also support this idea <cit.>. Furthermore, evaporation of the shocked cloud is also negligible due to the small thermal capacity of 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 $\sim1.4 \times 10^4$ yr [Sasaki et al., 2013], which survived the encounter. The surviving clouds, with a total mass of 63 $M_\odot$, 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 <cit.>. For N132D, clouds B–F have been partially evaporated because the bright H$\alpha$ emission and thermal X-rays are seen in the vicinity of them (see Figures <ref> and <ref>). In other words, almost all part of the clouds in N132D likely survived shock erosion or evaporation, considering the large total cloud mass of $\sim$750 $M_\odot$ and young age of N132D ($\sim$2500 yr). Therefore N132D may be considered to represent an early stage of the mixed-morphology SNR <cit.>.

The hard X-ray enhancements around the shocked clouds also provides alternative evidence for an inhomogeneous density distribution of the ISM. When a supernova shockwave propagates into an inhomogeneous ISM with a density fluctuation of roughly $10^5$, the shock-cloud interaction may generate turbulence that enhances the magnetic field up to $\sim$1 mG at the surface of the shocked clouds <cit.>. This can be observed as synchrotron X-rays or a radio continuum enhancement around the shocked gas clouds <cit.>. On the other hand, the shock-cloud interaction will develop multiple reflected shock structures which can heat the gas up to high temperatures <cit.>. For N132D, the observational trend in Figure <ref> — spatial correspondence between the 10–15 keV X-rays and molecular clouds B–F — shows possible evidence for the shock-cloud interaction with the magnetic field amplification and/or shock ionization. To confirm this idea, we need conclusive evidence of a synchrotron X-ray and/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 dataset (see Section <ref> and Figure <ref>). It is possible that the O-rich emission of optical and X-rays were efficiently produced by reverse shocks due to the shock interaction with the dense molecular clouds <cit.>. Future ALMA observations with $\sim$0.1 pc resolution will allow us to compare spatial and kinematic distributions of the ISM/circumstellar medium (CSM) and ejecta.

§.§ Is N132D the Energetic Accelerator of Cosmic-Ray Protons?

N132D is thought to be a promising candidate for a hadronic gamma-ray emitter because of its bright TeV 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_\mathrm{p}$, in N132D taking into account the target gas density. Previous studies measured the values of $W_\mathrm{p} \sim$10$^{50}$–$10^{51}$ 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 TeV gamma-rays are emitted from the direction of N132D itself, rather than from the surrounding GMCs and/or  cloud, even after taking into consideration the PSF of gamma-ray image (Figure <ref>). This implies that the surrounding three to four GMCs and the southern  cloud do not significantly contribute to the gamma-ray emission via the hadronic process[Note that future gamma-ray observatories with high-angular resolution and high sensitivity have a potential to detect gamma-rays from the escaped cosmic-ray protons from the surrounding GMCs and  cloud. The Cherenkov Telescope Array (CTA) can test the presence of high-energy particles escaped from N132D.]. 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 <ref>, the molecular clouds are located inside the wind bubble. Since the intercloud density in the bubble is thought to to be significantly low <cit.>, 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_\mathrm{p}$ can be described in case of the hadron-dominant model:

\begin{eqnarray}
W_\mathrm{p} \sim 10^{52} (n_\mathrm{H} / 1\; \mathrm{cm^{-3}})^{-1} \;\; \mathrm{(erg)},
\label{eq3}
\end{eqnarray}

where $n_\mathrm{H}$ is the number density of interstellar protons. Adopting proton densities of molecular clouds of 260–1980 cm$^{-3}$ (see Table <ref>), we then obtain $W_\mathrm{p} \sim$0.5–$3.8 \times 10^{49}$ erg. This is comparable to the values obtained for the Galactic gamma-ray SNRs <cit.>. 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_\mathrm{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_\mathrm{p}$() is calculated using the following equation [Dickey & Lockman, 1990];

\begin{eqnarray}
N_\mathrm{p}(\mathrm{H\,}{\textsc{i}}\xspace) = 1.823 \times 10^{18}  \cdot W(\mathrm{H\,}{\textsc{i}}\xspace) \;(\mathrm{cm}^{-2}),
\label{eq4}
\end{eqnarray}

where $W$() is the integrated intensity of  in units of K km s$^{-1}$. Since $W$() toward N132D has large uncertainty due to the radio continuum absorption (see Figures <ref> and <ref>), we derive it by referring to the  intensity surrounding the shell. The typical value of $W$() near the shell is to be $\sim$500 K km s$^{-1}$ (see Figure <ref>c), and hence the average column density of atomic hydrogen is derived to $\sim0.9 \times 10^{21}$ cm$^{-2}$. 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  wind shell are $\sim$25 pc and $\sim$5 pc, respectively. The former corresponds to an effective diameter of N132D, and the latter represents the typical thickness of wind shell surrounding a high-mass star or core-collapse SNR <cit.>. We finally obtain the atomic hydrogen density within the wind shell to be $\sim$30 cm$^{-3}$, corresponding to $W_\mathrm{p} \sim3 \times 10^{50}$ erg. This energy is significantly higher than the values which 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  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  with fine angular resolution and high sensitivity.

§ CONCLUSIONS

We have presented new $^{12}$CO($J$ = 1–0, 3–2) observations toward the LMC SNR N132D using ALMA and ASTE. The primary conclusions are summarized as follows.

* We have revealed the presence of diffuse CO emission inside the X-ray shell in addition to the previously known GMC at the southern edge of N132D. ALMA observations spatially resolved the diffuse CO emission into nine molecular clouds, whose sizes and masses are 1.2–2.2 pc and 30–240 $M_\odot$. High-intensity ratios of CO $J$ = 3–2 / 1–$0 > 1.5$ are seen toward the molecular clouds, indicating that shock-heating has occurred. The expansion  shell with an expanding velocity of $\sim$6 km s$^{-1}$ is also found toward N132D.

* Spatially resolved X-ray spectroscopy has revealed that the emission from the line of sight to cloud F can be well represented by a model with absorption in excess of the LMC absorption of $N_{\rm H}=1.04\times 10^{21}$ cm$^{-2}$ and two NEI thermal components (one of which approaches CIE conditions) and no thermal component for the forward shock emission. On the other hand, the fit to the X-ray spectrum of an adjacent region off of cloud F shows no additional absorption compared to the LMC value and requires a thermal component for the forward shock in addition to the two NEI components. The larger absorption and absence of a thermal component associated with the forward shock along the line of sight to cloud F suggests that cloud F has been engulfed by shocks and is located on the near side of remnant.

* We propose that the molecular clouds existed in the wind-blown bubble of the progenitor before the SNe explosion. The large $n_\mathrm{e}t$ value of one component of the plasma along the line of sight to cloud F is consistent with an elapsed time of $\ga 2000$ yr since the cloud was heated. The inhomogeneous density distribution inside the bubble — diffuse gas of $\sim$0.01 cm$^{-3}$ and dense clouds of $\sim$1000 cm$^{-3}$ — is also consistent with synchrotron X-ray and/or high-temperature plasma enhancement around the shocked clouds through the magnetic field amplification and/or shock ionization.

* If the hadronic process is the dominant contributor to the gamma-ray emission, the shock-engulfed molecular clouds play a role as targets for cosmic-rays. We estimate the total energy of cosmic-ray protons accelerated in N132D to be $\sim$0.5–$3.8 \times 10^{49}$ erg as a conservative lower limit, which is roughly the same value as seen in Galactic gamma-ray SNRs. The total energy could be as high as $\sim 3 \times 10^{50}$ erg if the shock front has reached the edge of the wind-blown cavity and the wind-shell of  has become a primary target for cosmic-ray protons. If the latter case is correct, N132D might be a very energetic accelerator of cosmic rays in the Local Group of galaxies.

CASA <cit.>, CIAO <cit.>, CALDB <cit.>, Xspec [Arnaud, 1996], SAO ds9 [Joye & Mandel , 2003].

§ ACKNOWLEDGMENTS

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.01042.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The ASTE radio telescope is operated by NAOJ. The Mopra radio telescope is part of the Australia telescope and is funded by the Commonwealth of Australia for operation as a National Facility managed by the CSIRO. The University of New South Wales Mopra Spectrometer Digital Filter Bank used for these Mopra observations was provided with support from the Australian Research Council, together with the University of New South Wales, the University of Adelaide, University of Sydney, Monash University, and the CSIRO. The scientific results reported in this article are based on data obtained from the Chandra Data Archive (Obs ID: 5532, 7259, and 7266). This research has made use of software provided by the Chandra X-ray Center (CXC) in the application packages CIAO (v 4.10). Based on observations made with the NASA/ESA Hubble Space Telescope, and obtained from the Hubble Legacy Archive, which is a collaboration between the Space Telescope Science Institute (STScI/NASA), the Space Telescope European Coordinating Facility (ST-ECF/ESA) and the Canadian Astronomy Data Centre (CADC/NRC/CSA). This research made use of the SAO/NASA Astrophysics Data System (ADS) bibliographic services. This work was supported by JSPS KAKENHI Grant Numbers JP16K17664 (H. Sano), JP19K14758 (H. Sano), JP19H05075 (H. Sano), and JP19K03908 (A. Bamba). H. Sano is also supported by the ALMA Japan Research Grant of NAOJ Chile Observatory (grant no. NAOJ-ALMA-244). This work is supported in part by Shiseido Female Researcher Science Grant (A. Bamba). P. Sharda is supported by the Australian Government Research Training Program (AGRTP) Scholarship. C. Law acknowledges funding from the National Science Foundation Graduate Research Fellowship under Grant DGE1745303. K. Tokuda was supported by NAOJ ALMA Scientific Research Grant Number of 2016-03B. M. Sasaki acknowledges support by the Deutsche Forschungsgemeinschaft through the Heisenberg professor grants SA 2131/5-1 and 12-1.

Spectral decomposition of the X-ray data cube into different spectral components with a novel deblending method. The position of the Ovii and Oviii lines at 0.574 and 0.654 keV respectively are shown by the dotted lines.

§ APPENDIX: AN X-RAY MAP FOR OXYGEN DOMINATED COMPONENT

Despite their multi-dimensional nature, X-ray data are most frequently analyzed as 2D images or 1D spectra independently therefore disconnecting the spatial and morphological information. To generate the oxygen X-ray map presented in Figure <ref>a we used a novel deblending technique recently adapted to X-ray data in [Picquenot et al., 2019] that takes the full advantage of the 3D (X, Y, E) information provided by X-ray spectro-imagers. This method (the Generalized Morphological Component Analysis ; GMCA) was initially developed to separate the cosmic micro-wave background image from the foregrounds in Planck data <cit.>. The general idea is to decompose the input X-ray data cube in a linear sum of images and associated spectra each component being different from the next one by its morphological and spectral signature. Note that the algorithm has a blind source separation approach and has no instrumental (instrument response) nor astrophysical (spectral emission) knowledge. Only the number of components to retrieve is fixed by the user. The main disentangling factor is the morphological diversity of each component in the wavelet domain[The wavelet transform is applied to each image slice of the data cube in order to enhance the contrast between small and large scale features.] and their associated spectral signatures.

Based on the same Chandra data set as presented in Section <ref>, a data cube was produced with the instrumental energy channel binning of 14.6 eV and a spatial bin size of 1.5 arcsec. The algorithm was applied in the 0.5 to 2.2 keV band and the number of components to retrieve was fixed to three. Figure <ref> shows the resulting spectral decomposition with one component dominating the low energies and exhibiting notable line emission at 0.574, 0.654 keV (dotted lines in Figure <ref>). Due to this spectral feature and the morphological similarities to the HST [Oiii] map (see Figure <ref>a), we associate this component with an oxygen rich component. The image associated with this spectral component is shown in Figure <ref>a (blue). As this component is dominating the 0.5–0.8 keV band, it is the most sensitive to absorption along the line of sight and the drops in flux in the image reflects the regions of highest absorption traced by the ALMA CO data.

[Acero et al., 2016] Acero, F., Ackermann, M., Ajello, M., et al. 2016, , 224, 8
[Ackermann et al., 2016] Ackermann, M., Albert, A., Atwood, W. B., et al. 2016, , 586, A71
[Aharonian et al., 1994] Aharonian, F. A., Drury, L. O., & Voelk, H. J. 1994, , 285, 645
[Arikawa et al., 1999] Arikawa, Y., Tatematsu, K., Sekimoto, Y., et al. 1999, , 51, L7
[Arnaud, 1996] Arnaud, K. A.  1996, Astronomical Society of the Pacific Conference Series, 101, 17A
[Badenes et al., 2010] Badenes, C., Maoz, D., & Draine, B. T. 2010, , 407, 1301
[Bamba et al., 2005] Bamba, A., Yamazaki, R., Yoshida, T., et al. 2005, , 621, 793
[Bamba et al., 2018] Bamba, A., Ohira, Y., Yamazaki, R., et al. 2018, , 854, 71
[Banas et al., 1997] Banas, K. R., Hughes, J. P., Bronfman, L., et al. 1997, , 480, 607
[Bell, 1978] Bell, A. R. 1978, , 182, 147
[Blair et al., 2000] Blair, W. P., Morse, J. A., Raymond, J. C., et al. 2000, , 537, 667
[Blandford & Ostriker, 1978] Blandford, R. D., & Ostriker, J. P. 1978, , 221, L29
[Bobin et al., 2015] Bobin, J., Rapin, J., Larue, A., et al. 2015, IEEE Transactions on Signal Processing, 63, 1199
[Bobin et al., 2016] Bobin, J., Sureau, F., & Starck, J.-L. 2016, , 591, A50
[Bolte et al., 2015] Bolte, J., Sasaki, M., & Breitschwerdt, D. 2015, , 582, A47
[Borkowski et al., 2001] Borkowski, K. J., Lyerly, W. J., & Reynolds, S. P. 2001, , 548, 820
[Borkowski et al., 2007] Borkowski, K. J., Hendrick, S. P., & Reynolds, S. P. 2007, , 671, L45
[Bozzetto et al., 2017] Bozzetto, L. M., Filipović, M. D., Vukotić, B., et al. 2017, , 230, 2
[Bykov et al., 2018] Bykov, A. M., Ellison, D. C., Marcowith, A. and Osipov, S. M.  2018, , 214, 41
[Cash, 1979] Cash, W. 1979, , 228, 939
[Celli et al., 2019] Celli, S., Morlino, G., Gabici, S., et al. 2019, , 487, 3199
[Celis Peña et al., 2019] Celis Peña, M., Paron, S., Rubio, M., et al. 2019, , 628, A96
[Chen et al., 2003] Chen, Y., Zhang, F., Williams, R. M., et al. 2003, , 595, 227
[Chen et al., 2013] Chen, Y., Zhou, P. and Chu, Y.-H.  2013, , 769, 16C
[Cornwell, 2008] Cornwell, T. J. 2008, IEEE Journal of Selected Topics in Signal Processing, 2, 793
[Cowie & McKee, 1977] Cowie, L. L., & McKee, C. F. 1977, , 211, 135
[Danziger & Dennefeld, 1976] Danziger, I. J., & Dennefeld, M. 1976, , 88, 44
[Desai et al., 2010] Desai, K. M., Chu, Y.-H., Gruendl, R. A., et al. 2010, , 140, 584
[Dickel & Milne, 1995] Dickel, J. R., & Milne, D. K. 1995, , 109, 200
[Dickey & Lockman, 1990] Dickey, J. M., & Lockman, F. J. 1990, , 28, 215
[Dopita et al., 2018] Dopita, M. A., Vogt, F. P. A., Sutherland, R. S., Seitenzahl, I. R., Ruiter, A. J. and Ghavamian, P.  2018, , 237, 10D
[Drury et al., 1994] Drury, L. O., Aharonian, F. A., & Voelk, H. J. 1994, , 287, 959
[Ezawa et al., 2004] Ezawa, H., Kawabe, R., Kohno, K., et al. 2004, , 763
[Fukuda et al., 2014] Fukuda, T., Yoshiike, S., Sano, H., et al. 2014, , 788, 94
[Fukui et al., 2008] Fukui, Y., Kawamura, A., Minamidani, T., et al. 2008, , 178, 56.
[Fukui et al., 2012] Fukui, Y., Sano, H., Sato, J., et al. 2012, , 746, 82
[Fukui et al., 2017] Fukui, Y., Sano, H., Sato, J., et al. 2017a, , 850, 71
[Foster et al., 2013] Foster, A. R., Ji, L., Yamaguchi, H., Smith, R. K., & Brickhouse, N. S. 2013, in American Institute of Physics Conference Series, Vol. 1545, American Institute of Physics Conference Series, ed. J. D. Gillaspy, W. L. Wiese, & Y. A. Podpaly, 252
[France et al., 2009] France, K., Beasley, M., Keeney, B. A., et al. 2009, , 707, L27
[Fruscione et al., 2006] Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, , 6270, 62701V
[Funk, 2015] Funk, S. 2015, Annual Review of Nuclear and Particle Science, 65, 245
[Graessle et al., 2007] Graessle, D. E., Evans, I. N., Glotfelty, K., et al. 2007, Chandra News 14, 33
[H.E.S.S. Collaboration et al., 2015] H.E.S.S. Collaboration, Abramowski, A., Aharonian, F., et al. 2015, Science, 347, 406
[H.E.S.S. Collaboration et al., 2018] H.E.S.S. Collaboration, Abramowski, A., Aharonian, F., et al. 2018a, , 612, A4
[H.E.S.S. Collaboration et al., 2018] H.E.S.S. Collaboration, Abdalla, H., Abramowski, A., et al. 2018b, , 612, A6
[H.E.S.S. Collaboration et al., 2018] H.E.S.S. Collaboration, Abdalla, H., Abramowski, A., et al. 2018c, , 612, A7
[HI4PI Collaboration et al., 2016] HI4PI Collaboration, Ben Bekhti, N., Floer, L., Keller, R., Kerp, J., Lenz, D., Winkel, B., Bailin et al.  2016, , 594A, 116H
[Hughes, 1987] Hughes, J. P.  1987, , 314, 103H
[Hughes et al., 1998] Hughes, J. P., Hayashi, I., & Koyama, K. 1998, , 505, 732
[Inoue et al., 2008] Inoue, H., Muraoka, K., Sakai, T., et al. 2008, ISSTT proceedings, in Ninteenth International Symposium on Space Terahertz Technology, ed. W. Wild (Groningen: Univ. Groningen Press), 281Space Terahertz Technology, ed. W. Wild (Groningen: Univ. Groningen Press), 281
[Inoue et al., 2009] Inoue, T., Yamazaki, R., & Inutsuka, S.-. ichiro . 2009, , 695, 825
[Inoue et al., 2012] Inoue, T., Yamazaki, R., Inutsuka, S.-i., & Fukui, Y. 2012, , 744, 71
[Joye & Mandel , 2003] Joye, W. A., Mandel, E.  2003, Astronomical Society of the Pacific Conference Series, 295, 489
[Kaastra, 2017] Kaastra, J. S. 2017, , 605, A51
[Kim et al., 2003] Kim, S., Staveley-Smith, L., Dopita, M. A., et al. 2003, , 148, 473.
[Kuriki et al., 2018] Kuriki, M., Sano, H., Kuno, N., et al. 2018, , 864, 161
[Lasker, 1978] Lasker, B. M. 1978, , 223, 109
[Lasker, 1980] Lasker, B. M. 1980, , 237, 765
[Law et al., 2020] Law, C. J., Milisavljevic, D., Patnaude, D. J., et al. 2020, , 894, 73
[Leahy et al., 2019] Leahy, D., Wang, Y., Lawton, B., Ranasinghe, S., Filipović, M. 2019 , 158, 149
[Maggi et al., 2016] Maggi, P., Haberl, F., Kavanagh, P. J., et al. 2016, , 585, A162
[Martin et al., 2019] Martin, P., Lu, C. -C., Voelk, H., Renaud, M., Filipovic, M., and Cherenkov Telescope Array Consortium, 2019, In Science with the Cherenkov Telescope Array, 7. KSP: Large Magellanic Cloud Survey, ed: Acharya, B. S. et al., p. 125-141, (World Scientific Publishing Co. Pte. Ltd.)
[McMullin et al., 2007] McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, adass XVI, 376, 127
[Milisavljevic & Fesen, 2015] Milisavljevic, D., & Fesen, R. A. 2015, Science, 347, 526
[Minamidani et al., 2011] Minamidani, T., Tanaka, T., Mizuno, Y., et al. 2011, , 141, 73
[Morse et al., 1995] Morse, J. A., Winkler, P. F., & Kirshner, R. P. 1995, , 109, 2104
[Morse et al., 1996] Morse, J. A., Blair, W. P., Dopita, M. A., et al. 1996, , 112, 509
[Ohira et al., 2012] Ohira, Y., Yamazaki, R., Kawanaka, N., et al. 2012, , 427, 91
[Okuno et al., 2018] Okuno, T., Tanaka, T., Uchida, H., et al. 2018, , 70, 77
[Picquenot et al., 2019] Picquenot, A., Acero, F., Bobin, J., et al. 2019, , 627, A139
[Pietrzyński et al., 2013] Pietrzyński, G., Graczyk, D., Gieren, W., et al. 2013, , 495, 76
[Rho & Petre, 1998] Rho, J., & Petre, R. 1998, , 503, L167
[Sano et al., 2010] Sano, H., Sato, J., Horachi, H., et al. 2010, , 724, 59
[Sano et al., 2013] Sano, H., Tanaka, T., Torii, K., et al. 2013, , 778, 59
[Sano et al., 2015] Sano, H., Fukui, Y., Yoshiike, S., et al. 2015a, Revolution in Astronomy with ALMA: The Third Year, 257
[Sano et al., 2015] Sano, H., Fukuda, T., Yoshiike, S., et al. 2015b, , 799, 175
[Sano et al., 2017] Sano, H., Yamane, Y., Voisin, F., et al. 2017a, , 843, 61
[Sano et al., 2017] Sano, H., Reynoso, E. M., Mitsuishi, I., et al. 2017b, Journal of High Energy Astrophysics, 15, 1
[Sano et al., 2018] Sano, H., Yamane, Y., Tokuda, K., et al. 2018, , 867, 7
[Sano et al., 2019] Sano, H., Matsumura, H., Nagaya, T., et al. 2019a, , 873, 40
[Sano et al., 2019] Sano, H., Rowell, G., Reynoso, E. M., et al. 2019b, , 876, 37
[Sano et al., 2019] Sano, H., Matsumura, H., Yamane, Y., et al. 2019c, , 881, 85
[Sasaki et al., 2006] Sasaki, M., Kothes, R., Plucinsky, P. P., et al. 2006, , 642, L149
[Sasaki et al., 2013] Sasaki, M., Plucinsky, P. P., Gaetz, T. J., et al. 2013, , 552, A45
[Schenck et al., 2016] Schenck, A., Park, S., & Post, S. 2016, , 151, 161
[Seok et al., 2013] Seok, J. Y., Koo, B.-C., & Onaka, T. 2013, , 779, 134
[Sharda et al., 2020] Sharda, P., Gaetz, T. J., Kashyap, V. L., and Plucinsky, P. P.  2020, , 894, 145
[Sorai et al., 2000] Sorai, K., Sunada, K., Okumura, S. K., et al. 2000, , 4015, 86
[Sutherland and Dopita, 1995] Sutherland, R .S., and Dopita, M. A.  1995, , 439, 365S
[Tappe et al., 2006] Tappe, A., Rho, J., & Reach, W. T. 2006, , 653, 267
[Tappe et al., 2012] Tappe, A., Rho, J., Boersma, C., et al. 2012, , 754, 132
[Tatematsu et al., 1990] Tatematsu, K., Fukui, Y., Iwata, T., et al. 1990, , 351, 157
[Verner et al., 1996] Verner, D. A., Ferland, G. J., Korista, K. T., et al. 1996, , 465, 487
[Vogt & Dopita, 2011] Vogt, F., and Dopita, M. A. 2011, , 331, 521
[Weaver et al., 1977] Weaver, R., McCray, R., Castor, J., et al. 1977, , 218, 377
[Williams et al., 2006] Williams, B. J., Borkowski, K. J., Reynolds, S. P., et al. 2006, , 652, L33
[White & Long, 1991] White, R. L., & Long, K. S. 1991, , 373, 543
[Wilms et al., 2000] Wilms, J., Allen, A., & McCray, R. 2000, , 542, 914
[Wong et al., 2011] Wong, T., Hughes, A., Ott, J., et al. 2011, , 197, 16
[Xiao & Chen, 2008] Xiao, X., & Chen, Y. 2008, Advances in Space Research, 41, 416
[Yamamoto et al., 2006] Yamamoto, H., Kawamura, A., Tachihara, K., et al. 2006, , 642, 307
[Yamane et al., 2018] Yamane, Y., Sano, H., van Loon, J. T., et al. 2018, , 863, 55
[Yoshiike et al., 2013] Yoshiike, S., Fukuda, T., Sano, H., et al. 2013, , 768, 179
[Zhang & Chevalier, 2019] Zhang, D., & Chevalier, R. A. 2019, , 482, 1602
[Zhu et al., 2019] Zhu, H., Slane, P., Raymond, J. and Tian, W. W.  2019, , 882, 135Z
[Zirakashvili & Aharonian, 2010] Zirakashvili, V. N. & Aharonian, F. A. 2010, , 708, 965