© The Authors 2025

1. Introduction

The central starburst (CSB) regions in intensely star forming galaxies are optimal for the study of relativistic (“cosmic ray”; CR) protons (CRp) and electrons (CRe) that are accelerated by supernova (SN) – driven shocks, and effectively couple to the dense matter, radiation, and magnetic fields. Measurements of the spectral energy distribution (SED) of the nonthermal (NT) radiative yields of CR interactions in CSB provide essential insights into the particle spectra and the physical conditions in star-forming environments.

Gas and CR densities and magnetic and radiation field distributions vary significantly across galaxy disks, generally requiring a spectro-spatial treatment based on a solution to the diffusion-advection equation in cases in which the spatial profiles of relevant key quantities may be specified (e.g., Rephaeli & Sadeh 2019, 2024; Heesen 2021; Strong et al. 2007, 2010; Taillet & Maurin 2003; Wallace 1980; Jones 1978; Syrovatskii 1959). However, in order to gain an insight into spatially averaged NT properties in the relatively small CSB region, we started from the CR yields and deduced the underlying CR steady-state spectra, making it possible (in principle) to deduce the injection spectra.

In previous works (Persic & Rephaeli 2022; Persic et al. 2024) we examined NT emission in nearby galaxies with low-to-moderate star formation rates – the Magellanic Clouds, M31 and M33. Adopting a one-zone model for the extended disk emission, we modeled disk radio emission as a combination of synchrotron and thermal bremsstrahlung and γ-ray measurements as neutral-pion (π⁰) decay (following CRp interactions with the ambient gas) and leptonic emission.

In this follow-up study, we consider the local starburst (SB) galaxy M82. The enhanced star formation occurs in the CSB region, which we model as a spheroid of semimajor (semiminor) axis R_s = 338 (h_s = 165) pc, where the NT emission relevant to this study (see below) is produced. Located at a distance of D = 3.6 Mpc, M82 is seen almost edge-on, with an inclination angle (to the line of sight) ofi ∼ 80°.

A member of the M81 group, M82 shows traces of past tidal encounters. One encounter with M81 less than 400 Myr ago (Hutton et al. 2014), shown by an intergalactic gas bridge, triggered a major SB at a rate of ∼10 M_⊙ for ∼50 Myr. Two subsequent SBs followed, the last of which took place ∼1 Myr ago (Barker et al. 2008). A prominent galactic wind (“superwind”) shoots out of both sides of the galactic disk in nearly perpendicular directions, out to several kiloparsecs, likely driven by the superposed winds of many massive stars and SN explosions. The CSB region is contained in the innermost (R ≲ 300 pc) disk, surrounded by a molecular ring that has a radius of 400 pc.

M82 emission levels are in accord with the well-known correlations between radio, infrared, and γ-ray emission for star-forming galaxies (radio–IR: Condon 1992, Bell 2003, Persic & Rephaeli 2007; γ–IR: Ackermann et al. 2012, Kornecki et al. 2020). [Infrared (IR) emission is variously estimated by the total-IR (TIR) band (8−1000 μm), the far-IR (FIR) band (40−120 μm), or a specific monochromatic flux (e.g., 60 μm).] These correlations are traced back to CR acceleration by SN explosions, and the interaction of particles and radiation with the dusty magnetized medium characterizing the sites where massive progenitor stars were formed (Lacki et al. 2010; Bell 2003).

This study aims to determine NT aspects (CR spectra, magnetic field) in the CSB region by spectrally modeling emission in the relevant energy ranges accessible to observations. Radio and γ-ray spectra from the central region of M82 have long been available (radio: Klein et al. 1988, Seaquist & Odegard 1991, Carlstrom & Kronberg 1991, Williams & Bower 2010, Peel et al. 2011; γ ray: VERITAS Collaboration 2009; Abdo et al. 2010). Several studies of the NT emission from the CSB were published, both before and after its γ-ray detection. Among the former, early estimates (Völk et al. 1989, 1996; Akyüz et al. 1991) and subsequent more detailed models (Persic et al. 2008; De Cea del Pozo et al. 2009) predicted that γ-ray emission from M82 was detectable by then-current and upcoming detectors. Among the latter, Lacki et al. (2011), Peretti et al. (2019), Krumholz et al. (2020), and Ambrosone et al. (2022, 2024) discussed calorimetric properties and CR physics in the dense CSB environment.

The motivation for the present study stems from an important recent observational development. Diffuse NT X-ray emission was detected in an ellipse of semimajor and semiminor axis $19\stackrel{″}{.}4\times 9\stackrel{″}{.}4$ (338 × 165 pc) morphologically consistent with the radio and far-infrared (FIR) emission maps of M82’s CSB region (Iwasawa et al. 2023). The co-spatial emission in the three bands and the similarity of the NT X-ray and radio-synchrotron spectral indices (α_X ≃ α_r ≃ 0.7) led Iwasawa et al. (2023) to interpret the detected NT X-ray emission as being due to Compton scattering of the FIR photons by radio-emitting CRe. The direct relevance of Compton scattering is obvious given the co-spatial presence in the CSB of CRe and FIR photons. The detected Compton/FIR flux enables a calibration of the CRe spectrum to be made given the obvious relation between Compton emissivity and CRe number density, j_Compt ∝ n_e. Once the CRe spectrum is fully determined, the mean (volume-averaged) magnetic field can be directly deduced. This procedure, which does not rest on equipartition arguments, is possible for the first time in this prototypical SB galaxy.

The paper is organized as follows. In Sect. 2 we review archival data on broadband NT interstellar emission from the CSB of M82, and we summarize our analysis of 16.3 yr of data collected by the Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope (fully described in Appendix A). In Sect. 3 we review the radiation fields permeating the CSB region. In Sects. 4 and 5 we describe and discuss our SED modeling. Our conclusions are summarized in Sect. 6.

2. Observations of NT emission from the CSB

Diffuse NT emission from M82 has been measured in radio/microwave, X-ray, and γ-ray (e.g., Gendre et al. 2013; Iwasawa et al. 2023; Ohm 2016), in addition to extensive measurements of diffuse (thermal) emission in the near-ultraviolet, optical, and IR (Hutton et al. 2014; Gurzadyan et al. 2015). In this section we review some of the published observations of NT CSB emission and describe our newly reduced 16.3 years of accumulated Fermi-LAT data.

2.1. Radio

M82 has been extensively observed over the years (e.g., Kellermann et al. 1969, Klein et al. 1988, Carlstrom & Kronberg 1991, Williams & Bower 2010, Peel et al. 2011). In the present study we focus on emission from the small CSB region from which NT X-ray and γ-ray emission has been detected. The spectrum of radio emission from this region was presented by Adebahr et al. (2013), based on combined archival Very Large Array (VLA) data (4.996 GHz, 8.327 GHz) and Westerbork Synthesis Radio Telescope (WSRT) data (0.326 GHz and 1.364 GHz). The VLA data (Seaquist & Odegard 1991; Reuter et al. 1992) had previously shown emission mostly confined to the CSB region of M82.

