Non-thermal emission in M31 and M33

Spiral galaxies M31 and M33 are Fermi/LAT-detected gamma-ray sources. We model the broadband non-thermal (NT) emission of the central region of M31 (R<5.5 kpc) and of the disk of M33 (R ~ 9 kpc). For either galaxy, we self-consistently model the broadband SED of the diffuse NT emission based on published radio and gamma-ray data. All relevant radiative processes involving relativistic and thermal electrons (synchrotron, Compton scattering, bremsstrahlung, and free-free emission and absorption), along with relativistic protons (neutral-pion decay following interaction with thermal protons), are considered, using exact emissivity formulae. We also use the Fermi/LAT validated gamma-ray emissivities for pulsars. We find that, in both sources, the radio emission is composed of primary and secondary electron synchrotron and thermal bremsstrahlung. The M33 gamma-ray emission appears to be mainly hadronic, similar to the Magellanic Clouds (Persic&Rephaeli 2022). In contrast, we find suggestions of a more complex situation in the central region of M31, whose emission could be a mix of pulsar emission and hadronic emission, with the latter possibly originating from both the disk and the vicinity of the nuclear black hole. The alternative modelling of the spectra of M31 and M33 is motivated by the different hydrogen distribution in the two galaxies: the hydrogen deficiency in the central region of M31 partially unveils emissions from the nuclear BH and the pulsar population in the bulge and inner disk. If this were to be the case in M33 as well, these emissions would be outshined by diffuse pionic emission originating within the flat central-peak gas distribution in M33.


Introduction
Non-thermal phenomena are mainly probed by the radiative yields of relativistic electrons and protons (hereafter simply referred to as electrons and protons, unless noted otherwise).Measurements of the emission spectral energy distribution (SED) provide direct insight on the particle energy spectra and on magnetic and radiation fields in star-forming environments.
In a previous paper (Persic & Rephaeli 2022: Paper I), we examined non-thermal emission in the Magellanic Clouds, satellite galaxies of the Milky Way.Adopting a one-zone model for the extended disk emission, we self-consistently modelled disk non-thermal emission at radio and γ-ray frequencies, interpreting radio data as synchrotron plus thermal bremsstrahlung, and γ-ray data as neutral-pion decay (following energetic proton interactions in the ambient gas) plus leptonic emission.
In this follow-up paper on the non-thermal properties of starforming galaxies, we consider the two Local Group member galaxies M31 and M33 -at distances of 0.78 and 0.85 Mpc, respectively. 1As it was for the Magellanic Clouds in Paper I, our main objective in this study is to determine relativistic electron and proton spectra in the disks of M31 and M33 by spectral modelling of their non-thermal emission in all the relevant energy ranges accessible to observations.In this paper, too, we focus on the extended disk emission of each galaxy because this emission traces the mean galactic properties of non-thermal particles and magnetic fields.Their integrated radio spectra have been measured long ago (Beck 2000 and references therein;Dennison et al. 1975), whereas it is only recently that γ-ray data have become available (Abdo et al. 2010, Xi et al. 2020) -but there are no published measurements of extended non-thermal X-ray emission.Adopting a one-zone model for the extended disk emission we calculate the radiative yields of electrons and protons and contrast these with available radio and γ-ray data.The latter are fully specified in Table 1 (M31) and Table 2

(M33).
In Section 2, we review published observations of broadband non-thermal interstellar emission from M31 and M33.In Section 3, we review the radiation fields permeating the two galaxies.In Section 4, we describe the SED modelling procedure.In Sections 5 and 6, we present and compare our SED models to the data for M31 and M33, respectively.We summarise our results in Section 7.

Observations of extended emission
M31 and M33 have been extensively observed over a wide range of radio-microwave, (soft) X-ray, and γ-ray bands; point-source, extended small-and large-scale emission has been detected over a wide spectral range.Their detected disk X-ray diffuse emission is thermal from hot (T e ∼ 10 6 K) ionised plasma.It is relevant to our SED analysis as it contributes to the ambient proton density that enters the calculation of the pionic γ-ray spectrum.Both galaxies were searched for γ-ray emission early on in the Fermi-LAT mission (Abdo et al. 2010): M31 first, and later M33, were detected as marginally extended sources (Abdo et al. 2010;Xi et al. 2020).
Emission from both galaxies is consistent with an empirical, quasi-linear L >0.1 GeV versus L 8−1000 µm correlation established for nearby star-forming galaxies (Ackermann et al. 2012), suggesting a SF-related origin of their γ-ray emission.In this section, we review the observations most relevant to non-thermal emission, referring for details to the cited papers.

