A Serendipitous Discovery of GeV Gamma-Ray Emission from Supernova 2004dj in a Survey of Nearby Star-forming Galaxies with Fermi-LAT

As the supernova (SN) ejecta expands in the circumstellar medium, a collisionless shock forms. The collisionless shock accelerates electrons to high energies, which emit photons from the radio–submillimeter through X-ray energies, as have been observed from a large number of supernovae. Protons may also be accelerated by the shocks and they can produce GeV–TeV gamma-rays via the neutral pion decay (Berezhko et al. 2011, 2015; Murase et al. 2011; Li 2019; Wang et al. 2019). Gamma-rays can also be produced by nonthermal bremsstrahlung and inverse Compton radiation. Indeed, gamma-ray emission from GeV to TeV energies have been detected from a dozen of historic supernova remnants (SNRs) in our Galaxy. Although SNRs are distinct from the early supernova as they have already swept up a large amount of circumstellar medium (CSM), the physics governing the production of gamma-rays is the same. Thus, gamma-ray emission from supernovae, particularly from the nearest supernova SN 1987A, has been predicted for a long while (Berezhko et al. 2011, 2015). Recently, Malyshev et al. (2019) reported a recent enhancement of the GeV emission from the SN 1987A region as observed with Fermi-LAT. But the location of this source overlaps with several other potential gamma-ray sources, so the nature of the GeV emission remains to be clarified.

Gamma-ray emission from SNe IIn and from superluminous SNe have been searched by Ackermann et al. (2015) and Renault-Tinacci et al. (2018), respectively. No evidence for a signal was found, but their observational limits start to reach interesting parameter ranges expected by the theory. Recently, Yuan et al. (2018) reported the detection of a variable gamma-ray source spatially and temporally consistent with a peculiar supernova, iPTF14hls. However, there is a quasar in the error circle of the Fermi-LAT source, which is a blazar candidate according to the infrared data. The lack of multiwavelength observations of this quasar makes it difficult to conclusively address its connection with the gamma-ray variable source.

SN 2004dj is the nearest and brightest SN IIP exploded on about 2004 July 28 (UT) in the galaxy NGC 2403 at a distance of about 3.5 Mpc (Nakano et al. 2004; Patat et al. 2004; Vinkó et al. 2006; Nayana et al. 2018). The progenitors of SNe IIP are thought to be red supergiants. The fast moving stellar ejecta interacts with the CSM created by the stellar wind of the red supergiants. SN 2004dj was detected in a wide range of wavelengths from radio through infrared to X-rays during the first several years after the explosion (Meikle et al. 2011; Chakraborti et al. 2012; Nayana et al. 2018).

In this Letter, we report a serendipitous discovery of a GeV gamma-ray source whose position is compatible within uncertainties with the position of SN 2004dj. This discovery is made when we search for GeV emission from nearby star-forming galaxies. By analyzing the 11.4 yr of Fermi-LAT data, we find a GeV source from the direction of the nearby galaxy NGC 2403. The flux of the this source, however, does not obey the well-known relation between the gamma-ray luminosity and infrared luminosity for star-forming galaxies (Ackermann et al. 2012), disfavoring the usual cosmic ray–interstellar medium (ISM) interaction origin. We also find the flux is decaying during the Fermi-LAT observation period from 2008 to the present. Motivated by this, we search for transient sources in the error region of the gamma-ray emission and find that SN 2004dj lies within the error region.

We perform a search for possible gamma-ray emission from galaxies in the IRAS Revised Bright Galaxies Sample (Sanders et al. 2003), using 11.4 yr of gamma-ray data taken by the Fermi-LAT telescope (Atwood et al. 2009). This is a complete flux-limited sample of all extragalactic objects brighter than 5.24 Jy at 60 μm, covering the entire sky surveyed by IRAS at Galactic latitudes $| b| \gt 5^\circ$. In the sample, 15 infrared (IR)-bright galaxies have been detected in gamma-rays with Fermi-LAT and listed in Fermi-LAT Fourth Source Catalog (4FGL; Abdollahi et al. 2020), including six galaxies with AGNs (i.e., Cen A, IC 4402, NGC 3067, NGC 3683, NGC 1275, NGC 3424; Peng et al. 2019), seven star-forming galaxies (SMC, LMC, M31, NGC 253, M82, NGC 2146, Arp 220; Abdo et al. 2010a; Ackermann et al. 2012; Tang et al. 2014; Griffin et al. 2016; Peng et al. 2016), and two star-forming galaxies with obscured AGNs (NGC 1068, NGC 4945; Ackermann et al. 2012). These galaxies are excluded, leading to a sample of 614 galaxies.