Improved data reduction and instrumental techniques enabled Adebahr et al. (2013) to analyze the radio maps in increased detail. The archival WRST data were reduced by Adebahr et al. (2013) using a new calibration technique able to reach the high dynamical ranges needed to resolve the complex source morphology and map the residual diffuse emission in the central region. The combined datasets yielded total power maps at 0.326, 1.364, 4.996, and 8.327 GHz. At each frequency the intensity was integrated over the angular size of the CSB region (assumed by Adebahr et al. (2013) to have a radius of 450 pc, at a distance of 3.5 Mpc: i.e., $26\stackrel{″}{.}5$). A fifth measurement at 1.67 GHz (by Braun et al. 2007) was included in their power-law (PL) spectral fit by means of free-free absorption by a screen of ambient ionized gas, S = S₀ (ν/ν₀)^−α e^−τ_ff, with α = 0.62 ± 0.01 and τ_ff = 8.2 × 10⁻²ν^−2.1EM/T_e^1.35, where EM = (3.16 ± 0.10)×10⁵ cm⁻⁶ pc is the emission measure and T_e ∼ 10⁴ K is the warm gas temperature. We adopted the spectral data of Adebahr et al. (2013) after rescaling each flux value by the appropriate factor to match the assumed size of the CSB region from which the NT X-ray emission was measured by Iwasawa et al. (2023). The rescaled data are listed in Table 1.

2.2. X-ray

X-ray measurements of M82 have shown early on that most of the emission comes from the central region (e.g., Watson et al. 1984), possibly due to either a variable accreting source (Ptak & Griffiths 1999) or a hot (kT ∼ 6−9 keV) thermal interstellar plasma with a solar metal abundance of ∼0.3 (Cappi et al. 1999). The presence of this hot ejecta of many SN explosions, whose superposed shocks drive a superwind into the intergalactic medium, was suggested by Griffiths et al. (2000). Later Strickland & Heckman (2007) found the CSB region to be filled by 6.7 keV He α line emission as well as a marginally significant 6.4 keV Fe Kα line and possibly NT continuum.

An improved insight into the origin of X-ray emission came recently from analysis of 570 ks of Chandra data by Iwasawa et al. (2023), who found evidence of a NT emission component, at a level of ∼70% of the 4−8 keV emission. They deduced that this emission is morphologically consistent with the radio map of the CSB, which they interpreted as arising from Compton scattering of local (i.e., SB-sourced) FIR photons by radio-emitting CRe. This result can be used to estimate the spectral X-ray luminosity at 5 keV: for an energy spectral index of α = 0.7 (Iwasawa et al. 2023), the photon spectrum is N_ph(ϵ) = N_1 keV (ϵ/keV)^−1.7 ph/(cm² s keV), with energy in kilo-electronvolts. The 4−8 keV spectral energy flux is ϕ_{4 − 8 keV} = ∫₄⁸ϵ N_ph(ϵ) dϵ = 1.17 N_1 keV keV/(cm² s). Since this should equal f_{4 − 8 keV} = L_{4 − 8 keV}/(4π D²) erg cm⁻² s⁻¹, where L_{4 − 8 keV} = 0.75 × 10³⁹ erg s⁻¹ (Iwasawa et al. 2023), it follows that N_1 keV = 2.21 × 10⁻⁴ and N_5 keV = N_1 keV(5/keV)^−1.7 ≃ 1.43 × 10⁻⁵. The monochromatic 5 keV energy density flux, f_5 keV = N_5 keVE_5 keV, is then (Table 2)

$$\begin{array}{c}\hfill f_{5\mathrm{keV}}=(1.14\pm 0.08)\times {10}^{-13}\mathrm{erg}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}.\end{array}$$(1)

2.3. γ ray

M82 was measured to be a point source with (non-varying) > 200 MeV emission detected at a 6.8σ significance based on 11 months of scanning-mode data since the inception of the Fermi mission (Abdo et al. 2010). The measured emission was interpreted as originating from the interaction of relativistic and ambient protons, suggesting a link between massive star formation and γ-ray emission. This interpretation was found to be valid for star-forming galaxies in general (Ackermann et al. 2012; also Mannheim et al. 2012 and Chen et al. 2024). Here we analyze 195.6 months (16.3 yr) of scanning-mode LAT data (Table 3), detailed in Appendix A.

Based on 137 h of high-quality Very Energetic Radiation Imaging Telescope Array System (VERITAS) data, collected in 2008–2009, a steady > 0.7 TeV flux of (3.7 ± 0.8_stat ± 0.7_syst)×10⁻¹³ cm⁻² s⁻¹ was detected (at a 4.8σ significance level) from M82 (VERITAS Collaboration 2009) – the first detection of a SB galaxy at tera-electronvolt energies. The differential photon γ-ray spectral index was Γ = 2.5 ± 0.6_stat ± 0.2_syst. Owing to a follow-up observational campaign in 2011−2022, 254 hours of good-quality data are now available. The expanded dataset and improved analysis techniques yielded a significant (6.5σ) point-like detection (VERITAS Collaboration 2025). The resulting flux is F(> 450 GeV) = (3.2 ± 0.6_stat ± 0.6_syst)×10⁻¹³ cm⁻² s⁻¹ (∼0.4% of the Crab Nebula flux above the same threshold), with a spectral slope of Γ = 2.3 ± 0.3_stat ± 0.2_syst, and no obvious spectral cutoff at energies ≤5 TeV. The 2008−2022 VERITAS data are summarized in Table 4 and Appendix B. In summary, the broadband γ-ray SED of M82 used in the present study extends from 50 MeV to 16.25 TeV, with data taken during the periods 2008−2022 (VERITAS) and 2008−2024 (Fermi-LAT).

3. Radiation fields

A detailed description of the ambient radiation fields is needed to predict γ-ray emission from Compton scattering of target photons by the radio-emitting CRe. The ambient radiation fields are extragalactic and local. The former include the cosmic microwave background (CMB), a Planckian with a temperature of T_CMB = 2.735 K and energy density of u_CMB = 0.25 eV cm⁻³, and the extragalactic background light (EBL). The EBL is the superposed emission from direct (optical) and dust-reprocessed (into IR) starlight integrated over the star formation history of the Universe. Its spectrum has two humps, optical and IR, peaking at ∼1 μm and ∼100 μm. These humps are described as diluted Planckians, characterized by a temperature, T, and a dilution factor, C_dil¹. We adopted the commonly used EBL model by Franceschini & Rodighiero (2017), numerically approximated as a combination of diluted Planckians (Eq. (1) of Persic et al. 2024).

The dominant radiation field in the CSB region arises from the local stellar population. Similar to the EBL in shape and origin, it is dominated by the IR and optical humps. Full-band luminosities of these thermal components are needed to determine n(ϵ), the spectral distributions of target IR/optical photons (of energy ϵ) that are Compton-upscattered to X-ray/γ-ray energies. We proceeded as follows:

(i) The total optical-band luminosity was computed from the narrow B-band luminosity (B = 8.76, B−V = 1.15; Blackman et al. 1979) by applying a suitable bolometric correction, B_bol = B + BC_B with BC_B = −0.85 − (B − V) (Buzzoni et al. 2006; Persic & Rephaeli 2019). The corresponding number density spectrum of optical photons is a (diluted, unmodified) Planckian. The optical energy density inside the source (assumed to be homogeneous), calculated from u_opt ≃ (9/4)L_IR;  opt/[4πR_s²c] (Ghisellini 2013)², is u_opt ∼ 6.8 × 10⁻¹⁰ eV cm⁻³. The optical-hump temperature, estimated applying Wien’s displacement law to the optical-peak frequency, 1.5 μm (Vaccari & Franceschini 2008), is $T_{\mathrm{opt}}=2000{\mathrm{K}}^{[c]}$. The corresponding dilution factor is log(C_dil^opt) =  − 8.248.

