Does or did the supernova remnant Cassiopeia A operate as a PeVatron?

Within the current Galactic CR paradigm, SNRs are considered the major contributors to the observed CR flux. For decades, this conviction has been based on phenomenological arguments and supported by the Diffusive Shock Acceleration (DSA) theory (for a review, see Malkov & Drury (2001)). The dominance of SNRs’ contribution up to the ”knee”, a distinct break in the CR spectrum around 3-4 PeV, would imply that at certain epochs of their evolution, SNRs should operate as PeVatrons. However, as has been recognised long ago (Lagage & Cesarsky, 1983), within the framework of the standard DSA applied to young SNRs, the energy of accelerated protons could hardly exceed 100 TeV. A possible solution was proposed by Bell (2004) by amplifying the magnetic field upstream of the shock through the instabilities driven by CRs. Combined with two other critical conditions, (1) shock speed of thousands of km/s, and (2) particle diffusion in the extreme (Bohm) regime, enhancing the magnetic field to $B\geq 100\mu\rm G$ allows acceleration of protons and nuclei to 1 PeV/nuc.

The interactions of the shock accelerated protons and nuclei with the ambient gas make young SNRs potentially detectable sources of $\gamma$-rays and neutrinos. Therefore, the detection of TeV $\gamma$-rays from more than a dozen young SNRs ¹¹1see the online catalogue for TeV $\gamma$-ray sources, http://tevcat.uchicago.edu/ is a direct proof of effective acceleration of protons and/or electrons to very high energies. In general, the available data do not allow us to distinguish between the hadronic (“$\pi^{0}$-decay”) and leptonic (“inverse Compton”) origin of TeV $\gamma$-rays. Even ignoring this ambiguity, i.e., assuming that $\gamma$-rays are produced by accelerated protons, we face a problem with the explanation of TeV $\gamma$-ray measurements. The reported steep $\gamma$-ray spectra with a differential power-law index $\Gamma\sim 2.4-2.7$ (e.g., as summerised in the supplementary of Aharonian et al., 2019) do not readily agree with predictions of the standard DSA theory for strong shocks of young SNRs. This result can be interpreted in three different ways:

(1) A steep power-law proton spectrum is mimicked by the result of the combination of a hard ($E^{-2}$ type) power-law spectrum and an ”early” exponential cutoff, $E_{0}\ll 100$ TeV. This simple explanation supports the predictions of the standard DSA and the conclusion of Lagage & Cesarsky (1983) for the case of a non-amplified magnetic field. This mathematical trick works in a relatively narrow (one decade or so) energy interval around and below the cutoff energy. If this explanation is correct, at higher energies the measurements should reveal a gradual spectral steepening in the deep exponential cutoff region.

(2) The claim that strong shock acceleration results in a hard acceleration spectrum is not always correct. Indeed, steep power-law proton spectra can be formed over a wide energy interval in certain realistic environments and scenarios (see, e.g., Bell et al., 2011, 2019; Malkov & Aharonian, 2019; Hanusch et al., 2019; Caprioli et al., 2020; Xu & Lazarian, 2022). The steep power-law spectra can extend over a wide energy interval, formally up to 1 PeV, assuming that such energetic protons are effectively confined in the shell. Correspondingly, the resulting steep power-law $\gamma$-ray spectra could continue up to 100 TeV. The detection of steep multi-TeV $\gamma$-ray spectra in both above scenarios is a rather difficult but feasible task for detectors like CTA and LHAASO.

(3) Acceleration of protons to PeV energies at the very early ($\ll 100$ yr) epoch of the SNR evolution is preferable (Bell et al., 2013; Zirakashvili et al., 2014), in particular, because of the shock speed exceeding 10,000 km/s and the high density of the progenitor wind. However, the CRs with highest energy have already escaped from the shock upstream and thus the remnant is currently emptied of multi-TeV CRs. This can naturally explain the steep spectra of $\gamma$-rays from SNRs, including from very young ones, in particular from SNR Cassiopeia A (Aharonian et al., 2001; Ahnen et al., 2017; Acciari et al., 2010) and Tycho (Acciari et al., 2011).