This work uses ∼11.4 yr (MET 239557417-595385929) of Fermi-LAT Pass8 SOURCE class events with reconstructed energies between 300 MeV and 500 GeV, excluding those with a zenith angle larger than 90° to avoid Earth limb contamination. We implement a standard sequence of analysis steps. For each IR-bright galaxy, we select the events in a 17° × 17° region of interest (ROI) centered at the galactic center and use the gtmktime tool to select time intervals expressed by (DATA_QUAL > 0) && (LAT_CONFIG==1). We bin the data in 20 logarithmically spaced bins in energy and in a spatial bin of 0025 per pixel. We make use of recent developments of the Science Tools for likelihood analysis. The background model for each celestial ROI contains all sources listed in the 4FGL along with the standard diffuse emission background, i.e., the foreground for Galactic diffuse emission ($\mathrm{gll}\_\mathrm{iem}\_{\rm{v}}7.\mathrm{fits}$) released and described by the Fermi-LAT collaboration through the Fermi Science Support Center (Acero et al. 2016) and the background for the spatially isotropic diffuse emission with a spectral shape described by $\mathrm{iso}\_{\rm{P}}8{\rm{R}}3\_\mathrm{SOURCE}\_{\rm{V}}2\_{\rm{v}}01.\mathrm{txt}$. We ignore the energy dispersion, which is not important for the analysis above 300 MeV. We employ the gttsmap tool to evaluate the 6° × 6° map of the test statistic (TS), defined as $\mathrm{TS}=-2(\mathrm{ln}{L}_{0}-\mathrm{ln}L)$, where L₀ is the maximum-likelihood value for null hypothesis and L is the maximum-likelihood with the additional point source with a power-law spectrum, and then survey the new gamma-ray sources when the significance of the gamma-ray excess above backgrounds is TS > 25. The best location and uncertainty of the new sources can be determined by maximizing the TS value and using the distribution of the localization test statistic (LTS), defined by twice the log of the likelihood ratio of any position with respect to the maximum. We claim a gamma-ray source associated with a target galaxy when the target galaxy lies within the 95% confidence location region of the gamma-ray source. Note that new background sources could be found and, if so, we redo the analysis using an updated background model including the new background sources. In the likelihood analysis, we allow all the sources that are separated from the ROI center by less than 65 to have a free normalization and fixed index.⁷ The normalizations of both Galactic and extragalactic diffuse emission models are left free.

Our analysis results in the detection of two new gamma-ray sources that are, respectively, spatially coincident with the galaxy Arp 299 and the galaxy M33 (see Table 1; more details can be found in Xi et al. 2020). In addition, we find that the gamma-ray source reported in 4FGL (4FGL J0737.4+6535) is spatially coincident with the galaxy NGC 2403, since the best location of 4FGL J0737.4+6535 is within the optical disk (about 8 kpc in diameter) of NGC 2403. For the rest of the galaxies in our sample, we do not find any significant emission.

Table 1.  New Sources Detected in the Survey

Name	dL	Optical Center	Gamma-Ray Location	${F}_{0.1\mbox{--}100\mathrm{GeV}}$	Γ	TS
	(Mpc)	(R.A., Decl.)	(R.A., Decl.)	(10−12 erg cm−2 s−1)
(1)	(2)	(3)	(4)a	(5)b	(6)	(7)
SN 2004dj	3.49	(114321, 65599)	(114395, 65573) ± 0041c	2.29 ± 0.51c	1.92 ± 0.07c	73.1c
Arp 299	47.74	(172136, 58561)	(172050, 58526) ± 0111	1.08 ± 0.28	2.07 ± 0.20	27.8
M33	0.84	(23475, 30669)	(23609, 30784) ± 0089	1.28 ± 0.42	2.23 ± 0.24	25.1

Notes. (1) Source name; (2) optical R.A. J2000; (3) optical decl. J2000; (4) best location of the γ-ray source detected by Fermi-LAT; (5) 100 MeV–100 GeV γ-ray average fluxes; (6) power-law spectral photon index derived by broad band spectrum fitting; (7) TS value of the gamma-ray sources.