(ii) The total-infrared (TIR) energy density was derived from u_FIR ∼ 8.2 × 10⁻¹⁰ eV cm⁻³ (Iwasawa et al. 2023), multiplied by 1.65 to convert FIR to TIR fluxes (Persic & Rephaeli 2007), and further corrected for the mean photon-escape time in a (quasi) spherical homogeneous source (see above). This yielded u_TIR ∼ 3 × 10⁻⁹ erg cm⁻³. The calculation of the IR dilution factor, $C_{\mathrm{dil}}^{\mathrm{IR}}$, was less straightforward than in the optical case because the ambient IR emission is clearly not blackbody but “graybody.” The calculation is outlined below.

IR graybody emission

The spectral flux from a dusty cloud is described as a modified blackbody (graybody), i.e., the product of a Planckian function, B(ν, T_d), by the dust emissivity function, Q(ν) = 1 − e^−τ_gb, where τ_gb = (ν/ν₀)^β is the dust-related optical depth (Hildebrand 1983) and 0 < β ≲ 2 is the dust emissivity index (e.g., Yun & Carilli 2002)³. The reference frequency, ν₀, corresponds to τ_gb = 1. At frequencies of ν ≫ ν₀, the dust spectrum is optically thick and reduces to a blackbody, whereas at ν ≪ ν₀ it describes an optically thin graybody, j_ν(β, T_d) = B(ν, T_d) (ν/ν₀)^β ∝ ν^2 + βT_d (which reduces to the Rayleigh-Jeans law for β = 0).

The IR emission of M82 is represented by an optically thin graybody, with ν₀ corresponding to the peak emission frequency, ν₀ = ν_p = 2.5 × 10¹² Hz (Yun & Carilli 2002; Peel et al. 2011). In evaluating Compton scattering of IR photons by CRe in the CSB we assume that the spectral energy density of target photons is n_IR(ϵ) = 8π/(h³c³) ϵ²/(e^ϵ/k_BT_d − 1)(ϵ/ϵ₀)^β (which replaces Eq. (A.7) of Persic et al. 2024). As for the (β, T_d) pair, whose combination in the ν^βB(ν, T_d) (graybody) fit describes the dust spectral energy profile (e.g., Fig. 1 of Yun & Carilli 2002), we used several published values relative to various FIR datasets⁴ that quantify the overall FIR modeling uncertainty. For these (β, T_d) pairs, Eq. (12) of Elia & Pezzuto (2016) specifies values for the peak emission frequency, ν_p ∼ 2.5 × 10¹² Hz, consistent with the observed value (e.g., Peel et al. 2011). In principle these ν^β B(ν, T) fits are mutually consistent; nevertheless, we present SED models for all of them to check the specificity of the NT SED versus the IR SED.

To determine the appropriate dilution factor for an optically thin graybody, we proceeded as follows. The power radiated by an optically thin graybody per unit surface area, $W_{\mathrm{GB}}={\int }_0^{\infty }{\left(\frac{\nu }{{\nu }_0}\right)}^{\beta }{B}_{\nu }(T)\mathrm{d}\nu$, is

$$\begin{array}{c}\hfill W_{\mathrm{GB}}=\frac{2\pi k_B^{4+\beta }}{h^{3+\beta }{c}^2}\frac{1}{{\nu }_0^{\beta }}\mathrm{\Gamma }(4+\beta )\zeta (4+\beta )T^{4+\beta },\end{array}$$(2)

where Γ(z) and ζ(z) are Euler’s gamma function and Riemann’s zeta function, respectively (see Eqs. (30), (31) of Elia & Pezzuto 2016).

The graybody equivalent of the blackbody displacement law is ${\nu }_p\simeq \frac{k_{B}T}{h}(3+\beta )$ (Elia & Pezzuto 2016, Eq. (12)). If ν₀ = ν_p, then

$$\begin{array}{cc}\hfill W_{\mathrm{GB}}& =\frac{2\pi k_B^4}{h^3{c}^2}\frac{1}{{(3+\beta )}^{\beta }}\mathrm{\Gamma }(4+\beta )\zeta (4+\beta )T^4\hfill \\ \hfill & ={\sigma }^{\mathrm{GB}}(\beta )T^4,\hfill \end{array}$$(3)

where we introduced the quantity

$$\begin{array}{c}\hfill {\sigma }^{\mathrm{GB}}(\beta )\equiv \frac{2\pi k_B^4}{h^3{c}^2}\frac{1}{{(3+\beta )}^{\beta }}\mathrm{\Gamma }(4+\beta )\zeta (4+\beta ).\end{array}$$(4)

A graybody reduces to a blackbody for β = 0, i.e.,

$$\begin{array}{c}\hfill {\sigma }^{\mathrm{GB}}(0)={\sigma }_{\mathrm{sb}}=\frac{2\pi k_B^4}{h^3{c}^2}\mathrm{\Gamma }(4)\zeta (4),\end{array}$$(5)

where σ_sb is the Stefan-Boltzmann constant. Thus Eq. (4) can be written as

$$\begin{array}{c}\hfill {\sigma }^{\mathrm{GB}}(\beta )=\frac{1}{{(3+\beta )}^{\beta }}\frac{\mathrm{\Gamma }(4+\beta )\zeta (4+\beta )}{\mathrm{\Gamma }(4)\zeta (4)}{\sigma }_{\mathrm{sb}}.\end{array}$$(6)

Equation (6) yields the substitute to the Stefan-Boltzmann constant for an optically thin GB specified by the spectral index, β, with ν₀ = ν_p. The latter is generally valid in SB galaxies (Yun & Carilli 2002). Equation (6) is a monotonically decreasing function of β; clearly, lim_β → 0σ^GB(β) = σ_sb and lim_{β → ∞}σ^GB(β) = 0. Indeed, σ^GB ≈ 0 already at β = 10. The dilution factor clearly depends on the parameter pair (β, T) describing the graybody, $C_{\mathrm{dil}}^{\mathrm{IR}}=u_{\mathrm{IR}}/\left[\frac{4}{c}{\sigma }^{\mathrm{GB}}(\beta )T^4\right]$ (specified below)⁵.

4. SED modeling

We aim to determine NT quantities (CR, magnetic field) in the CSB region by comparing the SED dataset to the emission predicted by relevant radiative processes (well-known formulae are collected in Appendix A of Persic et al. 2024). Given the relatively small size of the CSB and lacking any spatial information on NT emission from this region, we assumed the CR spectral distributions to be time-independent, locally isotropic, and uniform.

Below are the CR spectral distributions:

• The CRp spectrum is a PL in energy, N_p(E_p) = N_p, 0E_p^−q_p cm⁻³ GeV⁻¹ for m_pc² < E_p < E_p^max (energies in giga-electronvolts). The normalization and slope, N_p, 0, q_p, are free parameters.

• Secondary CRe are produced from π^± decays following p−p interactions of CRp with the ambient gas. The source spectrum of freshly created secondary CRe, essentially a PL in energy, is modified by energy losses: radiative losses (Coulomb, bremsstrahlung, synchrotron/Compton) and, in the present case, also superwind advection losses (Longair 1981),