M31
• Radio.Battistelli et al. (2019) deduced an overall skyintegrated2 density spectrum of M31 from the radio to the infrared (IR), based on a compilation of Sardinia Radio Telescope (SRT) 6.7 GHz measurements, newly analysed Wilkinson Microwave Anisotropy Probe (WMAP) and Planck data, as well as other ancillary data -after subtracting point sources at ≤3σ above noise estimation.Decomposing the radio-IR SED into synchrotron, thermal free-free, thermal dust, and the 'anomalous' microwave emission (AME) arising from electric dipole radiation from spinning dust grains in ordinary interstellar conditions (Draine & Lazarian 1998a,b; see Dickinson et al. 2018 for a review), Battistelli et al. (2019) found that at ν < 10 GHz the SED is dominated by synchrotron emission with a spectral index α = 1.10 +0.10 −0.08 (in agreement with Berkhuijsen et al. 2003) combined with a subdominant ∝ ν −0.1 free-free emission at the highest frequencies.At ν ∼ 10 GHz AME becomes comparable to synchrotron and free-free emissions, at 20 GHz < ν < 50 GHz it dominates the emission, and at ν > 100 GHz, thermal dust emission is virtually the only emission component.Thus our analysis of non-thermal emission will use the Battistelli et al. (2019) data compilation at low frequencies only (ν < 10 GHz; Table 1), where AME and thermal dust emission are negligible.
A peculiar feature characterises the (atomic and molecular) hydrogen distribution in M31.The central galactic region is poor • X-ray.Presence of diffuse emission in the bulge was suggested by Primini et al. (1993) and Supper et al. (1997), based on Röntgen Satellit (ROSAT) soft X-ray data after point-source subtraction3 , with a luminosity L 0.1−2.4keV< 3 10 38 erg s −1 .Unresolved emission from the disk was clearly seen in residual maps after subtraction of point sources, instrumental, and extragalactic backgrounds (West et al. 1997, based on ROSAT PSPC and HRI data).Two components were distinguished, a spherical core plus an exponential disk with a radial profile similar to that of the optical light (i.e. the surface brightness is described as I(R) ∝ e −(R/R d ) where R d = 5 kpc is the radial scalelength).The disk diffuse emission makes up ∼45% of the total diffuse emission and corresponds to a hot plasma mass of M X < 7 10 6 M ⊙ .Takahashi et al. (2001: Advanced Satellite for Cosmology and Astrophysics [ASCA] data) interpreted spatially averaged (out to galactocentric radii 12 ′ ) spectra as integrated emission (L 0.5−10keV = 2.6 10 39 erg s −1 ) from a population of low-mass X-ray binaries 4 plus two diffuse thermal plasma models characterised by k B T = 0.3 and 0.9 keV both exhibiting L 0.5−10keV = 2 10 38 erg s −1 .
Bogdan & Gilfanov (2008; Chandra and XMM-Newton data) detected soft emission from ionised gas with a temperature of k B T ∼ 0.3 keV and M X ∼ 2 10 6 M ⊙ , significantly distributed along the minor axis of the galaxy, possibly indicating outward flowing gas perpendicular to the disk.The vertical extent of the gas is ∼ > 2.5 kpc, as suggested by the shadow (suppressed emission) seen on the 10 kpc radius star-forming ring at ∼ 10 ′ from the centre.
This review of the X-ray literature on M31 has revealed no non-thermal X-ray flux that, interpreted as Comptonised cosmic microwave background (CMB) radiation, could have been used to calibrate the relativistic electron spectrum (e.g.Persic & Rephaeli 2019a).The detected X-ray emission is thermal from hot (T e ∼ 10 6 K) ionised plasma, whose density contributes to the ambient density in the p-p interactions that lead to the production of neutral pions (π 0 , which quickly decay into γ rays) and charged pions (π ± , which ultimately produce secondary electrons and positrons) as well as non-thermal bremsstrahlung.
• γ ray.M31 was detected at 5σ significance with ∼ < 2 yr of Fermi-LAT data as a marginally extended source in the 0.2−20 GeV band (Abdo et al. 2010), and then again (at ∼ < 10σ) with >7 yr of LAT Pass 8 data as an extended source in the energy range 0.1−100 GeV (Ackermann et al. 2017)  5 .The γ-ray brightness distribution was found to be consistent with a 0 • .4-radiusuniform-brightness disk and no offset from the galaxy centre, 3 86 sources in the central ∼34 ′ down to L 0.2−4keV = 1.4 10 36 erg s −1 from a 48 ks observation with the High Resolution Imager (HRI; Primini et al. 1993); and 396 sources within 6.3 deg 2 in the luminosity range to 3 10 35 < L 0.1−2.4keV/erg s −1 < 2 10 38 from a 205 ks observation with the Position Sensitive Proportional Counters (PSPC; Supper et al. 1997), respectively. 4Hard (3−100 keV) X-ray data the emission shows two components, one broadly peaking at ∼5 keV, and another extending to >100 keV, that are interpreted as arising from NS binaries with, respectively, high/low accretion rates and luminosities above/below ∼10 37 erg s −1 (Revnitsev et al. 2014: Rossi X-ray Timing Explorer (RXTE)/Proportional Counter Array (PCA), INTErnational Gamma-Ray Astrophysics Laboratory (INTEGRAL)/Integral Soft Gamma-Ray Imager (ISGRI), and Neil Gehrels Swift Observatory/Burst Alert Telescope (BAT) data). 5Based on 33 months of High Altitude Water Cherenkov (HAWC) telescope data, an upper limit was set on the 1−100 TeV flux from the galactic disk of M31 (Albert et al. 2020).
Article number, page 2 of 19 Massimo Persic,Yoel Rephaeli and Riccardo Rando : Relativistic particles in M31 and M33 with a spectrum compatible with either a power-law (PL) profile (index 2.4 ± 0.1) or a hump -more likely to originate from π 0 decay than from Compton scattering of starlight.The emission appeared to originate from the inner regions of the galaxy, and seemed not to be correlated with regions rich in gas or ongoing star-formation activity, suggesting (according to Ackermann et al. 2017) a non-interstellar origin unless the energetic protons originated in previous episodes of star formation.On the other hand McDaniel et al. (2019), using multi-frequency data, suggested that the emission most likely has a lepto-hadronic (LH) origin, namely, a combination of pionic and Compton scattering of primary and secondary electrons off ambient radiation fields, with deduced parameter values consistent with previous studies of M31 and cosmic-ray physics.We note also that Di Mauro et al. (2019) argued against dark-matter decay as a source for the measured emission.
In a spectro-morphological analysis employing templates for the distribution of the stellar mass of the galaxy and updated astrophysical foreground models for its sky region, Zimmer et al. (2022) constructed maps of the old stellar population in the bulge that mimicked the distribution of a population of millisecond pulsars (msPSR).Analysing LAT data using such templates they obtained a 5.4σ detection, which is comparable with the detection by Ackermann et al. (2017) using a disk template.Zimmer et al. (2022) argued that ∼106 unresolved msPSR could account for the measured γ-ray luminosity and spectrum; the latter then interpreted as arising from Compton upscattering of ambient photons by electrons in the bulge that were accelerated by msPSRs in the disk.Eckner et al. ( 2018) compared the extended emission detected by Ackermann et al. (2017) with that of a msPSR population in M31 obtained from the one detected in the local Milky Way disk via rescaling by the respective stellar masses of the systems: with no free parameters, they estimated a contribution of ∼1/4 of the M31 emission, as well as the energetics and morphology of the Galactic Center 'GeV excess'.Fragione et al. (2019) estimated that msPSR in old globular clusters (that were likely dragged to the centre by dynamical friction) may contribute ∼1/8 of the emission deduced by Ackermann et al. (2017).
Using 7.6 yr of LAT data, Karwin et al. (2019) carried out a detailed study of the 1−100 GeV emission toward M31's outer halo with a total field radius of 60 • centred on M31.Accounting for foreground γ-ray emission from the Milky Way, they suggested the existence of an extended excess unrelated to the Milky Way and having a total radial extent of ∼ < 200 kpc from the centre of M31.While not ruling out an extended cosmic-ray halo underlying such purported emission, citing GALPROP predictions Karwin et al. (2019) argued against it based on the radial extent, spectral shape, and intensity of the large-scale emission.However, Do et al. (2021) and Recchia et al. (2021) pointed out that the emission radial extent and intensity strongly depend on assumptions on cosmic-ray diffusion outside the plane and in the halo of M31 and that the spectral shape is affected by foreground Milky Way emission and by the intrinsic weakness and limited statistics of the signal, so the emission out to ∼ < 200 kpc may indeed originate from an extended cosmic-ray leptonic and/or hadronic halo if cosmic-ray propagation is more effective than assumed by Karwin et al. (2019).Based on a model-dependent approach, Zhang et al. (2021) used the purported halo emission to estimate the total baryon mass contained in the circumgalactic medium within M31's the ∼250 kpc virial radius of M31 and suggested it could account for <30% of the missing (with respect to the cosmic abundance) baryons.
Most recently Xing et al. (2023) carried out an analysis of the expanded ∼14 yr LAT database, finding that the extended source claimed by Ackermann et al. (2017) actually consists of two separate regions: a central emission region, with a log-parabola spectrum, coincident with the optical centre of the galaxy, and a second region, with a PL spectrum, 0 • .4southeast of the centre (note: it is unclear whether the latter is a background extragalactic source or a local source associated with M31).Xing et al. (2023) argued that the newly revealed point-like nature of the central source necessitates a revised interpretation since the Ackermann et al. (2017) conclusion was based on the assumption that the central emission has an extended distribution (cosmic rays, pulsars).
This review of the γ-ray literature on M31 has shown that previous γ-ray observations were affected by source confusion.This motivates our new interpretation of the M31 central emission in Section 5.