The $\gamma$-ray spectrum in scenario (2) is predicted to have a steeper power law shape without cutoff, which is significantly different from the other two scenarios. The $\gamma$-ray spectra from the remnants in scenarios (1) and (3) could be similar, although for different reasons. On the other hand, scenario (3) differs principally from scenario (1) by the $\gamma$-radiation outside the remnant. While in scenario (1), PeV particles are not produced at all, scenario (3) allows proton acceleration up to 1 PeV, but assumes that presently the highest energy particles already left the remnant and presently occupy regions beyond the shell. Any outcome of $\gamma$-ray measurements at tens to hundreds TeV performed with detectors of adequate sensitivity - either detection of positive signals or upper limits - would provide definite answers to whether specific SNRs could accelerate protons 1 PeV at any epoch of their evolution.

The detection or upper limits from SNRs of different types and ages could finally establish whether SNRs, as a source population, can be responsible for the CR fluxes up to the ”knee”. The results would be relevant to both scenarios (2) and (3).

The historical SNR Cas A represents a key target for such studies. It has a well-defined age and distance. It also belongs to a class of core-collapse supernova that explodes in a dense red supergiant wind. Such systems are believed to be able to accelerate ions to the highest energy due to the efficient amplification of magnetic fields by CR streaming in dense environments.

In this paper, using LHAASO KM2A data, we derive $\gamma$-ray flux upper limits towards Cas A and test its ability to act as a PeVatron. In Sec.2, we discuss the status of previous gamma-ray observations of Cas A and the theoretical challenges of acceleration of protons to ultrahigh energies at different epochs of its evolution. In Sec.3, we review the physical properties of Cas A, including the gas distribution in proximity to the remnant. In Sec.4, we describe the procedure of extraction of flux upper limits based on the LHAASO KM2A data. In Sec.5, we discuss the astrophysical implications of the obtained results.

Cassiopeia A (Cas A, SNR G111.7-02.1) is one of the youngest SNRs in our Galaxy. This $\approx 340$ year old (Reed et al., 1995) structure is the remnant of a type IIb supernova (Krause et al., 2008) located at a distance of 3.4 kpc. It is one of the brightest galactic radio sources with a shell structure of angular radius 2.5’ corresponding to the physical size of 2.5 pc (Kassim et al., 1995). The synchrotron radiation extends from radio (Tuffs et al., 1997) to hard X-rays of about $100\leavevmode\nobreak\ \rm keV$ (Grefenstette et al., 2015). Fermi LAT observations reveal a hint of a hadronic origin of the $\gamma$-ray emissions (Yuan et al., 2013; Zirakashvili et al., 2014). TeV $\gamma$-ray emission from Cas A has been discovered by the HEGRA IACT array (Aharonian et al., 2001), and later confirmed by the MAGIC (Ahnen et al., 2017) and VERITAS (Acciari et al., 2010) collaborations.

The hadronic models of GeV-to-TeV radiation require a very high, but still acceptable, acceleration efficiency of about 25% (Zirakashvili et al., 2014), while in the two-zone models the acceleration efficiency required can be lower (Liu et al., 2022). The significant steepening or a cutoff in the spectrum reported around a few TeV by Ahnen et al. (2017) and Abeysekara et al. (2020) implies a corresponding steepening in the spectrum of parent protons at energies of tens of TeV. This constrains but does not exclude the presence of PeV protons in the nebula or outside the shell. Due to the limited age of Cas A, even if PeV protons have been accelerated at the early ($\ll 100$ years) epoch of the SNR, they could not propagate too far from the SNR. Thus, for robust conclusions, one should probe both the SNR itself and the $\approx 100$ pc environment surrounding Cas A, for $\gamma$-rays of energy exceeding 100 TeV.