$$\begin{array}{c}\hfill b_{\mathrm{adv}}=2.16\times {10}^{-17}\gamma \left(\frac{V_{\mathrm{wind}}}{\mathrm{km}/\mathrm{s}}\right){\left(\frac{r_h}{\mathrm{kpc}}\right)}^{-1}{\mathrm{s}}^{-1}.\end{array}$$(7)

With the CRp spectrum and the environmental parameters (density, magnetic field, superwind speed) in the CSB region specified (see Tables 5 and 6), the secondary electron spectrum has no free parameters. For computational ease, the secondary spectrum can be analytically approximated as in Eq. (A.29) of Persic et al. (2024).

• The primary CRe spectrum is a PL in Lorentz factor γ, N_e(γ) = N_e, 0γ^−q_e cm⁻³ (per unit γ) for γ_min < γ < γ_max (with γ_min = 100). The normalization and slope, N_e, 0 and q_e, are free parameters.

Our modeling procedure began with fitting a pionic emission profile to the high-frequency (ν ≳ 10^24.5 Hz) γ-ray data with free normalisation, a slope, and a high-end cutoff; photon–photon opacity is not important, as is discussed in Sect. 6 and shown in Fig. 3. The emission spectrum from π⁰-decay has constraining power owing to its characteristic “shoulder” at ∼100 MeV; if apparent in the data at higher energies, a spectral cutoff corresponds to E_p^max. For the set of models developed in this study, corresponding to different (β, T_d) pairs (see Sect. 5), we deduced CRp spectral slopes of q_p ≃ 2.3, energy cutoffs of E_p^max ≈ 7 TeV, and energy densities of u_p ≃ 385 eV cm⁻³ (Table 6).

Next we used the secondary- and primary-CRe spectra to calculate the monochromatic 5 keV Compton/starlight emission. With the secondary CRe spectrum fully determined, unlike the primary CRe spectrum, which was assumed, fitting the 5 keV point allowed us to determine N_e, 0. We assumed a primary CRe injection index of q_e ≤ 2.2. This assumption was based on observational evidence that the distribution of radio synchrotron indices for Galactic SN remnants (SNRs), i.e., the putative sites of Galactic CR acceleration, peaks at 0.50 ≲ α_r ≲ 0.55 (Klein et al. 2018), i.e. 2 ≲ q_e ≲ 2.1. While there are many SNRs in the CSB region, their filling factor in the CSB is low, so the CRe index averaged over the region must be steeper than that measured locally in SNRs.

Radio emission was modeled as CRe synchrotron in a disordered (spatially averaged) magnetic field, B. Deducing B directly from spectral modeling, rather than assuming particle-field energy equipartition, is a much-needed development in NT spectral analyses of M82 and other SB galaxies. The lack of an obvious high-ν cutoff in the radio spectrum makes it impossible to evaluate γ_max from radio data alone; however, the Compton/starlight γ-ray emission suggests values of γ_max ∼ 5.5 × 10⁴. At high radio frequencies, the ν^−0.1 component represents diffuse thermal free-free emission from a warm (∼10⁴ K) ionized plasma (e.g., Spitzer 1978), possibly concentrated in the CSB (Shopbell & Bland-Hawthorn 1998). This emission may be gauged by the Hα flux, F(Hα), if both come from the same HII regions, since in this case the relevant warm-plasma parameters (temperature, density, filling factor) are the same. In this case the measured F(Hα) can be used to predict the level of free-free emission. To do this we used Eq. (17) of Klein et al. 2018 with a warm-plasma temperature of T_e = 10⁴ K and (aperture-integrated) flux of F(Hα) = (0.79 ± 0.05)×10⁻¹⁰ erg cm⁻² s⁻¹ (Kennicutt et al. 2008). The radio spectral data and model are shown in Fig. 2.

As in our (NT) analyses of M31 and M33 (Persic et al. 2024), we consider thermal and NT electrons to be homogeneously mixed throughout the CSB region (Förster Schreiber et al. 2001; Seaquist et al. 1985). The emerging intensity is I_ν ∝ j_ν(1 − e^−τ_ff(ν))/τ_ff(ν), where j_ν is the (synchrotron plus thermal free-free) spectral emissivity and τ_ff(ν) is the optical depth for free-free absorption by thermal plasma (Persic et al. 2024, Eqs. (2) and (3)). We consider this “mixed absorption-emission” geometry to be more realistic then the “screened source” geometry adopted by Adebahr et al. (2013), who considered a synchrotron source absorbed by a foreground screen of 10⁴ K plasma. Our model matches the observed radio spectrum equally well by assuming EM ≈ 0.75 × 10⁶ cm⁻⁶ pc, which implies that ⟨n_e²⟩^1/2 ≈ 40 cm⁻³, consistent with the derived range of free electron densities⁶.

Having determined the CRe spectra, the full NT bremsstrahlung and Compton-scattered starlight γ-ray yields were calculated using standard emissivity formulae. The total radiative yields from CRe interactions were added to the pionic yield to fit the full Fermi-LAT range of γ-ray data (Fig. 1). An acceptable fit was obtained iteratively by varying the parametrized proton spectrum. We find that γ-ray emission is dominated by the pionic component over the high-frequency LAT and VERITAS frequency range, but Compton-scattered starlight dominates emission at the lowest LAT energies (50−100 MeV).

Fig. 1. X-ray/γ-ray range of the CSB SED. The emission model (thick solid curve), including total (i.e., primary plus secondary) Compton/M82-EBL-CMB radiation (dashed curve) and NT bremsstrahlung (dash-dotted curve), and a pionic component (dotted curve), is overlaid with data from Tables 2–4 (black dots). The leptonic component is dominated by the Comptonization of the local FIR radiation. The 5 keV flux, resulting from the NT 4−8 keV flux measured by Iwasawa et al. (2023), is fit by (basically) the local Compton/FIR yield, in which the secondary electron spectrum is uniquely set by the CRp spectrum fit to the γ-ray data and the primary-electron spectrum is normalized to match the 5 keV flux.

5. Results and discussion

Different FIR datasets analyzed by various authors have resulted in a range of (β, T_d) values. This observational uncertainty in the parameters of the FIR graybody model for the dusty CSB region results in a range of viable SED models presented in the previous section. We note that the models listed in Table 6 cover the observational range of values of the parameters β and T_d; clearly, their respective statistical likelihoods (quantified by the value of χ_ν²) do not imply a relative (model) preference.

Our SED models share basic similarities but also differ significantly, as is briefly discussed below. The models are similar in that the γ-ray spectrum is viably modeled as pionic emission for most of the measured range except at the lowest LAT energies, where the peak of the Compton-scattered intense local FIR radiation field dominates. The monochromatic 5 keV flux lies on the rising branch of such Compton/FIR starlight. Thus the CRe and CRp spectra are anchored to the X-ray and high-energy γ-ray emissions, respectively. The radio spectrum is synchrotron emission by primary and secondary CRe in comparable contributions (plus a very minor level of thermal bremsstrahlung), so the volume-averaged magnetic field is directly linked to the radio spectrum.

Model differences stem from the fact that, for a given T_d, different values of β imply different densities of the FIR target photons that are Compton-scattered by the radio electrons into the X-ray/γ-ray domain. This is exemplified by the models having β = 1, 1.3, 1.5: a higher β implies a higher Compton/(CMB + starlight) emission, so for consistency with the 5 keV flux constraint, the associated primary CRe spectrum must be renormalized lower. Consequently, the deduced B has to be higher in order for the computed synchrotron emission to match the radio data (Table 6).