M33
• Radio.A study of the large-scale radio emission over the wide frequency range 21-10,700 MHz was carried out by Israel et al. (1992), conflating original data (21-610 MHz) and published data (178-10,700 MHz).In light of the 18 ′ − 1 ′ beamsize range and point-source subtraction procedure, this database (Table 2) is still the most suitable for our analysis of the galaxy-wide integrated radio properties of M33 6 For example, the main focus of the recent deep 1.5 and 5 GHz Jansky Very Large Array (JVLA) survey of this galaxy (White et al. 2019) is a detailed mapping of point sources.Israel et al. (1992) used the method from Baars et al. (1977) to convert instrumental intensities to flux densities.The <100 MHz flux densities were obtained with the University of Maryland's (defunct) Clark Lake Radio Observatory: their uncertainties mostly stem from the low surface brightness of the radio disk and the relatively high confusion level for a large beam at low frequencies.The 151.5 MHz flux density was obtained from several observations carried out with Cambridge University's Mullard Radio Astronomy Observatory, and the 327 and 610 MHz flux densities were obtained from maps made with the Westerbork Synthesis Radio Telescope.
Virtually all the M33 emission analysed by Israel et al. (1992) is contained in a box of the following size: 45 ′ × 85 ′ .These authors noted a spectral turnover at ∼750 MHz, and proposed a double-PL (2PL) spectral model with low/high index of 0.18/0.86(compatible with previous works in the literature, based on much sparser data).To explain this spectral shape they suggested free-free absorption of the synchrotron emission by a cool (≤ 1000 K) ionised plasma to be a very likely possibility.
The circularly integrated HI surface distribution shows a broad (R ∼ < 6 kpc) flat central peak.At larger radii the distribution declines steeply and cuts off at ∼10 kpc (Wright et al. 1972).
• X-ray.Diffuse emission in M33 was first suggested based on Einstein Imaging Proportional Counter (IPC) and HRI data.In addition to a dominant nuclear source and 17 point-like sources located mostly in the spiral arms (Long et al. 1981;Markert & Rallis 1983), Trinchieri et al. (1988) identified a diffuse lowsurface-brightness emission from the galactic plane to which source confusion and hot (kT < 1 keV) gas could contribute.Using deeper ROSAT HRI data, Schulman & Bregman (1995) detected more point sources above their data sensitivity threshold and noted that the (residual) extended X-ray disk emission (with a luminosity of L 0.1−2.4keV∼ 10 39 erg s −1 ) is not due to (additional) point sources.
Finally, Long et al. (1996) used deep (50.4 ks) ROSAT PSPC data to expand the list of individual sources to 577 and confirmed the diffuse disk emission (possibly tracing the spiral arms) to be of thermal (kT ∼ 0.4 keV) origin.Using XMM-Newton data ∼ > 10 times deeper than earlier ROSAT observations, Pietsch et al. (2004) and Misanovic et al. (2006) detected 350 sources in the galaxy sky area (of which ∼200 are in M33) and again found evidence of a diffuse emission.A subsequent very deep (1.4 Ms) Chandra Advanced CCD Imaging Spectrometer (ACIS) survey of M33 (Plucinsky et al. 2008) revealed soft diffuse emission from the large-scale hot ionised plasma and specific extended regions (e.g. the star-forming HII regions NGC 604, with kT ∼ 0.2 keV; or IC 131 of enigmatic spectral interpretation : Tüllmann et al. 2009).
As in the case of M31, the conclusion of this survey of the Xray literature on M33 is that only thermal X-ray emission from hot plasma has been detected thus far.The plasma density is (obviously) used in the calculation of the outcome of p-p interaction and non-thermal bremsstrahlung.
• γ ray.M33 was not detected in the ∼ < 2 yr of Fermi-LAT data (Abdo et al. 2010), nor in the >7 yr of data (Ackermann et al. 2017).Karwin et al. (2019) saw positive residual emission in the direction to M33 while examining the γ-ray emission from the outer halo of M31.Di Mauro et al. (2019) claimed a point-like source detection of M33; however, this was challenged by Xi et al. (2020) on grounds that the Di Mauro et al. background model had missed three background sources positionally close to M33; if unaccounted for, these would be attributed to M33.Finally, using a state-of-the-art background emission model and 11.4 yr of LAT data, Xi et al. (2020) claimed a detection of M33: the LAT image shows an extended signal peaking on the northeast region of the galaxy, where its most massive HII region, NGC 604, is located.We use the Xi et al. (2020) data (Table 2).

Radiation fields
A reasonably precise determination of the ambient radiation field is needed to predict the level of γ-ray emission from Compton scattering of target photons by the radio-emitting electrons (and positrons).The total radiation field includes cosmic (background) and local (foreground) components.
Relevant cosmic radiation fields include the CMB, a pure Planckian with a local temperature of T CMB = 2.735 K and energy density of u CMB = 0.25 eV cm −3 , and the extragalactic background light (EBL).The latter originates from direct and dust-reprocessed starlight integrated over the star formation history the Universe, with two respective peaks that are referred to as the cosmic optical background (at ∼1 µm), and cosmic infrared background (at ∼100 µm).The two peaks are described as diluted Planckians, characterised by a temperature, T, and a dilution factor, C dil .The latter is the ratio of the actual energy density, u, to the energy density of an undiluted blackbody at the same temperature, T (i.e.u = C dil 4 c σ SB T 4 , where σ SB is the Stefan-Boltzmann constant).The widely used EBL model, based on galaxy counts in several spectral bands (Franceschini & Rodighiero 2017), can be numerically approximated as a combination of diluted Planckians, where A j and T j are suitable dilution factors and temperatures. 8ocal radiation fields in the two galaxies (galactic foregound light, GFL) arise from their stellar populations.Similar to the EBL by shape and origin, the GFL is dominated by two thermal humps, IR and optical.Full-band luminosities of these thermal components are needed to determine n(ǫ), the spectral distributions of the photons to be Compton scattered (shown in Eq.A.8).We compute the total IR (8 − 1000 µm) luminosities using In-fraRed Astronomical Satellite (IRAS) flux densities at 12µm, 25µm, 60µm, and 100µm (sampling the IR hump), from the relation f IR = 1.8 • 10 −11 (13.48 f 12 µm + 5.16 f 25 µm + 2.58 f 60 µm + f 100 µm ) erg cm −2 s −1 (Sanders & Mirabel 1996).The total optical luminosities are computed from narrow-band luminosities (H band for M31, B band for M33) by applying suitable bolometric corrections that restore the full-band luminosities.The corresponding energy densities, with k denoting either IR or optical, are: (a) M31: The bolometric IR luminosity, L IR = 1.2 10 43 erg s −1 , is computed from the IRAS flux densities reported in Table 3.The bolometric optical luminosity, L opt = 1.2 10 45 erg s −1 , is computed from the narrow-band 3.6µm IRAC luminosity, L IRAC = 1.17 10 11 L ⊙ (Courteau et al. 2011), multiplied by the bolometric correction, f bol = 2.7 (Rosenfield et al. 2012).Our region of interest is within R ⋆ = 1.09R d , within which only 30% of the galaxy-wide emission is produced and is relevant to our GFL computation. 9Given the relative thinness of the radio disk, only a negligible fraction of the radiation emitted at R > R ⋆ enters the inner region.We obtain u IR = 0.04 eV cm −3 and u opt = 4 eV cm −3 .The IR and optical hump temperatures reported in Table 3 imply dilution factors for the M31-sourced photon fields of c IR = 10 −4.925 , and C opt = 10 −10.917 ; (b) M33: The bolometric IR luminosity, L IR = 4.8 10 42 erg s −1 , is computed from the data reported in Table 3.The total (bolometric) optical luminosity, L opt = 3.1 10 43 erg s −1 , is computed from the de-reddened blue magnitude and color index, B 0 T = 5.92 and (B − V) 0 ∼ 0.47 (Rice 1988), by applying Buzzoni et al.'s (2006) bolometric correction as described in Persic & Rephaeli (2019a).Thus, u IR = 0.022 eV cm −3 and u opt = 0.14 eV cm −3 .The IR and optical hump temperatures (Table 3) imply c IR = 10 −5.267 , and C opt = 10 −12.436 .