^aThe uncertainties of the position correspond to 95% containment radius, which are derived by fitting the distribution of LTS to a 2D Gaussian function. ^bThe fluxes and errors were computed in the 0.1–100 GeV band, from our analysis in the 0.3–500 GeV band, taking into account the covariance matrix of the maximum-likelihood fit to the data. ^cThese values are derived for the first 5.7 yr period (2008 August 4 (UT)–2014 March 25 (UT)).

Download table as:  ASCIITypeset image

It has been found that there is an empirical correlation between the gamma-ray luminosity in (0.1–100) GeV and total IR luminosity (8–1000 μm) for Local Group galaxies and nearby star-forming galaxies (Ackermann et al. 2012; Peng et al. 2016). This correlation is generally interpreted as that the gamma-ray emission arises from cosmic-ray (CR) protons interacting with the ISM via the proton–proton (pp) collisions (Ackermann et al. 2012; Zhang et al. 2019). In Figure 1, we show the relation between the gamma-ray luminosity and IR luminosity for the three gamma-ray sources found in our work. While Arp 299 and M33 are well consistent with the aforementioned correlation, NGC 2403 is obviously an outlier. The latter is even above the theoretical calorimetric limit, which is obtained by considering that all of the cosmic-ray energy is converted into secondary particles and may be achieved only in those gas-rich starburst galaxies with luminous IR radiation (e.g., L_{8–1000 μm} > 10¹¹ L_⊙). This suggests that the gamma-ray emission from NGC 2403 cannot primarily arise from the CR–ISM interaction in the galaxy.

Zoom In Zoom Out Reset image size