Our treatment is based on the attainment of a steady state in the CSB region; namely, that particle acceleration and energy loss rates are roughly comparable during the SB phase. Supporting evidence for the validity of this assumption can be seen from the (deduced) near-equipartition between CR and magnetic energy densities, by the appreciable role of secondary CRe in providing the requisite complementary radio synchrotron and X-ray/γ-ray Compton/FIR emission, and – more generally – by the effectiveness of an isotropic, time-independent modeling of the SED data.

Clearly, this study focuses on the radiative manifestations of NT degrees of freedom in the dense magnetized, likely turbulent CSB gas. As such, it is obviously quite distinct from the more micro-physically detailed (but necessarily parametrized) MHD studies of the strongly coupled media in the CSB. Insights gained from the latter could possibly have some relevance for CR propagation modes and the small-scale morphology of magnetic fields. Given that our analysis is merely spectral (not spectro-spatial), it only relates to (what are essentially) the coarse-grained quantities that characterize the gas, NT particles, and magnetic field. Insights from the results of our spectral analysis could be useful in limiting the parameter space of the more comprehensive MHD modeling of turbulent gas in gas-rich environments, such as in CSB regions, and throughout galactic disks (e.g. Poggianti et al. 1999; Bland-Hawthorn et al. 2024, 2025). In the following we discuss key aspects of our analysis and its results.

• Lepto-hadronic interplay. Although the γ-ray spectrum is naturally interpreted as a pionic emission at log(ν)≳25, a leptonic component (i.e., Compton-scattered starlight⁷) is necessary to model LAT data at lower frequencies. The interplay between leptonic and hadronic emissions in the model is best revealed at log(ν)∼22.7, where the two components crossover: Compton-scattered starlight declines steeply from its peak at log(ν)∼22 and pionic emission rises steeply to its own peak at log(ν)∼23.25 (Fig. 1). At about the same frequency, log(ν)∼23.5, the Compton/optical peak enhances the hadronic peak. As a consequence the model SED exhibits a local minimum between the leptonic and the hadronic peaks.

The LAT good-quality (spectrally resolved) data in the range log(ν) = 22.7−24 of the pionic hump enable the two different spectral components to be separated. However, at lower frequencies (i.e. ≲200 MeV) just below this range the data points are upper and lower limits, and thus can only bracket the viable range of the model SED. With the predicted (fully determined) Compton/FIR shape of the SED at energies just below ∼207 MeV (log(ν) = 22.7), it seems likely that the last two upper limits must closely trace the SED in order for the latter to hit the first statistically sound flux point (where the upper and lower limits converge), at log(ν) = 22.87. Should the curve be even slightly lower, it would violate the quickly rising lower limits. Thus this highly structured SED curve and the arrowhead pattern of ≲200 MeV limits converging to the first statistically detected flux point enable the (model-based) conclusion that these limits meaningfully constrain the SED model, so much so that the LAT-data SED can be safely traced down to log(ν) = 22.4 (i.e. ∼100 MeV). Higher-sensitivity measurements are needed to better constrain the emission in the crucial 10−200 MeV band; this could be achieved by a dedicated mega-electronvolt–giga-electronvolt orbiting telescope such as e-Astrogam (De Angelis et al. 2017; Rando et al. 2019).

• Uncertainty in the level of optical luminosity. The largest uncertainty in the intrinsic luminosity of the CSB is on the optical emission, i.e., the level of de-reddening necessary to recover the intrinsic optical flux from the measured one. The lower the de-reddening – i.e., the dimmer the deduced intrinsic flux of the optical hump – the less prominent the Compton/optical hump (at log(ν/Hz)≃23.5): as a consequence, at log(ν/Hz)≈23.3 our SED models would no longer exhibit a local minimum (which is indeed shown by the LAT data) but rather an inflection point. Thus, adopting a lower optical luminosity would worsen the match between model and data in this frequency range. In summary, the secondary peak shown by the LAT data at log(ν/Hz)≈23.5 suggests a relatively high intrinsic optical luminosity, L_B ≳ 7.7 × 10⁴³ erg s⁻¹, similar to that derived from Blackman et al. (1979) who argued for a high color index – especially in the dusty CSB. As the optical emission of the optical component is described as a blackbody (the graybody model only applies to the FIR emission), its dilution factor is the standard blackbody one (see fn. 6).

• CRp. Our fitting CRp spectrum has an index, q_p ≈ 2.3, consistent with the one derived by the VERITAS Collaboration (2025). The limiting CRp energy resulting from our fits, E_p^max ≈ 7 TeV, is comparable to the maximum energies attained by CR accelerated by shocks in a uniform medium (with a mild dependence on density for n_H = 1 ÷ 10² cm⁻³) during the Sedov-Taylor phase of SN evolution. Such energies are an outcome of semi-analytical models of particle acceleration, which are based on kinetic simulations for a wide range of astrophysical shocks and account for the interplay between particle acceleration, magnetic-field amplification, and shock evolution (Diesing 2023; Bell et al. 2013).

The CRp energy density, u_p, in our models is in the range of 360−400 eV cm⁻³ (for n_H = 200 cm⁻³). An independent estimate of u_p may provide a consistency check. Combining the SN frequency, ${\mathcal{R}}_{\mathrm{SN}}$, with the CRp residency timescale, τ_res, in the source volume, $\mathcal{V}=7.9\times {10}^7$ pc³, and assuming a nominal value of the fraction η = 0.1 (e.g., Higdon et al. 1998; Tatischeff 2008) of the total SN energy, E_ej = 10⁵¹ erg per SN (Woosley & Weaver 1995), that is channeled to particle acceleration, the estimated CRp energy density is $u_p={\tau }_{\mathrm{res}}{\mathcal{V}}^{-1}{\mathcal{R}}_{\mathrm{SN}}\eta E_{\mathrm{ej}}$. The CRp energy losses mainly occur through superwind advection and p–p interactions, with a characteristic residency timescale of τ_res = (1/τ_wind + 1/τ_pp)⁻¹. In the latter expression, ${\tau }_{\mathrm{wind}}=r_s/v_{\mathrm{wind}}={10}^6(r_{\mathrm{s}}/0.1\mathrm{kpc}){(v_{\mathrm{wind}}/100{\mathrm{km}\mathrm{s}}^{-1})}^{-1}$ (Longair 1981), which, for v_wind = 500 km s⁻¹ (Shopbell & Bland-Hawthorn 1998), gives τ_wind = 6.8 × 10⁵ yr. Also, τ_pp = (σ_ppcn_p)⁻¹ which, for σ_pp ∼ 40 mb at a few tera-electronvolts (Kelner et al. 2006) and n_p = 200 cm s⁻¹ (Table 5), is τ_pp ∼ 1.3 × 10⁵ yr. Thus, τ_res ∼ 1.1 × 10⁵ yr. The SN rates ${\mathcal{R}}_{\mathrm{SN}}\sim 0.13$ yr⁻¹ (consistent with X-ray and radio estimates, see Iwasawa 2021 and Lacki et al. 2011) yield u_p ∼ 385 eV cm⁻³, in accordance with our SED modeling results.