Xi et al. (2020)
⋆ Considered to be a lower limit (Israel et al. 1992).Rice et al. (1990).Temperatures T IR and T opt estimated using Wien law from the IR peak at 100 µm and the optical peak at 1.1 µm (Rice et al. 1990).
and nearby galaxies where the spatial profiles of key quantities, such as the particle acceleration sources, gas density, magnetic field, and generally also the particle diffusion coefficient may be specified (e.g.Strong et al. 2007Strong et al. , 2010;;Jones 1978;Taillet & Maurin 2003;Heesen 2021).However, if the aim is gaining insight on spatially-averaged non-thermal properties in the disk, then a parameter-intensive spectro-spatial modelling approach is not warranted.While 'top-down' diffusion-based modelling of the non-thermal emission in the disks and halos of galaxies is possible (Syrovatskii et al. 1959;Wallace et al. 1980;Rephaeli & Sadeh 2019), here we adopted a 'bottom-up' approach.Thus, we started from the non-thermal yields of relativistic particles and we deduced their steady-state spectra, making it is possible (in principle) to deduce the injection spectrum.
Radio emission in disks of spiral galaxies is mostly (electron) synchrotron in a disordered magnetic field whose mean total value, B, is taken to be spatially uniform in the regions of interest10 plus thermal free-free from a warm (∼10 4 K) ionised plasma (Spitzer 1978).This emission may be absorbed, at low frequencies, by cold ( ∼ < 10 3 K) plasma (Israel & Mahoney 1980;Oster 1961).11A significant proton component contributes additional radio emission by secondary e ± produced by π ± decays and γ-ray emission from π 0 decay (Stecker 1971).In addition, γray emission is produced by electron scattering off photons of local and background radiation fields (Blumenthal & Gould 1970).The calculations of the emissivities from all these processes are standard.
We just briefly expand here on thermal absorption of radio emission since it was not discussed in Paper I. A key feature in modelling the radio data is accounting for free-free ab-sorption of non-thermal synchrotron by a warm thermal ionised plasma, as suggested by a high value of the Hα emission measure, EM ≡ n 2 e dl.This was originally pointed out for a sample of highly inclined spiral galaxies, including M33 and M31, based on the low ratio between the observed radio intensities and the intensities extrapolated from measurements at higher frequencies, with the ratio being smallest for edge-on galaxies (Israel & Mahoney 1980).The frequency-dependent optical depth for free-free absorption by thermal plasma is (Israel & Mahoney 1980): whereby where the electron temperature T e is in K, the frequency ν is in GHz, and the emission measure of the absorbing gas, EM, is in cm −6 pc.We assume that the radio emission and absorption are well mixed throughout the galaxy, so the emerging intensity is , where j ν is the (synchrotron plus thermal free-free) spectral emissivity and τ f f (ν) curves the spectrum at low frequencies.We assume the particle spectral distributions to be timeindependent and locally isotropic; then we have: (a) The proton spectrum is: (energies in GeV).The quantities N p,0 , q p are free parameters whose values will be determined by fitting the model to the data.(b) Secondary electrons are produced from π ± decays following p-p interactions of the protons with ambient gas: their spectrum has no free parameters once the proton spectrum is specified.A summary of secondary electron production is given in Section A.2.2.The secondary electron spectrum will be analytically approximated as where q 1 , q 2 and γ b1 , γ b2 are the low-and high-energy spectral indices and breaks; and η gauges the steepness of the high-end cutoff.Their numerical values will be obtained by fitting Eq. 5 to the actual spectrum (Eq.A.29) which does not have any degrees of freedom.Thus, in our treatment, these are not free parameters.
(c) Primary electrons, required to fit the synchrotron data together with the secondary electrons, are chosen such that they approximate the steady-state electron spectrum following injection.Electron spectra are expressed as functions of the Lorentz factor, γ.
The emission spectrum from the π 0 -decay has more constraining power than a generic PL given its characteristic 'shoulder' at ∼ < 100 MeV and its spectral cutoff corresponding to E max p .Thus, as in Paper I our modelling procedure begins with fitting a pionic emission profile to the γ-ray data with free normalisation, slope, and high-end cutoff.We then use the deduced (fully determined) secondary-electron spectrum together with a (modelled) primary electron spectrum to calculate the combined synchrotron emission using an independently estimated value of B. The fit to the radio data also includes, at high frequencies, a thermal bremsstrahlung component computed with literature-deduced values of the gas density and temperature, and thermal absorption at low frequencies.Finally, the full nonthermal bremsstrahlung and Comptonised-starlight γ-ray yields are determined using standard emissivities (see Appendix A) and added to the pionic yield to fit the γ-ray data.In principle, if the resulting model exceeds the LAT data, the proton spectrum is revised accordingly and the procedure is repeated until convergence: in practice, for M31 and M33, the leptonic γ-ray yields turn out to be quite small compared to the pionic yields, so no iteration is needed.

M31
In Ackermann et al (2017), based on 7.3 yr of Fermi-LAT data, M31 was detected as extended with a significance of 4σ, with the favored model based on a disk of 0 • .4radius.This result provides an estimate of the LAT sensitivity to source extension given the LAT instrumental response, the source spectrum, the diffuse backgrounds in the region, and the available statistics.
In Xing et al. (2023), based on 14 years of LAT data, the measured emission was resolved into two regions, but now with a statistically significant evidence to reveal or rule out spatial extension.It is reasonable to assume that the γ-ray source associated with the M31 core has an extension (radius) of (at most) 0 • .4-otherwise it would appear as extended in Xing et al. (2023), as found in Ackermann et al. (2007).In this subsection, we test the hypothesis of an extended γ-ray emission in the unresolved central region of M31, arising from truly diffuse LH emission and possibly an unresolved population of pulsars.

γ-ray yields
The average ambient proton density is n p = n HI +2 n H 2 +n i ≃ 0.53 cm −3 (based on values listed in Table 4).A pionic fit, described by the parameters q p , E max p , u p (reported in Table 5) matches the log-parabola best-fit from Xing et al. (2023) to the LAT data.In a fully pionic model, the unusually high value of u p stems from the low hydrogen content in the central region of M31 (Fig. 12 of Nieten et al. 2006; the 'central deficiency' noted early on by, for instance, Davies &Gottesman 1970 andEmerson et al. 1974): the proton energy-loss time by p-p interactions is ∝ n −1 p , so for a given proton injection rate a low ambient gas density implies a high u p .
The closely related π ± -decay secondary electron spectrum, N se (γ), is computed as outlined above.Once the secondary synchrotron radiation has been computed (assuming B = 7.5 ± 1.5 µG; Fletcher et al. 2004, from particles-field equipartition), its low-frequency residuals from the radio data are modelled as synchrotron by a PL primary electron spectrum: in the range of γ min < γ < γ max (with γ min >> 0), where the quantities N e,0 , q e are free parameters.The ensuing pri- The appropriate model for the LAT-measured γ-ray spectrum of M31 might be not just LH.The galaxy bulge and disk are characterised by a half-light radius of R e = 1 kpc and an exponential scale length of R d = 5.3 kpc (Courteau et al. 2011), respectively.The upper limit to the central source radius, R ⋆ = 5.45 kpc, includes the bulge and ∼30% of the (thin exponential) disk.So we may expect a γ-ray contribution from a fairly large unresolved population of γ-ray pulsars (PSR) located in the bulge and the inner disk 12 .The γ-ray PSR in M31 are expected to be similar to those in our Galaxy.The latter average spectral parameters are Γ = 1.6, E cut = 2.5, and b = 1 (Abdo et al. 2013; Lopez et al.   12 If the LAT data need to be modelled as pure pulsar emission, a decent match is obtained for N psr = 3500 Galactic-like pulsars emitting, through a photon spectrum dN/dE ∝ E −Γ e −(E/Ecut) b (Abdo et al. 2013), where Γ = 2, E cut = 3 GeV, and b = 0.6 at a luminosity level L >0.1 GeV = 5 10 34 erg s −1 .2018).These values define a spectral shape inconsistent with the M31 LAT data.We incorporate the PSR spectrum and the LH spectrum into a mixed pulsar-plus-LH (PSR+LH) scenario.Possible examples of PSR+LH models involve N psr = 500 (1000) PSR (Fig. 1, middle and bottom; and Table 5); since the M31 central source has L >0.1 GeV = 1.4 10 38 erg s −1 (Xing et al. 2023), PSR contribute 18% (36%) of the source luminosity in these models.Eyeball fits suggest that this family of PSR+LH models can nicely match the data for N psr∼ < 1200 (and u p∼ > 5 eV cm −3 ).If the unresolved PSR population is composed of msPSR, models using average Galactic msPSR γ-ray parameters, namely, L >0.1 GeV = 3 10 33 erg s −1 , Γ = 1.3, E cut = 2, and b = 1 (Abdo et al. 2013) can accommodate up to ∼10 4 msPSR.