We generate the light curves of the target gamma-ray source with eight and four time bins, respectively. In each time bin, all sources (including the target source) within 65 of the ROI center have spectra fixed to the shapes obtained from the full data set analysis, and only normalizations are allowed to vary. As shown in Figure 2, the fluxes appear to decay with time in both cases. We then use a likelihood-based statistic to test the significance of the variability. Following the definition in 2FGL (Nolan et al. 2012), the variability index from the likelihood analysis is constructed, with a value in the null hypothesis where the source flux is constant across the full time period, and the value under the alternate hypothesis where the flux in each bin is optimized: ${\mathrm{TS}}_{\mathrm{var}}={\sum }_{i=1}^{N}2\times (\mathrm{log}({L}_{i}({F}_{i}))-\mathrm{log}({L}_{i}({F}_{\mathrm{mean}})))$, where L_i is the likelihood corresponding to bin i, F_i is the best-fit flux for bin i, and F_mean is the best-fit flux for the full period assuming a constant flux. The statistic TS_var is expected to be distributed, in the null case, as ${\chi }_{N-1}^{2}({\mathrm{TS}}_{\mathrm{var}})$. We find that the gamma-ray emission is variable at a confidence level of 3.3σ for the analysis with four time bins and of 2.7σ for the analysis with eight time bins, respectively. The variability can be also seen from Figure 3, which shows an obvious gamma-ray excess over the background using the first 5.7 yr of the Fermi-LAT observation, but no significant excess in the second 5.7 yr. In this case, the most significant variability is obtained at the 3.5σ (TS_var = 13.0) significance level. We can estimate a post-trial significance level at 3.1σ considering the number of trials is three as we test the variability using two, four, and eight time bins respectively. Thus, we suggest that the gamma-ray source 4FGL J0737.4+6535 in the region of NGC 2403 is fading with time.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We use the catalogs from the NASA/IPAC Extragalactic Database (NED) and SIMBAD Database to search for possible counterparts of the gamma-ray source. Due to the nondetection of gamma-ray emission during the second half of the interval of Fermi observations, we use the data of the first 5.7 yr only to relocalize the gamma-ray emission, which is shown in Figure 3. We do not find any Galactic sources (i.e., novae, pulsars, gamma-ray binaries) or promising extragalactic gamma-ray emitters such as blazars located within the 95% error region of the gamma-ray emission. Two faint radio sources (i.e., NVSS J073724+653628 and NGC 2403:[ECB2002] alpha) with an unknown nature are located within the 95% error region. The 1.4 GHz flux densities are 4 ± 0.6 mJy for NVSS J073724+653628 (Condon et al. 1998) and ∼1.9 mJy for NGC 2403:[ECB2002] alpha (Eck et al. 2002). If the radio sources are background radio galaxies, which constitute the dominant mJy radio source population at 1.4 GHz (Norris 2017), we can estimate the gamma-ray flux assuming that the gamma-ray emission arises from inverse Compton scattering of the radio-emitting electrons off the cosmic microwave background (CMB) photons, a reasonable expectation in light of the conclusion reached by Abdo et al. (2010b). We use the simple scaling F_C ≈ F_Sρ₀/ρ_B, where F_C is the total Compton flux, F_S is the radio flux, and ρ₀ ≃ 4 × 10⁻¹³(1 + z)⁴ erg cm⁻³ and ρ_B = B²/8π are the CMB energy density at redshift z and magnetic field energy density, respectively. The expected gamma-ray fluxes of the two radio sources (assuming z = 1) are only at the level of ∼10⁻¹⁴ erg cm⁻² s⁻¹ for B ∼ 1 μG, which is about two orders of magnitude lower than the observed flux. The expected gamma-ray flux would be more consistent with the measured one if these two radio sources are blazars. However, the latter source has a mid-infrared colors (Cutri et al. 2013) not compatible with those of blazars (D'Abrusco et al. 2014). It also has an X-ray counterpart (Binder et al. 2015) with a flux of ∼3 × 10⁻¹⁶ erg cm⁻² s⁻¹ in 0.35–8 keV, which suggests that it might be an X-ray binary in NGC 2403 with an X-ray luminosity ∼4 × 10³⁵ erg s⁻¹. The measured gamma-ray luminosity ∼10³⁹ erg s⁻¹ would then be too high to come from this X-ray source, since the GeV luminosities of X-ray binaries detected in our Galaxy are ≲10³⁶ erg s⁻¹ (Dubus 2013). We cannot rule out completely the possibility of the former source being a blazar. However, the probability is quite small because the amount of blazars satisfying the measured gamma-ray spectral feature is small.⁸ Therefore, it is unlikely that the gamma-ray emission comes from the two radio sources.

On the other hand, we find a recent supernova explosion, SN 2004dj, occurred in the error region of the gamma-ray source. SN 2004dj is peculiar in that it is the second nearest supernova to Earth (after SN 1987A). The probability of chance coincidence between the gamma-ray source and SN 2004dj is estimated to be 0.0022 (see Appendix A for more details). Combining with the transient nature of the gamma-ray source, we suggest that it is physically associated with SN 2004dj.

The GeV emission may originate from the dissipation of the kinetic energy of the SN ejecta, either internally (which may be analogous to the case of a nova; e.g., Li et al. 2017) or externally. In the latter scenario, which is more common for a supernova, the ejecta may encounter a dense gas shell (or clump) in the vicinity of the SN progenitor and drive a strong shock. The gas shell may be formed by the stellar wind of the progenitor star (Dwek 1985). It is not clear at which epoch the shock is driven, but apparently it could be as early as a time period shortly after the SN explosion, so we assume that the shock has been launched somewhere around the beginning of the Fermi mission, i.e., at t_LAT,0 = 1465.6 days (approximately 4 yr), with respect to the SN explosion. Assuming a mass of M_ej = 10 M_⊙ for the ejecta and a total kinetic energy of E_k = 10⁵¹ erg, the bulk of the ejecta moves with a velocity of ${v}_{\mathrm{ej}}=\sqrt{2{E}_{k}/{M}_{\mathrm{ej}}}=3000\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and leads to a shock radius of at least ${R}_{\mathrm{sh}}={v}_{\mathrm{ej}}{t}_{\mathrm{LAT},0}\simeq 4\times {10}^{16}\,$ cm. The distance of the interior shell to the SN progenitor is not supposed to be larger than R_sh. Considering that generally ≲10% of the dissipated kinetic energy of the SN ejecta could be converted to nonthermal particles (Aharonian et al. 2004), the energy in nonthermal particles can reach as high as ∼10⁵⁰ erg. Given that the total gamma-ray energy detected by Fermi-LAT is ∼6 × 10⁴⁷ erg,⁹ we find that ≳0.6% of the ejecta's kinetic energy is needed to be dissipated. As the shock energy is proportional to the ratio between the swept-up material mass and the SN ejecta mass, this implies that a total amount of material of ≥0.06 M_⊙ is swept up by the shock. As the shock sweeps through the shell, it accelerates protons and electrons therein to relativistic energies. In principle, both electrons and protons can radiate gamma-rays, through inverse Compton (IC) scatterings off ambient photons or the decay of π⁰ produced in pp collisions, respectively.

In the leptonic scenario, the target photon field for the IC radiation of electrons is a ∼500 K thermal radiation with a total luminosity of L_IR ∼ 4 × 10³⁸ erg s⁻¹ from the dust formed in the SN ejecta (Meikle et al. 2011; Szalai et al. 2011; Fox et al. 2013), inferred from the observation of Spitzer overlapping with the first one and a half year of Fermi-LAT's observation (i.e., until 2010 May 23). The photon energy density at the shock downstream is given by ${u}_{\mathrm{IR}}\simeq {10}^{-7}{({R}_{\mathrm{sh}}/{10}^{17}\mathrm{cm})}^{-2}\,\mathrm{erg}\,{\mathrm{cm}}^{-3}$, yielding an IC cooling timescale of ${t}_{\mathrm{IC}}\simeq 100{({E}_{e}/50\mathrm{GeV})}^{-1}{({R}_{\mathrm{sh}}/{10}^{17}\mathrm{cm})}^{2}\,$ yr. In the meantime, the accelerated electrons will also radiate in the downstream magnetic field of the shock via the synchrotron process. The ratio between the IC emissivity and the synchrotron emissivity of an electron is equal to the ratio between the energy density of the target photon field and that of the magnetic field. On the other hand, the synchrotron radiation flux should not be larger than 0.227 mJy at 1.4 GHz as observed by the Very Large Array (VLA) after 1578 days of the explosion (Nayana et al. 2018). Therefore, as the gamma-ray flux is ascribed to the IC radiation of accelerated electrons, it constrains the downstream magnetic field to be B < 12 (R_sh/10¹⁷ cm)^−4/3μG (see Appendix B for details). Given that R_sh ≃ 4 × 10¹⁶ cm at the operation starting time of Fermi-LAT, we obtain a more conservative upper limit for the downstream magnetic field B < 41 μG, noting that it is supposed to decrease with time. This value is probably too low considering that the inferred magnetic fields of some well-studied SNRs, such as Cassiopeia A, Tycho, Kepler, SN 1006, and G347.3-0.5, at a much later evolution stage (i.e., from several hundred to several thousand years after the SN explosion) is ≳0.1 mG (Reynolds et al. 2012, and reference therein), and therefore the leptonic interpretation is not favored.

On the other hand, the energy-loss timescale of a proton via the pp collision is t_pp = 7 × 10⁷(n/1 cm⁻³)⁻¹ yr, where n is the target gas density. About one-third of the lost energies will go into pionic gamma-rays. We calculate the pionic gamma-ray spectrum and the secondary electron/positron pair spectrum from the pp collision with the parameterized formulae given by Kamae et al. (2006). The result is shown in Figure 4. The total energy of the accelerated protons is assumed to be W_p = 2 × 10⁴⁸ erg with a differential spectrum in a power-law distribution with slope −2 in the range of 1 GeV–1 TeV. The atom density of the gas shell is set to be n = 1.2 × 10⁷ cm⁻³, so that protons will cool in 5.7 yr. With this assumption, the nondetection of gamma-ray emission in the second 5.7 yr (2014 March 25 (UT)–2019 November 14 (UT)) could be explained as that the proton acceleration has stopped provided the shell's thickness Δ < v_ej × 5.7 yr = 5.4 × 10¹⁶ cm, while protons accelerated earlier have already depleted their energies in pp collisions. Assuming that the shell is located at R = 4 × 10¹⁶ cm and that the shell's thickness is 0.1 of the radius, the total mass of the shell is found to be 0.8M_⊙ with the proposed density. This is consistent with the mass loss of a star in the red supergiant phase of its evolution (Smith et al. 2009; Beasor & Davies 2018). In order to suppress the synchrotron radiation of the secondary pairs to the level consistent with the VLA measurement, we find the downstream magnetic field needs to be weaker than 1 mG.

Zoom In Zoom Out Reset image size

To summarize, using 11.4 yr of Fermi-LAT observations, we searched for possible gamma-ray emission from nearby galaxies in the IRAS Revised Bright Galaxies Sample. In additional to the new detection of gamma-ray emission from two nearby galaxies, Arp 299 and M33, we find a gamma-ray source in the region of the galaxy NGC 2403. The gamma-ray source cannot be explained as pionic gamma-ray emission from the CR–ISM interaction in NGC 2403, as it appears as an outlier of the ${L}_{0.1-100\mathrm{GeV}}\mbox{--}{L}_{8-1000\mu {\rm{m}}}$ correlation for star-forming galaxies and has a fading flux with time. On the other hand, we find that SN 2004dj, the second nearest supernova, is spatially and temporally consistent with the fading gamma-ray source. The probability of the chance coincidence of SN 2004dj with an unrelated background gamma-ray source is estimated to only 0.0022 based on the 4FGL catalog. We interpret this gamma-ray emission arising from the supernova ejecta interacting with a surrounding high-density shell, which converts ∼1% of the ejecta's kinetic energy to relativistic protons. This shows that supernovae can accelerate cosmic-ray protons at a very early stage of their lives.

We would like to thank the anonymous referee for the detailed and constructive report that improved the quality of this paper. The work is supported by the National Key R&D program of China under the grant 2018YFA0404203 and the NSFC grants 11625312 and 11851304.

Using a Poisson distribution, the chance probability of observing a background source in the position of SN 2004dj is ${P}_{\mathrm{ch}}=1-\exp [-\pi ({R}_{0}^{2}+4{\sigma }_{\gamma }^{2})\sum (\gt {F}_{\mathrm{th}})]$, where $\sum (\gt {F}_{\mathrm{th}})$ is the surface density of Fermi-LAT sources with fluxes higher than F_th, σ_γ is the 68% position uncertainties of Fermi-LAT counterpart, and R₀ is the angular distance between the gamma-ray location and the target source SN 2004dj. Since the density distribution of the 4FGL sources in the box defined by $| \sin (b)| \gt 0.25$ is uniform with respect to the angle, we estimate a number density of $\sum (\gt {F}_{\mathrm{th}})=579.4{\mathrm{sr}}^{-1}$ above the flux F_th = 1.18 × 10⁻¹² erg cm⁻² s⁻¹, using the power-law fitting of the cumulative numbers of sources as a function of the threshold fluxes for the 4FGL sources, as shown in Figure A1. Then the chance coincidence probability is estimated to be 0.0022 for R₀ = 0041 and σ_γ = 0024.

Zoom In Zoom Out Reset image size

Leptonic scenario. The average gamma-ray luminosity in 0.1–100 GeV is L_{γ, 0.1–100 GeV} = 3.3 × 10³⁹ erg s⁻¹. Approximating the photon index to be −2, the differential luminosity (i.e., ${E}_{\gamma }^{2}{dN}/{{dE}}_{\gamma }{dt}$) is ${L}_{\gamma }({E}_{\gamma })\,=\,{L}_{\gamma ,0.1-100\mathrm{GeV}}/\mathrm{log}(1000)\,=4.8\times {10}^{38}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, independent of the photon energy. As we ascribe the gamma-ray emission to the IC radiation of accelerated electrons, L_γ(E_γ) should be equal to the total power of electrons at energy E_e, where

$$\begin{eqnarray}&&{E}_{e}=50{({T}_{\mathrm{IR}}/500{\rm{K}})}^{-1/2}{\left({E}_{\gamma }/1\mathrm{GeV}\right)}^{1/2}\,\mathrm{GeV}.\end{eqnarray} \tag{ B1 }$$

The differential energy spectrum of electrons can be found by ${E}_{e}^{2}{dN}/{{dE}}_{e}\,=\,{L}_{\gamma }({E}_{\gamma }){t}_{\mathrm{IC}}({E}_{e})$ $=\,1.5\,\times {10}^{48}{({E}_{e}/50\mathrm{GeV})}^{-1}{({R}_{\mathrm{sh}}/{10}^{17}\mathrm{cm})}^{2}\,$ erg. According to Equation (B1), we find the energy of electrons responsible for 0.1 GeV gamma-rays to be 16 GeV. To calculate the synchrotron luminosity in the radio band, we need to extrapolate the electron spectrum to energies below 16 GeV. Note that to explain the photon index of −2, the electron spectral index above 16 GeV needs to be −3. We assume a break at 16 GeV and a flat spectrum (${dN}/{{dE}}_{e}\propto {E}_{e}^{-2}$) below 16 GeV to reduce the predicted radio flux. Such a flat spectrum is also consistent with the prediction of the canonical diffusive shock acceleration theory for strong shocks (Bell 1978). As such, we have ${E}_{e}^{2}{dN}/{{dE}}_{e}=4.7\times {10}^{48}{({R}_{\mathrm{sh}}/{10}^{17}\mathrm{cm})}^{2}$ for E_e < 16 GeV. The synchrotron cooling timescale of electrons radiating at a characteristic frequency ν is

$$\begin{eqnarray}&&{t}_{\mathrm{syn}}=4.2\times {10}^{7}\,{\left(\displaystyle \frac{\nu }{1.4\mathrm{GHz}}\right)}^{-1/2}{B}^{-3/2}\,{\rm{s}}.\end{eqnarray} \tag{ B2 }$$

Then we can estimate the synchrotron luminosity by

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{syn}}(\nu ) & = & \displaystyle \frac{{E}_{e}^{2}{dN}/{{dE}}_{e}}{{t}_{\mathrm{syn}}}\\ & = & 1.1\times {10}^{41}\,{\left(\displaystyle \frac{{R}_{\mathrm{sh}}}{{10}^{17}\mathrm{cm}}\right)}^{2}{\left(\displaystyle \frac{\nu }{1.4\mathrm{GHz}}\right)}^{1/2}{B}^{3/2}\,\mathrm{erg}\,{{\rm{s}}}^{-1}.\end{array}\end{eqnarray} \tag{ B3 }$$