• Radio spectrum. Given that the CRp spectrum is fully determined by the measured γ-ray spectrum, and that consequently its secondary yield is also, the synchrotron radio spectrum mainly constrains the primary CRe spectrum. In the frequency range 0.3−9 GHz, the emission from the CSB region is fit by a thermally absorbed α_r = 0.62 PL spectrum (Adebahr et al. 2013); thus, the emitting-CRe population has a spectral index of 2.24. The total emission is the superposed contributions by primary and secondary electrons in which the latter gets steeper than the total at ν ≳ 2 GeV (Fig. 2). This feature motivates our assumption of q_e ≤ 2.2 for primary CRe. Star formation is much less intense outside the CSR (Shopbell & Bland-Hawthorn 1998), so the influx of CR into this region is expected to be relatively week.

Fig. 2. Model radio spectra for the values of β considered in this study (thick solid curves), each including primary and secondary electron synchrotron emission (short-dashed and dotted curves, respectively) and thermal free-free radiation (long-dashed curves), are overlaid with data from Table 1 (black dots).

Owing to the varying lepton-to-hadron ratio (and therefore the primary-to-secondary electron ratio) as a function of the graybody model describing the FIR emission, the low-frequency thermal absorption, quantified by the emission measure (EM), also varies (Table 6). A quite minor ν^−0.1 thermal free-free emission also enters the model spectrum, calibrated using the measured Hα emission (Kennicutt et al. 2008) as a proxy.

• CRe. The deduced primary CRe spectrum, with PL index q_e ≲ 2.2, may appear surprisingly flat. Actually, the fact that it is very close to the injection spectrum is understandable: irrespective of details of the acceleration mechanism, the electron injection timescale can be gauged by the inverse of the SN rate, ${\tau }_{+}={\mathcal{R}}_{\mathrm{SN}}^{-1}$ where ${\mathcal{R}}_{\mathrm{SN}}\sim 0.1$ yr⁻¹ (Iwasawa 2021). The CRe energy loss time is ${\tau }_{-}=\gamma /\dot{\gamma }$, where $\dot{\gamma }$ is the sum of loss rates by electronic excitations, bremsstrahlung, synchrotron, Compton (detailed, e.g., in Appendix A of Persic et al. 2024, and using the matter and radiation/magnetic-energy densities in the CSB reported in Tables 5, 6), and advection (Eq. 7). The comparison indicates that the injection timescale is much shorter than the loss time at all energies, even at γ ∼ γ_max. It is not surprising then that the primary electron component of the radio spectrum reflects the injection spectrum.

The production timescales of secondary CRe depend on the timescale of π^± production following p−p interactions, t_π^±. The latter is t_π^± = 1/(κ_π^±n_Hσ_ppc), where κ_π^± = 0.25 is the fraction of kinetic energy of the incident proton transferred to π^±; assuming σ_pp = 35 mb as appropriate at a few tera-electronvolts (Kelner et al. 2006) and n_H = 200 cm⁻³ (Table 5), it is t_π^± ≃ 1.9 × 10¹³ s. This implies that the spectral distribution of secondary electrons mimics that of their parent CRp spectrum.

• Magnetic field. The value of the magnetic field (B) averaged over the CSB volume, i.e., co-spatial with the NT SED analyzed in this study, was derived by fitting the radio synchrotron emission generated by CRe whose spectrum was normalized to the NT 5 keV flux interpreted as Compton/FIR radiation (j_Compt ∝ n_e). This direct estimate of B is the first of its kind for any CSB. The deduced B values for our six models span the range B = 97–143 μG (Table 6).

Our estimates of B can be compared to earlier estimates that were typically based on the assumption of energy equipartition between CR and the magnetic field. As an example, in our model with β = 1, the magnetic field is B = 97 μG and the CRp-to-CRe number density ratio is p/e(1 GeV) = 120. Interestingly, Adebahr et al. (2013) found a similar result based on the radio spectrum alone: assuming particle-field equipartition according to Beck & Krause (2005) Eq. (A18) with p/e(1 GeV) = 120, they obtained B_eq ≈ 98 μG. Persic & Rephaeli (2014), too, estimated B = 100 μG from radio emission, using a revised equipartition formula that includes theoretically based p/e and u_p/u_e ratios (including secondary CRe) as function of the CR spectral indices, as well as environment-dependent energy losses. Using their own revised equipartition formula that included pionic secondaries and energy losses, Lacki & Beck (2013) estimated B = 240 μG. On the other hand, accounting for a single population of particles, Völk et al. (1989) estimated B ∼ 50 μG from minimum energy in CR plus the magnetic field. In their models of M82 multifrequency emission, Peretti et al. (2019) adopted B = 210 μG, while De Cea del Pozo et al. (2009) found B = 120 ÷ 290 μG accounting for uncertainties in the environmental parameters.

In all our models, the deduced values of B are very close to its equipartition value, $B\approx B_{\mathrm{eq}}=\sqrt{8\pi (u_{\mathrm{CRp}}+u_{\mathrm{CRe}})}$, where u_CRe includes primary and secondary contributions. Thus, rather than assuming equipartition ab initio, we deduced its approximate validity in the CSB region.

• Photon-photon absorption. The intense local IR photon field raises the question of whether any photon–photon absorption affects the emerging γ-ray emission. The relevant cross section (Heitler 1960) is σ_γγ(E_γ, ϵ, ϕ) = (3/16) σ_T (1 − β²) [2β(β² − 2)+(3 − β⁴) ln[(1 + β)/(1 − β)]], where $\beta =\sqrt{1-[2m^2{c}^4]/[E_{\gamma }ϵ(1-\cos \varphi )]}$, with E_γ the energy of the incoming photon, ϵ the energy of the local photon, and ϕ the interaction angle between their trajectories. Integrating over ϕ and the spectrum of target photons, n_IR(ϵ), and multiplying by the characteristic source size, r_s, we obtained the photon–photon optical depth, τ_γγ(E_γ). Figure 3 shows that in situ photon–photon absorption is negligible for incoming energies of ≲10 TeV, well above the energies sampled by the current VERITAS data.

Fig. 3. Photon–photon absorption, described as exp( − τγγ).

• Neutrino detectability. It is of interest to assess the detectability of the π^±-decay neutrino emission from M82 with current and upcoming neutrino observatories. With an apparent dominant π⁰-decay origin of γ ray emission in the CSB, π^±-decay neutrinos are (obviously) also produced (e.g., Ambrosone et al. 2021). Our calculation of the predicted ν_μ and ν_e spectra of M82, using (Kelner et al. 2006) formalism (Sect. A.2.3 of Persic et al. 2024), indicates that the broadly peaked giga-electronvolt–tera-electronvolt neutrino flux (Fig. 4) is too low for detection by current and upcoming ν projects. This conclusion is based on the following calculation of the estimated observation time needed to detect M82 with an experiment with detection sensitivity comparable to the Antarctica-based IceCube + DeepCore Observatory, the most sensitive current ν-detector at neutrino energies E_ν > 10 GeV (e.g., Bartos et al. 2013). We also considered the Mediterranean-based KM3NeT observatory, set to become a sensitive detector of E_ν ∼ 10–100 GeV atmospheric neutrinos and of E_ν > 1 TeV cosmic neutrinos (Adrián-Martinez 2016).