Since we are searching for $\gamma$-rays produced by PeV CR protons through interactions with the ambient gas, the information about the gas is crucial for our studies. Cas A is located within a dense gas environment. Zhou et al. (2018) have detected molecular clouds of total mass $200\leavevmode\nobreak\ \rm M_{\rm\odot}$ in front of Cas A. Within 200 pc around Cas A, Kim et al. (2008) has derived the average gas density of about $3\leavevmode\nobreak\ \rm cm^{-3}$ by investigating the infrared emission. In this region, compact molecular clumps with a density as high as $10^{5}\leavevmode\nobreak\ \rm cm^{-3}$ have also been detected (Reynoso & Goss, 2002). Ma et al. (2019) have performed a detailed J = 1-0 survey of CO in a large field near Cas A using the Purple Mountain Observatory (PMO) 13.7 m millimetre telescope. The total mass of molecular gas was $\approx 9.5\times 10^{5}\leavevmode\nobreak\ \rm M_{\odot}$, which corresponds to an average density of $10\leavevmode\nobreak\ \rm cm^{-3}$ within 100 pc proximity of Cas A. The line profiles have asymmetric or broadened shapes, which indicates a possible interaction between the SNR and the molecular clouds. In this work, we used the J = 1-0 ${}^{12}\rm CO$ data as described in Ma et al. (2019). The derived molecular gas column densities are shown in Figure.1 by assuming $2.0\times 10^{20}\leavevmode\nobreak\ \rm(K\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1})^{-1}$ for the $\rm CO-to-H_{2}$ conversion factor (Bolatto et al., 2013).

As mentioned above, the average gas density is estimated as $10\leavevmode\nobreak\ \rm cm^{-3}$ in the vicinity, so to estimate the gas mass we need to calculate the volume that is occupied by CRs. The latter depends on the CR propagation in the vicinity of Cas A. Since we deal with relatively small spatial scales and ultra-high energies, we cannot a priori assume that the propagation proceeds in the diffusive regime. One cannot exclude that close to the accelerator, highest energy CRs propagate ballistically near the source and enter the diffusive regime when $r>D/c$, where $r$ is the distance from the source, $D$ is the diffusion coefficient and $c$ is the speed of light. The diffusion coefficient in the Interstellar Medium (ISM) is estimated as $D(E)\sim 2\times 10^{28}(\frac{E}{1\leavevmode\nobreak\ \rm GeV})^{0.5}\leavevmode\nobreak\ \rm cm^{2}/s$ (Yuan et al., 2017). Its extrapolation to $E=1\leavevmode\nobreak\ \rm PeV$ gives $D\sim 2\times 10^{31}\leavevmode\nobreak\ \rm cm^{2}/s$. This implies $r\gtrsim D/c\simeq$ 200 pc. For comparison, rectilinear propagation length of protons over 340 years is $\simeq$ 110 pc. Thus, even in the case of acceleration of PeV protons at the early epochs of the SNR, they would still continue to move ballistically. Due to the beaming effects, only the protons propagating towards us can produce visible $\gamma$-rays. Then the $\gamma$-ray source should have a similar angular size as the accelerator, i.e., the SNR shell. This is about $2.5^{\prime}$, so it can be detected by LHAASO-KM2A as a point-like source.

The diffusion of CRs in the vicinity of the accelerators could be strongly suppressed due to enhanced turbulence introduced, e.g., by the streaming instability of injected CRs (see e.g. (Malkov et al., 2013)). The interaction of the cosmic ray precursor and dense clumps can also drive turbulence. By either following the turbulent tangling field lines or mirror diffusion in compressible magnetic fluctuations (Lazarian & Xu, 2021), high-energy CR protons might undergo very slow diffusion. In this case, PeV CRs could propagate in the diffusive regime and the region occupied by VHE CRs would be of the order $r\sim\sqrt{2DT}\sim 20\leavevmode\nobreak\ \rm pc(\frac{\chi}{10^{-2}})^{0.5}$, where $\chi$ is the suppression factor of the diffusion coefficient near the source compared with the Galactic diffusion coefficient, $D_{source}(E)=\chi D(E)$. The angular size of the $\gamma$-ray emission is estimated as $r/d\sim 0.33^{\circ}(\frac{\chi}{10^{-2}})^{0.5}$. The maximum diffusion length, in any case, should be smaller than 110 pc, the ballistic propagation length over 340 years. And the maximum angular extension of the $\gamma$-ray emission when CRs are in the diffusive regime should be smaller than $1.8^{\circ}$. To estimate the total CR energy budget from the $\gamma$-ray luminosity, only the average density matters. As estimated in Ma et al. (2019), the average density in the nearby 100 pc of Cas A is $10\leavevmode\nobreak\ \rm cm^{-3}$. In this case, due to the lower average density and larger integration area in deriving the $\gamma$-ray flux, the derived CR energy budget would be more conservative than the ballistic propagation case we considered above. To derive the robust and model independent estimations we consider only the diffusive regime in the following discussion.