The 1.4 GHz luminosity observed by VLA is 4.6 × 10³³ erg s⁻¹ at 1578 days after the explosion (Nayana et al. 2018). Since there could be other processes contributing to the flux, L_syn should be smaller than this value, and this results in B < 12(R_sh/10¹⁷ cm)^−4/3 μG. Such a strong limit on the magnetic field might in principle be relaxed, if one argues the possible existence of an extra intense radiation field around the SN during the period of Fermi-LAT's detection, which in the meantime is hidden to observers, or if one argues a low-energy cutoff present in the electron spectrum above the energy responsible for radio emission. These two possibilities are, however, not supported either observationally or theoretically, so we do not discuss them further in this work.

Hadronic scenario. In pp collisions, approximately one-third of produced pions are neutral pions and the other two-thirds are charged pions that will decay into neutrinos and electrons/positrons (hereafter we do not distinguish positrons from electrons for simplicity) with a ratio of 3:1. On average, the energy of a secondary electron is half of the energy of a gamma-ray photon produced by protons of the same energy, while the secondary electron energy production rate is also half of that of gamma-rays, i.e.,

$$\begin{eqnarray}&&{E}_{e}^{2}\displaystyle \frac{{{dN}}_{2\mathrm{nd}}}{{{dE}}_{e}{dt}}=\displaystyle \frac{1}{2}{L}_{\gamma }({E}_{\gamma })=2.4\times {10}^{38}\,\mathrm{erg}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ B4 }$$