Fig. 4. Predicted neutrino spectral flux from M82 (β = 1.5 model). The neutrino spectra for all the models considered in this study look similar, due to the models’ very similar hadronic characteristics. The curves represent total-ν, νμ, νe spectra (thick solid, thin dotted, thin dashed lines, respectively).

– The IceCube + DeepCore μ_ν effective area (Abbasi et al. 2012, Fig. 8-right) can be approximated as A_eff(E_νμ) = 10^{4 − [0.65x3 − 4.5x2 + 6.95x + 1.3]} cm², where x = log(E_νμ) (energies in giga-electronvolts). (A_eff is ∼2 times smaller for ν_e.) Only in the energy range 10 GeV–2.2 TeV do the IceCube + DeepCore sensitivity and the M82 predicted diffuse spectral ν-flux effectively overlap. In this band the ν flux can be approximated as $\frac{\mathrm{d}{N}_{\nu }}{\mathrm{d}{E}_{\nu }}\sim {10}^{-8.34}{E}_{{\nu }_{\mu }}^{-q_p}$ cm⁻² s⁻¹ GeV⁻¹. The corresponding number of detected ν per year, $N_{\nu }=t_{1\mathrm{yr}}{\int }_{10\mathrm{GeV}}^{2.2\mathrm{TeV}}\frac{\mathrm{d}{N}_{\nu }}{\mathrm{d}{E}_{\nu }}{A}_{\mathrm{eff}}(E_{{\nu }_{\mu }})\mathrm{d}{E}_{{\nu }_{\mu }}$ (where 2.2 TeV is the cutoff of the neutrino spectrum, Fig. 4), is N_νμ = 2.45 yr⁻¹.

The detector background (for upward-going events) is dominated by atmospheric neutrinos produced by CR in the northern hemisphere. To compute such a noise background, N_ν, we started from the atmospheric ν_μ spectrum (Fig. 12 of Aartsen et al. 2015; the spectrum is ∼20 times lower), which in the energy range 10 GeV−2.2 TeV (see Fig. 4) can be approximated as ϕ(E_ν) = E_ν⁻² 10^{−(ax2 + bx + c)} Gev⁻¹ s⁻¹ cm⁻² sr⁻¹, where a = −0.203457, b = 0.17507, c = 2.03243, and x = log(E_ν) with E_ν in giga-electronvolts. So the number of vertical neutrinos detected by IceCube + DeepCore is $N_{\nu }=t_{1\mathrm{yr}}{\int }_{10\mathrm{GeV}}^{2.2\mathrm{TeV}}\frac{\mathrm{d}{N}_{\nu }}{\mathrm{d}{E}_{\nu }}{A}_{\mathrm{eff}}(E_{\nu })(\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{M}82}/1.6\mathrm{sr})\mathrm{d}{E}_{\nu }\simeq 24$ yr⁻¹, for a solid angle, ΔΩ_M82 ≃ 2.4 × 10⁻⁴, corresponding to a 0.5 deg radius uncertainty in the direction reconstruction of IceCube + DeepCore ν_μ toward M82⁸.

Based on this rough estimate, observations of diffuse giga-electronvolt–tera-electronvolt muon neutrinos from M82 would imply S/N ∼ 0.1. This estimate indicates that the likelihood of detecting M82 with an IceCube + DeepCore–like detector is very low.

– The upcoming KM3NeT observatory is composed of two Cherenkov neutrino telescopes under construction in the Mediterranean sea, ARCA and ORCA (Astroparticle/Oscillation Research with Cosmics in the Abyss, respectively). ARCA, to be located 100 km southeast of Sicily’s southernmost tip at 3500 m depth, is optimized for detection of ≳1 TeV astrophysical neutrinos. ORCA, to be located off the coast of Toulon (France) at a depth of 2500 m, is optimized for energies of ≈10 GeV in order to study neutrino properties exploiting neutrinos generated in the Earth atmosphere. The ARCA sensitivity for sources with an unbroken E⁻² spectrum for an observation time of 6 years is ≈10⁻⁹ GeV cm⁻² s⁻¹ for the full declination range −1 ≤ sin (δ)  ≲ 0.8, similar to the results, for a similar exposure, reported by IceCube for northern hemisphere sources (Aiello et al. 2019). Clearly, ARCA would be the choice KM3NeT instrument to observe M82. Indeed, Aiello et al. (2024) have shown that individual galaxies with ongoing star formation, with concurrent AGN activity (NGC 1068 and the Circinus galaxy) or without (the Small Magellanic Cloud), can be detected by ARCA in 10 years of observation if their CRp spectra are hard (q_p = 2) and unbroken to very high energies (500 TeV). However, our modeling of the LAT and VERITAS γ-ray data for M82 implies q_p ≃ 2.3 and E_p = 7 TeV for the CRp spectrum, which translate into a ≈1 TeV cutoff in the M82 neutrino spectrum. At E_ν ≈ 1 TeV the ARCA effective area is as low as A_eff ≈ 1 m² (Adrián-Martinez 2016), making ARCA unsuitable for detecting M82.

6. Conclusion

The recent publication of diffuse NT 4−8 keV emission from M82’s CSB motivates a different interpretation of the NT emission in this region. We complement the X-ray data with 50 MeV−16.25 TeV γ-ray data and 0.3−8.3 GHz radio data.

This analysis hinges on two points: first, identifying the high-energy γ-ray data as pionic emission: this fixes the CRp andsecondary CRe spectra; and second, combining secondary and primary CRe to normalize the latter by fitting the predicted Compton/starlight emission to the NT X-ray flux and determine its spectral slope and cutoff from the low-energy γ-ray data. Observational uncertainties on the FIR graybody emission result in a range of viable SED models.

The SED models are quite similar in these models, the NT parameters only showing moderate variations. The CRp and primary-CRe spectra are PL with spectral indices of q_p ≈ 2.3 and q_e ≈ 2.2; the CR energy densities are u_p ≈ 385 eV cm⁻³ and u_e ≈ 5 eV cm⁻³, respectively. The magnetic field, deduced from a synchrotron fit to the radio spectrum, is B ≈ 120 μG, in energy equipartition with the CR. More precise measurements of the FIR emission would narrow down the range of viable NT SED models.

Acknowledgments

We thanks Dr Lab Saha of the VERITAS Collaboration for providing the tabulated 2008−2022 VERITAS spectral data (plotted in VERITAS Collaboration 2025). MP acknowledges hospitality from the Physics & Astronomy Department of Padova University. We acknowledge useful and enlightening conversations with Jean Ballet and Philippe Bruel on some aspects of the LAT data analysis. The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supportedboth the development and the operation of the LAT as well as scientific data analysis. These include NASA and the Department of Energy in the United States, the Commissariat à l’Énergie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et the Physique des Particules en France, the Agenzia Spaziale Italiana (ASI) and the Istituto Nazionale di Fisica Nucleare (INFN) in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT) and the High-Energy Acceleration Research Organization (KEK) in Japan, and the K.A. Wallenberg Foundation and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica (INAF) in Italy and the Centre National d’Études Spatiales in France.

References

Appendix A: Fermi-LAT data analysis

The Large Area Telescope (LAT) is the main instrument on the Fermi Gamma-Ray Space Telescope. It comprises a silicon microstrip tracker, a cesium-iodide calorimeter, and a plastic anti-coincidence detector. The LAT covers an energy range from 20 MeV to ∼2 TeV with a field-of-view of 2.4 sr.