Radio yields
The radio emission measured in Battistelli et al. (2019) was integrated over the whole galaxy.To determine the primary electron densities in the region appropriate for calculating the γ-ray leptonic yields within R ⋆ , we need to estimate the radio flux from this region.To do so, we should know (in principle) the radial behavior of all the flux densities used in our analysis in order to build the radio SED integrated out to R ⋆ .From the literature, we took the radial profiles of the 408 and 842 MHz flux densities (Gräve et al. 1981, their Fig.9).Both profiles are point-sourcesubtracted and averaged over circular rings centred on the nucleus of M31 and have a comparable overall appearance.As the 842 MHz profile appears to have a more regular, Gaussian-like shape (with a standard deviation of σ ≃ 1 • ), we assume that it is an adequate representative flux density profile.Integrating over a circular R ⋆ region, we estimate that ∼15% of the total flux emanates from this region.Accordingly, our spectral fitting will be based on these reduced fluxes.Given the thinness of the disk, radiation produced in the outer disk is hardly captured by the <R ⋆ region.
Reducing the flux according to the 842 MHz profile is a firstorder correction to the galaxy-integrated flux.We expect fluxdensity radial profiles to steepen with increasing frequency because the corresponding higher-energy electrons diffuse out to shorter distances from their acceleration sites during their shorter energy loss time. 13Moreover, the injection rate of freshly accelerated electrons is not sufficient to offset for their shorter radiative lifetimes because star formation rate declines fast with radius -even faster than the gas density (based on the Schmidt-Kennicutt law Σ ∝ σ N with N > 1). 14The flux-density radial profiles at different frequencies have been measured in several nearby galaxies (e.g.M33: Tabatabaei et al. 2007;M51: Mulcahy et al. 2014) and a systematic (albeit usually moderate) radial steepening of the flux density profiles with increasing frequencies is clearly seen.We also expect this to be the case in M31 even though we are not aware of any available highfrequency data for this galaxy: if so, the emission fraction within R ⋆ will be higher at GHz frequencies than the value we obtained at 842 MHz.So to compensate for the overcorrection inherent in the low-frequency-based first-order correction, which (as noted) is progressively more severe with increasing frequency, a second-order correction is needed.Lacking suitable spectral data, the strength of such second-order correction is uncertain: examples from M33 and M51 suggest a ∼ < 25% increase at the highest frequencies.Thus, as the (unknown) second-order correction is probably small, we limited our correction to the firstorder and calibrated the ≤R ⋆ radio emission using the 842 MHz flux density profile.The main results of this paper will not be appreciably affected.For example, a 25% second-order correction can be compensated by a ∼40% increase in the thermal freefree emission leaving the non-thermal components essentially unchanged; alternatively, a flatter (q e = 3) and lower (by a factor of ∼4) electron spectrum would leave the pionic component virtually unaffected.
We wish to reproduce the M31 radio spectral model presented by Battistelli et al. (2019) self-consistently in the context of our approach.Synchrotron emission by secondary electrons (no degrees of freedom) dominates the radio emission up to ∼2 GHz.Fitting to the synchrotron by primary electrons results in a spectral index, α r ≃ 1.1, matching the total synchrotron fitting value of 1.1 +0.10 −0.08 from Battistelli et al. (2019).With no evidence of a spectral cutoff in the data, we set the primary electron maximum energy at γ max = 6 10 4 .The synchrotron spectrum from Battistelli et al. turns over at ν ∼ 48 MHz; as we do for M33 (Section 4.2.2),we interpret this feature as absorption by a cold plasma cospatial and homogeneously mixed with the electrons (Spitzer 1978); this is a common feature in highly inclined disk galaxies (Israel & Mahoney 1980).M31's perpendicular emission measure EM 0 ∼ 30 cm −6 pc (Walterbos & Braun 1994) implies, for a galaxy inclination of 77 • , a line-of-sight value of EM= EM 0 /cos(77 • ) = 133 cm −6 pc.Assuming a temperature T c e = 550 K for the cold plasma (similar to M33, see Section 4.2.2), from Eqs.( 2) and (3), we can derive an optical depth, which matches the model from Battistelli at al. at ν < 100 MHz.A ν −0.1 component in the Battistelli et al. model accounts for a slight spectral concavity at ν > 1 GHz.We modelled it as thermal free-free emission in the parametrisation reported in Eq.A.12, which gauges this flux with the Hα flux assuming that both emissions come from the same emitting (HII) regions with the same plasma parameters (temperature, density, filling factor), so the measured Hα flux may be used to predict the free-free emission.The radio spectrum is shown in Fig. 2, overlaid with the Battistelli et al. ( 2019) fit (top) and our model (bottom).

Broad-band results and discussion
The broadband SED models of the M31 central source are shown in Fig. 3.The fitting parameters are given in Table 5.
The small statistics of the available γ-ray data (five LAT spectral points) do not allow for a clear statement about the pulsar contribution, which however appears to marginally improve the fit at ∼2 GeV.Indeed, a consistent population of γ-ray pulsars is very likely to exist in M31: in the ("twin") Milky Way galaxy, about 300 pulsars have been detected by Fermi-LAT as of August 202315 and thousands of unresolved msPSR have been proposed as candidates to interpret the Galactic Center "GeV excess" (Brandt & Kocsis 2015, Bartels et al. 2016, Arca-Sedda et al. 2018, Fragione et al 2018a,b, Miller & Zhao 2023, Malyshev 2024; for a different view see Crocker et al. 2022).
The values of u p obtained from the SED modelling are much higher than the local Galaxy value of 1 eV cm −3 (e.g.Webber  2019) fit (single electron population synchrotron: dotted curve; thermal free-free emission: dot-dashed curve), and (bottom) our LH emission model.The latter includes primary and secondary synchrotron radiation (dotted curves of increasing flux at their peaks) and a thermal free-free component (dot-dashed curve).

1987
).An independent estimate (see Section 5) suggests u p ∼ 5 cm −3 for the average proton energy density within R ⋆ .If so, the PSR+LH model with 1000 PSR may be more realistic.
We note a similarity between the values of the proton spectral index in M31, deduced from our SED modelling, and in the Milky Way, from the measured cosmic-ray flux at the Sun position (e.g.Gaisser 1990).The local Galactic value, q p = 2.7, is understood as resulting from injection and diffusion, q p = q inj +δ with q inj ∼ 2.2 (typical of supernova remnants like the Crab Nebula) and δ ∼ 0.5 (e.g.Ptuskin 2006).
Any current view on the nature of the central engine may change if the Fermi/LAT spatial resolution may improve, for instance, as a consequence of a long future operating lifetime.The central source could either be resolved into multiple sources or confirmed to be pointlike at the centre.As to the latter possibility, in the next paragraph we briefly discuss the possibility that the pointlike source at the centre is identified with M31's nuclear supermassive black hole (BH), M31 * .
A comparison with Sagittarius A * (Sgr A * , the supermassive BH at the Galactic centre) may be useful.Cafardo & Nemmen (2021) found that the centroid of the 0.1-500 GeV emission, measured by Fermi/LAT (∼11 yr of data), approaches Sgr A * as the energy increases (at 10 GeV the two positions coincide within 1σ): at the Galactic Center distance, Cafardo & Nemmen's (2021) point-like source has L 0.1−500GeV = 2.6 10 36 erg s −1 , a value consistent with the Sgr A * bolometric luminosity.The mass of M31 * is 7.5 10 7 M ⊙ (Ford et al. 1994), whereas that of Sgr A * is 4.3 10 6 M ⊙ (Ghez et al. 2008).As for the X-ray luminosities, that of M31 * is 2 10 36 erg s −1 (Garcia et al. 2010, Li et 2019) have been renormalised to the region encompassed by R ⋆ (i.e.∼15% of the total emission; see Section 4.1.2).Emission components are denoted by these line types: radio synchrotron, dotted; thermal free-free emission, dot-dashed; total radio, solid; Comptonised CMB, short-dashed; Comptonised starlight (EBL+FGL), long-dashed; non-thermal bremsstrahlung, dotted; pionic, solid.Leptonic emission includes that from primary and secondary electrons(secondary components are dominant).In the X-ray/γ-ray range, on top of the separate components, a thick solid line indicates the total emission.Same for a broad-band PSR+LH SED (middle).The 500 pulsars emission component (see text) is ahown as a dot-dashed line.Same as above but for 1000 pulsars (bottom).Note: in this case the secondary leptonic components are subdominant.al. 2009), whereas that of Sgr A * fluctuates between ∼10 33 (baseline) and ∼10 35 during the hr-timescale daily flares (Sabha et al. 2010).Thus, both M31 A * and Sgr A * emit at level ∼ < 10 −10 L Edd , with L Edd = 1.26 10 38 (M/M ⊙ ) erg s −1 the Eddington luminosity.This suggests that both M31 * and Sgr * are currently dormant supermassive BHs.Therefore, if the pointlike source identified by Cafardo & Nemmen (2021) at the Galactic Centre is the γ-ray counterpart of Sgr A * , with the γ rays produced by relativistic particles accelerated by (or in proximity of) the BH (Cafardo & Nemmen 2021;Ajello et al. 2021) and if a scaling relation of yet unspecified nature (Aharonian & Neronov 2005) exists between the γ-ray luminosity and the BH mass, then for the γray counterpart of M31 * we may expect L 0.1−500GeV = 5 10 37 erg s −1 ; this value that corresponds to ∼1/2 of the luminosity of the central source identified by Xing et al.'s (2023).In this case, the pointlike central source in M31 may be related to M31 * .A study of the flux variability (t var = R s /c = 10 4 M BH /(10 9 M ⊙ ) s) may clarify the issue.Lacking a clear prediction on the γray spectral signature from a supermassive BH, we generally expected the γ-ray emission of M31 * to be partly pionic, from relativistic protons accelerated near the BH (perhaps by winds; see Ajello et al. 2021, Kimura et al. 2021) and interacting with the inflowing gas, and partly leptonic from the Comptonisation of the optical/X-ray radiation in the immediate proximity of the BH (coming from, e.g. the accretion disk; Dermer et al. 1992).If so, the BH γ-ray emission may be similar in spectral shape to the diffuse pionic emission tracing the (low-density) hydrogen distribution within R ⋆ .These two emission components would not be separated based on current LAT data, and the high proton (and ensuing secondary electron) density derived from the purely pionic fit may be an artifact due to spectral confusion.