for 50 MeV < E_e < 50 GeV. The radio-emitting electrons cannot be cooled during Fermi-LAT's detection period, so the accumulated secondary electron spectrum can be given by ${E}_{e}^{2}({{dN}}_{2\mathrm{nd}}/{{dE}}_{e})={E}_{e}^{2}({{dN}}_{2\mathrm{nd}}/{{dE}}_{e}{dt}){\tau }_{\mathrm{GeV}}$, with τ_GeV being the duration of the GeV emission. The duration should last at least 5.7 yr as Fermi-LAT detected, while it can be at longest 9.7 yr if the GeV emission starts shortly after the SN explosion. Similar to Equation (B3), we estimate the synchrotron luminosity of secondary electrons by

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{syn},2\mathrm{nd}}(\nu ) & = & \displaystyle \frac{{E}_{e}^{2}{{dN}}_{2\mathrm{nd}}/{{dE}}_{e}}{{t}_{\mathrm{syn}}}\\ & = & 1.7\times {10}^{39}\left(\displaystyle \frac{{\tau }_{\mathrm{GeV}}}{9.7\,\mathrm{yr}}\right){\left(\displaystyle \frac{\nu }{1.4\mathrm{GHz}}\right)}^{1/2}{B}^{3/2}\,\mathrm{erg}\,{{\rm{s}}}^{-1},\end{array}\end{eqnarray} \tag{ B5 }$$

and we get B < 1.9(τ_GeV/9.7 yr)^−2/3 mG by requiring Equation (B5) be smaller than the measured 1.4 GHz luminosity 4.6 × 10³³ erg s⁻¹. This result is consistent with the one obtained in Section 4 (see also Figure 4) by a more accurate numerical calculation.