For the analysis described in this work we used 16.3 years of LAT data of the P8R3_SOURCE class. We considered events from 50 MeV to 300 GeV, within a region of interest centered on M82 and 15° in radius. We reduced the data set with the usual zenith angle cut < 90°, to limit the contamination from the bright Earth limb. The details of the LAT data selection are reported in Table A.1.

The global fit is performed with a binned likelihood analysis. The source model is based on the fourth LAT source catalog, data release 4, 4FGL-DR4 (Abdollahi et al. 2022; Ballet et al. 2023); the CSR cannot be resolved by the LAT and is modeled as a point source. Energy dispersion ⁹ and likelihood weighting (Bruel 2021) are taken into account in the analysis.

We started the analysis using a LogParabola spectrum for M82, following the spectral model used in the 4FGL. The spectral parameters α and β we derived were compatible with 4FGL within uncertainties. If the analaysis is restricted to energies > 200 MeV the agreement is excellent. To investigate the presence of a rising spectrum at lower energies, located close to the expected pionic shoulder, we tried including a second spectral component with a PL spectrum (modelled as a second point source, coincident with M82). The additional component resulted very steep, with spectral index ∼5, and significant, with a TS = 29.5. The new component is non-negligible only at energies ≲200 MeV. Fixing the spectrum for M82 to the catalog values, test-statistics maps indicate that the excess is approximately point-like and centered on top of M82. The parameters for this component are listed in Table A.2.

There are several sources of systematic uncertainties to consider. Energy dispersion corrections were applied to the spectra. We proceeded by increasing the number of edisp_bins in the global fit to −3, after which the fit become insensitive to an additional increase. The uncommonly steep spectrum of theadditional component at low energy magnifies any remaining uncertainty to an unprecedented level, though, and the usual estimates on the residual systematic uncertainties may be too optimistic. M82 is located close to the north celestial pole. In proximity of the celestial poles, there is an increased risk of contamination from photons from the Earth limb, especially at low energies due to the broader PSF. A loop of interstellar gas, part of the North Celestial Pole Loop, passes very close to M82 and a relatively brighter globular structure, much smaller than the LAT point spread function at low energies, lies very close to M82. It is worth noting that the gas structure in question may be a nearby SNR (Marchal & Martin 2023; Schmelz et al. 2023). Uncertainties in the emission model can affect the spectral estimate for M82 at low energy; uncertainties in modelling the diffuse models were addressed in part by applying log-likelihood weights.

To try and address the aforementioned issues we performed an additional analysis of the low-energy part of the spectrum, selecting only events belonging to the PSF3 class: LAT events are subdivided into quartiles depending on the quality of the direction reconstruction, and we restrict ourself to those with the best PSF. We divided the events in the energy range 50–213 MeV in four logarithmic bins, and proceeded with the binned likelihood analysis. In the first two bins, M82 is now not detectable with any significance, but the upper limits we derive are still relatively high. The difference with the previous estimates indicate severe contamination from sources other than M82. In bin 3, M82 is detected with TS = 20, and an excess with respect to the 4FGL spectrum is significant with TS = 6. In bin 4, M82 is detected with TS = 17 and the derived spectrum is compatible within uncertainties with the 4FGL value.

In this work we consider the 4FGL spectrum above 213 MeV. Below that energy, we have indications of an additional component, but systematics make it impossible to derive a spectral shape. We consider the 4FGL spectrum as a lower bound, and our best estimate (PSF3 U.L.) as an upper bound.

The SED was evaluated by dividing the whole energy range in 12 equally spaced logarithmic energy bins, and the first two low-energy bins were further divided in two, and finally by performing a dedicated binned likelihood fit in each bin. In the broader bins, we used the broader selection (SOURCE event class). The spectra of M82 was described as a PL with the fixed spectral index corresponding to the 4FGL spectrum considering the bin energy centroid. The SED data points resulted perfectly compatible with the ones distributed with 4FGL and with previous analyses (Ajello et al. 2020). In the finer, low-energy bins, we used the restrictive selection (SOURCE-PSF3 event class). We evaluated upper limits with the standard binned likelihood approach. The results are shown in Fig. A.1. The dashed line represents the 4FGL spectrum, the filled gray band indicates the 1-σ uncertainty. The data points are the SED we derived in this work, assuming the spectral shape above but using 16.3 years of data; the error bars indicate 1-σ uncertainties. The upper bounds are the PSF3 90% C.L. upper limits. The lower bounds at low energy are the 4FGL spectral values, evaluated at the bin energy centroids. The resulting LAT flux values and upper/lower limits are reported in Table A.3.

Fig. A.1. M82 spectral model. Dashed line: 4FGL spectrum with 1σ uncertainties (gray band). Data points: SED points derived assuming 4FGL spectrum, with 1σ uncertainties. Dotted line: 4FGL spectrum plus the low-energy PL component. Downward triangles: low-energy upper limits from the PSF3 analysis. Upward triangles: low-energy lower bounds, flux values from 4FGL model.

Appendix B: VERITAS data

VERITAS is a stereoscopic system of four 12-m–diameter imaging atmospheric Cherenkov telescopes, located at Whipple Observatory in southern Arizona, USA at an elevation of 1268 m above sea level. Its energy threshold is ∼100 GeV, and its sensitivity (i.e., 5σ detection) is ∼0.6% of the Crab Nebula flux in 50 hours of observation at small (< 30°) zenith angles (Park & VERITAS Collaboration 2015). Table B.1 lists the VERITAS spectral points described and plotted in VERITAS Collaboration (2025; see their Fig. 1).

All Tables

All Figures

Fig. 1. X-ray/γ-ray range of the CSB SED. The emission model (thick solid curve), including total (i.e., primary plus secondary) Compton/M82-EBL-CMB radiation (dashed curve) and NT bremsstrahlung (dash-dotted curve), and a pionic component (dotted curve), is overlaid with data from Tables 2–4 (black dots). The leptonic component is dominated by the Comptonization of the local FIR radiation. The 5 keV flux, resulting from the NT 4−8 keV flux measured by Iwasawa et al. (2023), is fit by (basically) the local Compton/FIR yield, in which the secondary electron spectrum is uniquely set by the CRp spectrum fit to the γ-ray data and the primary-electron spectrum is normalized to match the 5 keV flux.

In the text

	Fig. 2. Model radio spectra for the values of β considered in this study (thick solid curves), each including primary and secondary electron synchrotron emission (short-dashed and dotted curves, respectively) and thermal free-free radiation (long-dashed curves), are overlaid with data from Table 1 (black dots).
In the text

	Fig. 3. Photon–photon absorption, described as exp( − τγγ).
In the text

	Fig. 4. Predicted neutrino spectral flux from M82 (β = 1.5 model). The neutrino spectra for all the models considered in this study look similar, due to the models’ very similar hadronic characteristics. The curves represent total-ν, νμ, νe spectra (thick solid, thin dotted, thin dashed lines, respectively).
In the text

	Fig. A.1. M82 spectral model. Dashed line: 4FGL spectrum with 1σ uncertainties (gray band). Data points: SED points derived assuming 4FGL spectrum, with 1σ uncertainties. Dotted line: 4FGL spectrum plus the low-energy PL component. Downward triangles: low-energy upper limits from the PSF3 analysis. Upward triangles: low-energy lower bounds, flux values from 4FGL model.
In the text