Summary on M31
The broad-band radio-γ-ray non-thermal emission of M31 studied in this paper originates in the central (R ≤ R ⋆ ≃ 5.5 kpc) region.This radial scale is dictated by the Fermi-LAT sensitivity to source extension within which unresolved (i.e.point-like) emission has been detected by Xing et al. (2023).The radio (synchrotron and thermal free-free) emission from this region has been estimated starting from the all-disk integrated emission, using the radial profile of the 842 MHz emission line.
We suggest that the emission from M31's central region can be spectrally modelled as diffuse pionic with likely contributions from an unresolved population of ∼1000 Galactic-like pulsars plus the nuclear BH (M31 * ).In fact, the central 'hydrogen hole' may allow for emission from the enclosed pulsar population and M31 * to pass through.While the Galactic-like-pulsar emission used in our model differs in spectral shape from the pionic emission, the (putative) BH emission may be quite similar to it hence the two would not be separated based on current data.Thus, the proton component of the non-thermal interstellar gas, derived from a purely diffuse pionic fit to current LAT data, would be biased towards high values due to contamination from M31 * emission.

M33
Based on 11.4 yr of LAT data, Xi et al. (2020) detected a statistically significant emission from M33.Although it peaks on the northeast region of the galaxy, where the giant HII region NGC 604 16 is located, the LAT signal appears to originate from the whole galaxy.This is suggested by the fact that spatially extended templates based on IR maps of M33 (IRAS at 60µm, and Herschel Photoconductor Array Camera and Spectrometer [PACS] at 160µm) give almost equally good fits as the pointsource model centred at the (NGC 604-related) best-fit location 16 NGC 604 is a very active star-forming region that, within a cavity near the centre of the nebula, contains >200 massive blue stars -some with masses exceeding 120 M ⊙ and surface temperatures of ∼7 10 4 K).The cavity itself is blown by the collective stellar winds from these stars, whose strong ultraviolet radiation makes the surrounding gas fluorescent.

Pionic yields
The ambient proton density is n p ≃ 1.92 cm −3 (see Table 4 for densities in the different phases of the interstellar gas).The parameters describing the proton spectrum that matches the LAT data are listed in Table 6; in particular, the PL slope agrees with the values reported by Xi et al. (2020) for the the nominal 0.1-80 GeV LAT spectrum (Fig. 4).The respective relativistic proton energy density is u p = 0.46 eV cm −3 .The closely related π ± -decay secondary electron spectrum, N se (γ), is computed as described in Section 4. The parameters of its analytical fitting function, N se (γ), are given in Table 6.The corresponding synchrotron emissivity is specified in Eq.A.6.

Leptonic yields
To compute the electrons synchrotron emission we assume the relevant average total magnetic field in which primary and secondary electrons are immersed to be B = 6.5 µG (Tabatabaei et al. 2008).Once the secondary electron synchrotron yield has been computed, its low-frequency residuals from the data are modelled by attributing them to primary electrons.
To find the spectrum of the primary electrons, we adopt a phenomenological approach.We first try a PL spectrum but its resulting synchrotron emission (see Eq. A.4) grossly misses the data.We then use a smooth double-PL (2PL) spectrum: where N e0 is the normalisation, α and β are the low/high-energy slopes, γ b is the break energy, and γ max is the cutoff energy.Comparison with the steady-state electron spectrum in the presence of radiative losses (see Section 4) suggests that α = q inj − 1 and β = q inj + 1, where q inj is the injection spectral index.So the adopted 2PL spectrum is finally written as where N e0 , q inj , γ b are free parameters.The corresponding synchrotron emissivity is specified in Eq.A.5.The high measured value of the Hα emission measure (Deharveng & Pellet 1970) and the sub-GHz turnover of the radio spectrum (Israel et al. 1992) suggest free-free absorption of the radio emission at low frequencies by a cool ionised plasma.
The Galactic-absorption corrected EM in the direction to the observer (galaxy inclination is i = 57 • ) is EM= 64 cm −6 pc (Monnet 1971).By definition, the EM is an integral quantity, so a detailed knowledge of the diffuse ionised gas along the line of sight is not needed.In agreement with Israel et al. (1992), we find that, independent of the electron injection index, absorption of the non-thermal emission by a T e ∼ 500 K plasma is required to fit the ∼ < 100 MHz data.
At high frequencies, the ν −0.1 component represents diffuse thermal free-free emission.This emission may be gauged to the Hα flux 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 Hα flux may be used to predict the free-free emission.Following the latter approach (Klein et al. 2018), we modelled the thermal free-free emission from Eq. A.12 using T w e = 10 4 K for the warm plasma and F(Hα) = 2 10 −10 erg cm −2 s −1 (only 40% of the integrated Hα flux, 5 10 −10 erg cm −2 s −1 , is diffuse: Verley et al. 2009; also Hoopes & Walterboos 2000). 17The radio measurements and model, which includes primary and secondary synchrotron and thermal emission, are shown in Fig. 4 (top).
With the electron spectra determined, we can calculate the Compton and non-thermal-bremsstrahlung yields from electrons scattering off, respectively, CMB-EBL-GFL photons and thermal plasma nuclei using Eqs.(A7)-(A8).Although the shape of the γ-ray spectrum is distinctly pionic, the Comptonised starlight peaks near the pionic peak (ν ∼ 2 10 23 Hz), where it contributes 5% of the measured flux, and the non-thermal-bremsstrahlung peaks at ν ∼ 3 10 23 Hz, where it contributes 1%.

Results and discussion
While comparison of the LH SED model to the dataset is not formally a best-fit, the agreement is reasonable, given the independently set priors on the values of B and q 1 , and two of its features are well based.The first is the pionic fit to the LAT data: the γ-ray leptonic contribution is very minor, even with the appreciable level of he uncertainty in the electron spectral parameters (see below).The LAT data, however, show no spectral cutoff; the central value of the last LAT datapoint implies E max p = 100 GeV, 17 For consistency we note that computing the free-free emissivity using Z 2 n i = n e = 0.03 cm −3 (derived from the Hα emission measure; Tabatabaei et al. 2008) for the warm plasma in Eq.(A9) we obtain a spectral flux consistent with the one derived from the integral Hα flux in Eq.(A11).