The data used in this analysis were collected from December 2019 to September 2022 by the half KM2A, three quarters KM2A and full KM2A. Quality selections have been applied to ensure a reliable data quality. We select the periods when at least 95$\%$ of working sub-detectors were in good condition. The total effective observation time is 932.74 days after data quality selection. To achieve a better performance of the array, several event cuts are also implemented. The cuts on data and event selections are the same as listed in the KM2A performance paper (Aharonian et al., 2021).

Considering the energy resolution and statistics, one decade of energy is divided into 5 bins with a bin width of $log_{10}E=0.2$. The sky in celestial coordinates (right ascension and declination) is divided into cells of size $0.1^{\circ}\times 0.1^{\circ}$ and filled with detected events according to their reconstructed arrival directions for each energy bin. To extract the excess of $\gamma$-rays from each cell, the “direct integration method” (Fleysher et al., 2004) is adopted to estimate the number of cosmic ray background events. The test statistic (TS) used to evaluate the significance of the source in this work is ${\rm TS}=2\log(\lambda)$, where $\lambda={\mathcal{L}_{1}/\mathcal{L}_{0}}$, $\mathcal{L}_{1}$ is the maximum likelihood value for the alternative hypothesis, and ${\mathcal{L}_{0}}$ is the maximum likelihood value of the null hypothesis. The TS map of the ROI is shown in the left panel of Fig.1. We found that there is no significant excess in the vicinity of Cas A.

As mentioned in the last section, we consider both the ballistic propagation and diffusive regime. In the former case, the expected $\gamma$-ray emission produced by CRs escaping from Cas A is point-like. In KM2A analysis, taking into account the point spread function (PSF), we integrate all photons in the 90% contamination radius from a point source. For diffusive regime, the intrinsic size of the source varies with the diffusion coefficient; however, the maximum angular radius is 1.8^∘ as estimated in the previous section. Taking into account of the angular resolution of the detector, we use an integration radius of 2.2^∘ for the lower energy bin and 1.85^∘ for energies higher than 100 TeV, to estimate the upper limits. No significant excess has been detected. We derived an upper limit for each energy bin at the 95% confidence level using the method proposed by Helene (1983). For a more conservative estimate of the upper limit, we use the CERN ROOT tools following Feldman & Cousins (1998) to calculate the 95% upper limit of the signal when the background count is less than 20. The results are shown in Fig.2 and Fig.3. The main uncertainties come from the uncertainty in the differential spectral index of the $\gamma$-ray emission, which we assume to derive the upper limit. In order to estimate such systematic uncertainties, the analysis is separately performed for the $\gamma$-ray spectral index of -2 or -3. The uncertainties found from this test are small: about 10% and 3% for the first two energy bins, and less than 1% at higher energies.

The systematic uncertainties affecting the flux upper limits are similar to those studied in Aharonian et al. (2021). The number of operating units varied with time due to debugging a few percent of detector units. The variable layout of the array, affecting the $\gamma$-ray/background separation, has a slight effect on the flux. The systematic uncertainty also comes from the atmospheric model used in the simulation, which affects the detection efficiency. The overall systematic uncertainty affecting the flux is estimated to be $\sim$7$\%$.

In principle, the diffuse galactic $\gamma$-ray emission (DGE) should be subtracted before we derive the upper limit for sources. But the inclusion of such a component in deriving the upper limit only has a marginal (10%) impact on the final results. Thus to avoid further systematic errors, we do not subtract the DGE in this work and use the most conservative upper limits.