N e0
q inj γ b γ max u p q p E max p N se0 q 1 q 2 γ b1 γ b2 η F(Hα) 10 −11 cm −3 10 4 10 4 eV cm −3 GeV 10 −13 cm −3 10 2 10 4 10 −10 erg/(cm 2 s) 3.5 +4.1 −2.0 2.25 ± 0.10 1.3 ± 0.1 3.7 0.46 2.25 100 2.8 0.2 2.5 4 3 2.5 5 Fig. 4. SED of M33.Radio spectrum (top): the emission model (solid curve), including primary and secondary synchrotron radiation (dotted curves of decreasing flux) as well as a thermal free-free component (dot-dashed curve), is overlaid with data from Table 2 (black dots).The empty squares represent the estimated radio flux after removing the contribution from ≤200 pc size sources (Tabatabaei et al. 2022).The total radio emission, resulting if the 'nominal' primary electron spectrum is used, is indicated by the dashed curve.Broad-band SED (bottom): data points are shown as dots, model predictions as curves.Emission components are plotted by the following line types: radio synchrotron, dotted; thermal free-free emission, dot-dashed; total radio, solid; Comptonised CMB, long-dashed; Comptonised starlight (EBL+FGL), short-dashed; non-thermal bremsstrahlung, dotted; and pionic, solid.For each type of non-thermal leptonic yield, the total, primary, and secondary emissions are denoted as curves of progressively lower flux.In the X-ray and γray range, on top of the separate components, a solid line indicates the total emission.The empty squares and dashed curve in the radio are as described in the top panel.
and even with the uncertainty reflected in the width of the energy band, it is still clear that E max p > 30 GeV.The second apparent feature is absorption of the synchrotron emission by warm thermal gas.It is clear that even for the hardest possible primary electron spectrum (q inj = 2) and in spite of the substantial scatter of data, the spectrum at ≤100 MHz can only be reproduced with absorption by a low-temperature (500 K) thermal gas, in accordance with a suggestion made by Israel et al. (1992), and similarly (for NGC 891) by Mulcahy et al. 2018.)The relatively flat γ-ray slope of M33 suggests that, unlike for M31, a substantial PSR contribution is unlikely.The PSR spectral cutoff (at several GeV) results in a spectral hump in the relevant energy range (see Fig. 1).This in contrast to the relatively flat PL-like pionic model that appears to naturally match the data.
The uncertainty in the total (HI, H 2 , HII) gas density affects the energy-loss coefficients used to compute the secondary electron spectrum (Paper I) and the ensuing normalisation and shape of the primary electron spectra and their yields.However, the electron density is effectively unconstrained because it depends on assuming a magnetic field value rather than using a value derived by modelling the non-thermal-X-ray emission as Comptonised CMB (e.g.Persic & Rephaeli 2019a,b) or the hard Xray-MeV as Comptonised infrared (IR) radiation.
A further source of uncertainty in modelling the primary electron spectrum comes from the large scatter in the radio data.A strong coupling among N e0 , q inj , γ b implies a broad range of acceptable values of these parameters (Table 4) because a nonlinear combination of q inj and γ b determines the effective synchrotron slope at ∼ > 1 GHz which, for an assumed B (derived assuming equipartition between the energy densities of the magnetic field and relativistic particles; Tabatabaei et al. 2008), determines the primary electron normalisation N e0 .That said, we may surmise q inj = 2.25 is realistic for M33 based on the following argument.Electron indices measured in Galactic supernova remnants peak at 2.15 ≤ q e ≤ 2.3 (Mandelartz & Becker Tjus 2015), and values in the same range are also found in the quietly star-forming Magellanic Clouds (Paper I).As for the proton indices, which closely reflect the injection index given the low energy losses suffered by protons via p-p interactions with the local medium, values in a similar range are found, 2.2 ≤ q p ≤ 2.6 (Mandelartz & Becker Tjus 2015;Paper I).Thus, a value q inj = 2.25 is quite reasonable for both protons and electrons.The corresponding energy density ratio between proton and (primary) electron, ≃22, is higher than its value for PL injection spectra, ≃17 (Persic & Rephaeli 2014), as a consequence of the higher energy losses of electrons than of protons.Finally, the electron cutoff energy γ max is constrained low enough as not to overproduce the flux at ≥2700 GHz and high enough not to violate the highest-frequency point (interpreted as an upper limit by Israel et al. 1992).
The data from Israel et al. (1992), while free from the major individual sources that could have been identified when the observations were made, are not free from the integrated source emission of fainter sources.Tabatabaei et al. (2022) pointed out that radio sources with sizes ≤200 pc constitute ∼36% (46%) of the total radio emission at 1.5 (6.3) GHz in the inner 4 kpc × 4 kpc disk of M33.Indeed, the nominal steady-state primary electron spectrum, N pe (γ) = k inj /(q inj − 1) γ −(q inj −1) /(b 0 + b 1 γ + b 2 γ 2 ), where Ṅinj (γ) = k inj γ −q in j is the electron spectral injection rate and the b j terms denote Coulomb, bremsstrahlung, and synchrotron and Compton energy losses, as in Eqs.(A23)-(A25), while neglecting diffusion and advection losses).While it ends up under-reproducing the Israel et al. (1992) data, it is still com-A&A proofs: manuscript no.47976_corr_corr patible with the suggestion from Tabatabaei et al. (2022).If so, the X-ray and γ-ray leptonic yields would be even lower, making the LAT-measured emission even more likely to be of pionic origin.
Less noisy radio data are needed to clarify several issues.At low frequencies more precise ∼ < 100 MHz data would help to constrain the cold-plasma temperature in the range of 250 K ∼ < T c e ∼ < 2000 K.This is essential because a more accurate determination of T c e is key to our knowledge of the cold gas in galactic halos (which are best studied in inclined galaxies by measuring the absorption of radio emission from their disks; e.g.Dettmar 1992).For example, in the framework of the present analysis, the high-resolution 74 MHz VLSS 18 data (Cohen et al. 2007) may be suitable to this aim after integrating the M33 map over the whole disk and suitably subtracting point sources in order to achieve consistency with the data in Israel et al. (1992).At higher frequencies, the data between 100 MHz and 2 GHz, which is free from plasma absorption and located in the bremsstrahlung energy-loss regime 19 , where the electron spectrum is N e (γ) ∝ γ −q inj , would constrain q inj and thereby allow for a more reliable determination of the primary electron spectrum.

Conclusion
In this paper, we present a spectral fitting of electron and proton radiative yields to radio and γ-ray measurements of the Local Group galaxies M31 and M33, namely: the central region ( ∼ < 5.5 kpc) for M31 and the whole disk for M33.For M31, our analysis suggests that a combination of diffuse pionic, pulsar, and nuclear-BH-related emissions could explain the LAT data.For M33, our analysis clearly indicates that the LAT-measured emission is mostly of pionic origin in agreement with Xi et al. 2020).While uncertainties in the 3D galaxy structure directly affect the overall normalisations of the particles spectral distribution functions, the model SEDs are essentially unaffected by variations of the galaxy structure parameters.
Future measurement of diffuse non-thermal X-ray emission from the galactic disks will substantially improve the SED fits and result in improved determinations of the particle spectral parameters.For example, measurements of non-thermal 1 keV flux that result (mostly) from Compton scattering of the electrons by the CMB fix the overall normalisation of the electron spectrum (Persic & Rephaeli 2019a,b).This follows from the fact that the good pionic fit to the γ-ray emission firmly determines the secondary electrum spectrum, which then effectively leads to the primary electron normalisation, N e0 .Modelling the radio spectrum as synchrotron would then provide the (volumeaveraged) magnetic field, B. A spectral determination of B would remove the need to assume particle-field equipartition or to use line-averaged fields, such as those derived by Faraday rotation studies of extragalactic polarised sources seen through the disk.
Given the importance of magnetic fields in modelling nonthermal SEDs, a comment on the values we have adopted in this study is in order.Mean values of the magnetic fields of both galaxies adopted from radio studies of M31 (Fletcher et al. 2004) and M33 (Tabatabaei et al. 2008) rely on the assumption of energy density equipartition between relativistic particles 18 VLA (Very Large Array) Low-frequency Sky Survey. 19For the assumed B = 6.5 µG (Tabatabaei et al. 2008) the characteristic synchrotron frequency, ν s = 3γ 2 B MHz implies a range of electron energies 2.3 10 3 ∼ < γ ∼ < 10 4 emitting between 100 MHz and 1 GHz.This range is well inside the bremsstrahlung-loss-dominated regime, 600 ≤ γ ≤ 25000.and magnetic field.In estimating relativistic electron and proton densities from radio measurements, the electron densities can be readily derived but assumptions must be made on the associated proton number (energy) densities that are also required to estimate the equipartition magnetic field, B eq .Fletcher et al. (2004) used u p /u e = 100 in the standard equipartition formula, B eq ∝ (1 + u p /u e ) 2/(5+q e ) (e.g.Persic & Rephaeli 2014), whereas Tabatabaei et al. (2008) used N p /N e = 100 in the 'modified' equipartition formula: B eq ∝ (1 + N p /N e ) 2/(5+q e ) (Beck & Krause 2005).Our SED model for M31 exhibits u p /u e ∼ 100, in agreement with Fletcher et al. (2004).As for M33, the quantity N p /N e can be estimated as follows.Assuming radiative losses only, particle number densities are nearly constant (when integrated over energy) over time, thus, the particle number ratio can be calculated at injection: for an electrically neutral non-thermal plasma: N p /N e = (m p /m e ) (q inj −1)/2 (Schlickeiser 2002).Our M33 SED model has q inj = 2.25 (Table 6), so N p /N e = 110.This numerical coincidences suggest that equipartition between relativistic particles and magnetic field (i.e.minimum energy in the nonthermal fluid ;Longair 1981;Govoni & Feretti 2004), is achieved in the central region of M31 and the disk of M33; these serve as a posteriori proof that our models use the B-values from Fletcher et al. (2004) and Tabatabaei et al. (2008) in a self-consistent way.Finally we note that because B eq ∝ (1 + ρ p /ρ e ) 2/(5+q e )≃0.27 , where ρ j denote energy (number) density for M31 (M33), the derived equipartition fields are relatively stable against variations of these quantities.
Concerning the magnetic field in M31, we should note that the approximately constant field B ≃ 7.5 µG derived by Fletcher et al. (2004) spans galactocentric radii of 6 kpc < R < 14 kpc.In the present analysis, we extrapolated that value inwards into the region R < R ⋆ = 5.5 kpc relevant to our study.The justification to do so on the fact that the average plasma density measured out to R ∼ 10 kpc is n i ∼ 0.05 (Trudoljubov et al. 2005), based on modelling the interaction between supernova remnants and the diffuse hot ionised plasma with a non-equilibrium collisional ionisation model.This value is close to that reported in Table 4 for the central region R ≤ 5.5 kpc.Since n e = Z 2 n i , the scaling relation B ∝ n 2/3 e (Rephaeli 1988) substantiates our decision to extrapolate the Fletcher et al. (2004) B-value into the central region of M31.
An independent estimate of u p may provide a consistency check of the preferred (mostly) hadronic models of the γ-ray data of M31 and M33.Combining the SN frequency with the residency timescale, which represents the characteristic proton residence time in the disk, τ res , (and assuming a nominal value of the fraction of total SN energy that is channeled to particle acceleration) enables an estimation for u p (e.g.Persic & Rephaeli 2010).Assuming no appreciable mechanical energy loss (e.g. in driving a galactic wind), the value of τ res is determined by the energy-loss timescale for p-p interactions, τ pp = (σ pp cn p ) −1 .At the most relevant energy range (σ pp ≃ 40 mb), the characteristic value is τ pp ∼ 3 10 7 n −1 p yr.During τ res , a number ν SN τ res of SN explode in a region of volume, V, thus depositing the kinetic energy of their ejecta, E ej = 10 51 erg per SN (Woosley & Weaver 1995), into the interstellar medium.Arguments based on the cosmic-ray energy budget in the Galaxy and SN statistics indicate that a fraction η ∼ 0.05 − 0.1 of this energy is available for accelerating particles (e.g., Higdon et al. 1998;Tatischeff 2008).Accordingly, we express the proton energy density as: where the various quantities are normalised to typical Galactic values.The results are as follows: (a) In the M31 disk region, we are considering (R ≤ R ⋆ = 5.45 kpc), the average gas density is n g ≃ 0.5 cm −3 hence τ res ∼ τ pp = 6 10 7 yr.The galaxy-wide SN rate is ν SN ≃ 0.01 yr −1 (type Ia: 0.006 ± 0.004, van den Bergh 1988; type II: 0.004, from Eq.14 of Persic & Rephaeli 2010 using L IR ).Assuming the local SN rate to scale as the stellar surface brightness, the SN rate integrated within R ⋆ (corresponding to 1.09 R d ) is 0.3 ν SN .Then, Eq.( 9) yields u p ≃ 7 eV cm −3 , which is compatible with the value suggested by our PSR+LH model with 10 3 PSR.We note that a similar value, ∼6 eV cm −3 , was deduced for the inner Galactic plane from γ-ray observations (Aharonian et al. 2006).(b) In M33, the galaxy-wide average gas density is n p ∼ 2 cm −3 ; hence, τ res ∼ τ pp ∼ 10 7 yr.The total (mostly core-collapse, Type II) SN rate is ν SN ≃ 0.003 yr −1 (van den Bergh & Tammann 1991;Pavlidou & Fields 2001).Then, Eq.( 9) yields u p ≃ 0.6 eV cm −3 (galaxy average), in fair agreement with the SED modelling result. 20 We conclude with a remark on the detectability of the p-p related neutrino emission via π ± decay from M31 and M33.With the proton spectrum fully determined by our fit to the measured γ-ray emission, we estimated (using Kelner et al.'s 2006 formalism, reported in Section A.2.3) that the broadly peaked (∼0.1-10GeV) π ± -decay neutrino fluxes from the two galaxies are ∼60 (M31) and ∼100 (M33) times lower than the flux from the Large Magellanic Cloud reported in Paper I. These fluxes are too low to make detections using any of the current and upcoming neutrino projects.where N e (γ) is a generic electron spectral density (in our case, see Eqs. 6, 7, A.29), ǫ and ǫ 1 are the incident and scattered photon energies, and r 0 = (e 2 /m e c 2 ) 2 is the electron classical radius.In Eq. (A.8), we set ǫ min = 0 and ǫ max /h = 10 12 , 10 13 , and 10 15 Hz for the CMB, IR, and optical components, respectively.
The probability distribution functions in Eqs.A.30, A.32, and A.36 are normalised, as per 1 0 f (x) dx = 1.The corresponding injection spectra for µ-decay muon and electron neutrinos are, respectively, (A.40) and .41)In Eqs.A.40 and A.41, the factor of 2 accounts for the contributions from both the π + and π − pions.

Fig. 1 .
Fig. 1. γ-ray spectrum of emission from the central source of M31, R ⋆ ≤ 0 • .4.LAT data points (Xing et al. 2023) are shown as dots, with errorbars indicating frequency and flux uncertainties, and model predictions (LH: top; PSR+LH: middle and bottom) are shown as curves.The leptonic components used in the full calculation include primary and secondary electron yields.

Fig. 2 .
Fig. 2. Radio spectrum of M31.The data points (from Battistelli et al. 2019) denote the whole-galaxy radio flux reduced to the estimated flux expected from within the R ⋆ radial region as described in the text.Also shown are: (top) Battistelli et al.'s (2019) fit (single electron population synchrotron: dotted curve; thermal free-free emission: dot-dashed curve), and (bottom) our LH emission model.The latter includes primary and secondary synchrotron radiation (dotted curves of increasing flux at their peaks) and a thermal free-free component (dot-dashed curve).

Fig. 3 .
Fig. 3. Broad-band lepto-hadronic SED of the M31 central source: data points are shown as dots (with errorbars indicating frequency and flux uncertainties), model predictions as curves (top).The galaxy-wide flux densities ofBattistelli et al. (2019) have been renormalised to the region encompassed by R ⋆ (i.e.∼15% of the total emission; see Section 4.1.2).Emission components are denoted by these line types: radio synchrotron, dotted; thermal free-free emission, dot-dashed; total radio, solid; Comptonised CMB, short-dashed; Comptonised starlight (EBL+FGL), long-dashed; non-thermal bremsstrahlung, dotted; pionic, solid.Leptonic emission includes that from primary and secondary electrons(secondary components are dominant).In the X-ray/γ-ray range, on top of the separate components, a thick solid line indicates the total emission.Same for a broad-band PSR+LH SED (middle).The 500 pulsars emission component (see text) is ahown as a dot-dashed line.Same as above but for 1000 pulsars (bottom).Note: in this case the secondary leptonic components are subdominant.
⋆ For derivations and references of these quantities, see Appendix B.