Exploring A Cosmic-Ray Origin of the Multi-wavelength Emission in M31

TL;DR

This paper investigates whether cosmic rays can explain the multi-wavelength emission observed in M31, proposing a multi-component model that fits the data but suggests higher cosmic-ray power than typical sources.

Contribution

It introduces a comprehensive cosmic-ray based model for M31's emission, considering primary electrons, secondary particles from protons, and their combined effects.

Findings

Multi-component cosmic-ray model fits M31's spectrum.

Model requires higher CR power than typical supernovae.

Both primary and secondary CR processes are relevant.

Abstract

A recent detection of spatially extended gamma-ray emission in the central region of the Andromeda galaxy (M31) has led to several possible explanations being put forth, including dark matter annihilation and millisecond pulsars. Another possibility is that the emission in M31 can be accounted for with a purely astrophysical cosmic-ray (CR) scenario. This scenario would lead to a rich multi-wavelength emission that can, in turn, be used to test it. Relativistic cosmic-ray electrons (CRe) in magnetic fields produce radio emission through synchrotron radiation, while X-rays and gamma rays are produced through inverse Compton scattering. Additionally, collisions of primary cosmic-ray protons (CRp) in the interstellar medium produce charged and neutral pions that then decay into secondary CRe (detectable through radiative processes) and gamma-rays. Here, we explore the viability of a CR origin for multi-wavelength emission in M31, taking into consideration three scenarios: a CR scenario dominated by primary CRe, one dominated by CRp and the resulting secondary CRe and gamma rays from neutral pion decay, and a final case in which both of these components exist simultaneously. We find that the multi-component model is the most promising, and is able to fit the multi-wavelength spectrum for a variety of astrophysical parameters consistent with previous studies of M31 and of cosmic-ray physics. However, the CR power injection implied by our models exceeds the estimated CR power injection from typical astrophysical cosmic-ray sources such as supernovae.

Figures

Tables

Table 1. Table 1: Free parameters and their values in our best-fit model for a power law with exponential cutoff primary electron source. For reference, we have included in the bottom row the experimental values for the parameters as described in the text. Reference values for the source normalization are discussed in the context of CR power output in section V.1.1 .

${\alpha }_e$	$E_{cut}$ (GeV)	$N_{C{R}_e}$ (GeV-1cm-3s-1)	$B_{\mu }$ ($\mu$G)	$N_{SL}$	${\chi }_{min}^2/d.o.f$
$2.14$	$514$	$1.1\times {10}^{-25}$	$1.7$	$1.02\times {10}^{-12}$	3.32/7
		Observational Values
$2.0-2.3$	${10}^3-{10}^4$	—	$5-10$	${10}^{-12}-{10}^{-11}$	—

Table 2. Table 2: Normalization factors for various magnetic field strengths in the case of CRe primaries, normalized to the radio emission.

$B_{\mu }$ ($\mu$G)	$N_{C{R}_e}$ (GeV-1cm-3s-1)
3	$3.39\times {10}^{-26}$
5	$1.28\times {10}^{-26}$
10	$3.85\times {10}^{-27}$

Table 3. Table 3: Parameters and their values in our best-fit model for a power law primary proton source.

${\alpha }_p$	$N_{C{R}_p}$ (GeV-1cm-3)	${\chi }_{min}^2/d.o.f$
$2.66$	$8.89\times {10}^{-8}$	2.25 / 5

Table 4. Table 4: Parameters and their values in a selection of well fitting models for the multi-component model, along with the corresponding χ m ​ i ​ n 2 subscript superscript 𝜒 2 𝑚 𝑖 𝑛 \chi^{2}_{min} . The best-fit model parameters are given in the top row and the corresponding SED is plotted in figure 7 .

${\alpha }_e$	${\alpha }_p$	$E_{cut}$ (GeV)	$N_{CRe}$ (GeV-1s-1cm-3)	$N_{CRp}$ (GeV-1cm-3)	$B_{\mu }$ ($\mu$G)	$N_{SL}$	${\chi }_{min}^2/d.o.f$
$2.04$	$2.75$	$1658$	$1.32\times {10}^{-26}$	$7.48\times {10}^{-8}$	$4.8$	$1.10\times {10}^{-12}$	$2.27/5$
$2.09$	$2.75$	$854$	$3.29\times {10}^{-26}$	$5.24\times {10}^{-8}$	${}_{}{}^{a}3$	$1.02\times {10}^{-12}$	$2.41/6$
$1.92$	$2.71$	$1633$	$6.17\times {10}^{-27}$	$7.75\times {10}^{-8}$	${}_{}{}^{a}7$	$1.12\times {10}^{-12}$	$2.42/6$
$1.71$	$2.67$	$1550$	$2.73\times {10}^{-27}$	$7.01\times {10}^{-8}$	${}_{}{}^{a}10$	$1.44\times {10}^{-12}$	$2.80/6$
$1.57$	$2.70$	$1353$	$1.77\times {10}^{-27}$	$7.13\times {10}^{-8}$	$12.4$	${}_{}{}^{a}5\times {10}^{-12}$	$3.25/6$

Equations

$$B(r)=B_{\mu }.$$

$$n_{CMB}(\nu )=\frac{8\pi {\nu }^2}{c^3}\frac{1}{e^{h\nu /kT}-1}.$$

$$n_{bulge}(r)\propto e^{-b_n\left[{\left(\frac{r}{r_b}\right)}^{1/n}-1\right]}.$$

$$n_{disk}(r)\propto e^{-\frac{r}{r_d}}.$$

$$n_{SL}(\nu ,r)=N_{SL}\frac{8\pi {\nu }^2/c^3}{e^{h\nu /kT}-1}\left[e^{-b_n\left[{\left(\frac{r}{r_b}\right)}^{1/n}-1\right]}+\frac{e^{-\frac{r}{r_d}}}{135}\right].$$

$$L=4\pi r^{2}c{\bar{u}}_{SL},$$

$$\begin{split}\frac{\partial}{\partial t}\frac{dn_{e^{\pm}}}{dE}=&\mbox{\boldmath$\nabla$}\left[D(E,\mathbf{r})\mbox{\boldmath$\nabla$}\frac{dn_{e^{\pm}}}{dE}\right]\\ &+\frac{\partial}{\partial E}\left[b(E,\mathbf{r})\frac{dn_{e^{\pm}}}{dE}\right]+Q(E,\mathbf{r}).\end{split}$$

$$\begin{array}{rl}b(E,\mathbf{r})& =b_{IC}(E,\mathbf{r})+b_{Synch.}(E,\mathbf{r})+b_{Coul.}(E)+b_{Brem.}(E)\\ & =b_{IC}^0{u}_{CMB}{E}^2+b_{IC}^0{u}_{SL}(r)E^2+b_{Synch.}^0{B}^2(r)E^2\\ & +b_{Coul.}^0{\bar{n}}_e\left(1+\log \left(\frac{E/m_e}{{\bar{n}}_e}\right)/75\right)\\ & +b_{Brem.}^0{\bar{n}}_e\left(\log \left(\frac{E/m_e}{{\bar{n}}_e}\right)+0.36\right).\end{array}$$

$$D(E)=D_0{E}^{\delta }$$

$$\begin{array}{rl}G(r,\Delta v)& =\frac{1}{\sqrt{4\pi \Delta v}}\sum _{n=-\infty }^{\infty }{\left(-1\right)}^n{\int }_0^{r_h}d{r}^{\prime }\frac{r^{\prime }}{r_n}\left(\frac{Q(E,r^{\prime })}{Q(E,r)}\right)\\ & \times \left[\exp \left(-\frac{(r^{\prime }-r_n{)}^2}{4\Delta v}\right)-\exp \left(-\frac{(r^{\prime }+r_n{)}^2}{4\Delta v}\right)\right],\end{array}$$

$$v(E)={\int }_E^{\infty }d\tilde{E}\frac{D(\tilde{E})}{b(\tilde{E})}.$$

$$\frac{d{n}_{e^{\pm }}}{dE}\left(E,r\right)=\frac{1}{b(E,r)}{\int }_E^{\infty }d{E}^{\prime }G\left(r,\Delta v\right)Q(E,r).$$

$$j_i(\nu ,r)=2{\int }_{m_e}^{\infty }dE{P}_i(\nu ,E,r)\frac{d{n}_{e^{\pm }}}{dE}\left(E,r\right)$$

$$S_i(\nu )=\frac{1}{4\pi d^2}\int dV{j}_i(\nu ,r)\approx \frac{1}{d^2}\int dr{r}^2{j}_i(\nu ,r)$$

$$P_{syn}\left(\nu ,E,r\right)={\int }_0^{\pi }d\theta \frac{\sin \theta }{2}2\pi \sqrt{3}{r}_0{m}_{e}c{\nu }_0\sin \theta F\left(\frac{x}{\sin \theta }\right),$$

$$x\equiv \frac{2\nu m_e^2}{3{\nu }_0{E}^2},$$

$$F(s)\equiv s{\int }_s^{\infty }d\zeta K_{5/3}\left(\zeta \right)\simeq 1.25s^{1/3}{e}^{-s}{\left[648+s^2\right]}^{1/12},$$

$$P_{IC}\left(E_{\gamma },E,r\right)=c{E}_{\gamma }\int dϵn\left(ϵ,r\right)\sigma \left(E_{\gamma },ϵ,E\right)$$

$$\sigma \left(E_{\gamma },ϵ,E\right)=\frac{3{\sigma }_T}{4ϵ{\gamma }^2}G\left(q,\Gamma \right),$$

$$G(q,\Gamma )=\left[2q\ln q+(1+2q)(1-q)+\frac{(2q{)}^2(1-q)}{2(1+\Gamma q)}\right],$$

$$\Gamma =\frac{4ϵ\gamma }{m_e{c}^2}=\frac{4{\gamma }^2ϵ}{E},q=\frac{E_{\gamma }}{\Gamma (E-E_{\gamma })}$$

$$F_{\gamma }=\frac{1}{d^2}\int dr{r}^2{E}^2{Q}_{\gamma }\left(E,r\right).$$

$$Q_{e^{\pm }}(E)=N_{C{R}_e}{\left(\frac{E}{\text{GeV}}\right)}^{-{\alpha }_e}{e}^{-E/E_{cut}}.$$

$$P_{C{R}_e}=\int dV{\int }_{m_e}^{\infty }dEE{Q}_e(E),$$

$$P_{SN,e}={\eta }_e\Gamma E_{SN}$$

$$\Gamma =SFR\times \frac{{\int }_{M_{min}}^{M_{max}}\psi (M)dM}{{\int }_{0.1M_{\odot }}^{100M_{\odot }}\psi (M)MdM}.$$

$$P_{PWN}=\frac{\eta W_0}{t_0{\left[1+\left(\frac{t_p}{t_0}\right)\right]}^2}$$

$$P_{PWN}=\frac{\eta W_0}{4t_0}\approx 8\times 10^{37}{\text{erg s}}^{-1}.$$

$$\begin{array}{rl}& {\pi }^0\to 2\gamma \\ & {\pi }^{\pm }\to {\mu }^{\pm }+{\nu }_{\mu }/{\bar{\nu }}_{\mu }\to e^{\pm }+{\nu }_e/{\bar{\nu }}_e\end{array}$$

$$n_{C{R}_p}(E)=N_{C{R}_p}{\left(\frac{E}{\text{GeV}}\right)}^{-{\alpha }_p}$$

Full text

††preprint: APS/123-QED

I Introduction

The study of gamma rays in galactic environments offers an intriguing probe of many physical phenomena including cosmic-ray production and transport, star formation rates, or new physics such as dark matter. The Andromeda galaxy (M31) is particularly enticing as a target of gamma-ray studies due its status as the nearest large spiral galaxy. M31 has been the focus of several previous gamma-ray searches Fichtel et al. (1975); Pollock et al. (1981); Sreekumar et al. (1994); Blom et al. (1999); Abdo et al. (2010); Ackermann et al. (2017); Karwin et al. (2019); Feng et al. (2018). Early observations Fichtel et al. (1975); Pollock et al. (1981); Sreekumar et al. (1994); Blom et al. (1999) were only able to place upper limits until the galaxy was first detected in gamma-rays using 2 years of Fermi-LAT data at $5.3\sigma$ significance, along with some evidence of a spatial extension at the $1.8\sigma$ confidence level Abdo et al. (2010). It has also been observed by high energy Cherenkov telescopes, though no detection has yet been made at energies above the TeV Albert (2018); Bird and VERITAS Collaboration (2015); Aharonian et al. (2003); Rubenzahl et al. (2017).

Recently, M31 was detected in gamma-rays by the Fermi telescope at a significance of nearly $10\sigma$ with a detection of spatially extended emission out to $\sim 5$ kpc at the $4\sigma$ significance level Ackermann et al. (2017). This emission resembles to some extent the well studied Galactic Center Excess (GCE) of gamma rays in the center of the Milky Way, and has led to comparisons in possible origins for the emission in the two galaxies. Proposed explanations for the GCE include signals of annihilating dark matter Abazajian et al. (2014); Calore et al. (2015); Daylan et al. (2016); Gordon and Macías (2013); Hooper and Goodenough (2011); Goodenough and Hooper (2009), an unresolved population of millisecond pulsars (MSP) Bartels et al. (2016); Brandt and Kocsis (2015); Eckner et al. (2018), or additional cosmic-ray sources Carlson et al. (2016); Gaggero et al. (2015); Cholis et al. (2015). Due to some similarities between both the two galaxies themselves and the observed emissions, it is natural to also consider whether these are viable explanations for the M31 detection. Although there is significant uncertainty in the dark matter density profile, the possibility of a dark matter signal in M31 has previously been studied for gamma-rays Fornengo et al. (2004); Albert (2018); Li et al. (2016), as well as other wavelengths McDaniel et al. (2018); Egorov and Pierpaoli (2013); Beck and Colafrancesco (2017); Chan (2016); Ng et al. (2019); Boyarsky et al. (2008); Watson et al. (2012). A brief argument is presented in the recent Fermi detection paper Ackermann et al. (2017) using the relative $J$-factors of the Galactic center and M31 to infer that the expected gamma-ray emission from dark matter annihilation in M31 is roughly a factor of $\sim 5$ below the observed emission. In a recent paper McDaniel et al. (2018), we studied the possibility of a dark matter origin of the M31 emission from a multi-wavelength perspective. We found that when assuming a dark matter only interpretation using favored GCE dark matter models, such models typically require annihilation cross-sections above current constraints, and have spectral shapes that are inconsistent with the M31 observations. However, this does not necessarily preclude these models contributing a subdominant component of the observed emission, with the shape of the spectrum predominantly determined by another emission source. Additionally,the dark matter particle models in that analysis that can reproduce the M31 gamma-ray emission also produce synchrotron emission that is in tension with observational radio data. There have also been efforts made to explore a millisecond pulsar (MSP) explanation for the M31 gamma-ray emission Eckner et al. (2018); Fragione et al. (2019). Ref. Fragione et al. (2019) studied MSPs originating from globular cluster disruption in the bulge of M31, whereas Ref. Eckner et al. (2018) considered MSPs formed in situ. The in situ model was found to fit the energetics and morphology of the excess well, however neither study could account for the full detected emission, with each providing only $\sim 1/4$ of the M31 observation.

In light of the lack of a definitive dark matter or exclusive unresolved MSP explanation, this work focuses on the scenario in which this emission can potentially be accounted for using a conventional astrophysical interpretation of cosmic ray induced emissions. In order to explore this possibility, we study the multi-wavelength emission in Andromeda to consider whether a cosmic ray explanation is consistent across the spectrum. We predominantly focus on the radio and gamma-ray emission, as these regimes provide the most insight given currently available observational data. The production of cosmic rays in astrophysical systems can lead to emissions at various wavelengths from radio, to X-ray and gamma-ray Blumenthal and Gould (1970); Longair (2011); Rybicki and Lightman (1979). When cosmic ray electrons (CRe) are injected into regions containing magnetic fields, they radiate synchrotron emission at radio wavelengths, whereas through inverse Compton scattering the electrons upscatter ambient photons, such as from the CMB or starlight, to X-ray and gamma-ray energies. Gamma rays can also be produced directly through the decay of neutral pions produced in cosmic-ray proton (CRp) collisions with the interstellar medium Dermer (1986); Mannheim and Schlickeiser (1994). However, as first pointed out in Pshirkov et. al. 2016 Pshirkov et al. (2016) and also found in the Ackermann et. al. (2017) Ackermann et al. (2017) analysis, the gamma-ray emission in M31 does not seem to be spatially correlated with neutral gas or regions of high star formation, which are typically the regions wherein CRp are expected to be produced and interact with the ISM to produce the $\pi^{0}$ gamma-rays. The magnitudes of these fluxes depend on the components of the astrophysical environment, such as magnetic field, target photon radiation field, and abundance of cosmic ray production mechanisms. The diffusion of the relativistic cosmic rays additionally plays a significant role in the expected signal, as cosmic rays diffuse and escape the system, thus suppressing the expected flux. The mechanism by which cosmic rays are produced and accelerated has been a topic of extensive interest. For galaxies, supernovae remnants (SNR) are considered to be the main sources of cosmic rays Ginzburg and Syrovatskii (1964); Berezhko (2014); Bell (2013); Thompson et al. (2012); Drury (2012). While SNRs are thought to provide the dominant contribution of cosmic rays in galaxies, other mechanisms such as pulsars and their nebula can also provide significant contributions to the total cosmic ray population Aharonian et al. (1995); Gendelev et al. (2010); Hooper et al. (2009). The expected power injection from these cosmic ray sources provides a benchmark to which we can compare the cosmic ray power necessary to produce the multi-wavelength emissions in M31.

This paper is organized as follows: in section II we detail the relevant physical modeling, including the magnetic field, interstellar radiation field, and diffusion model, as well as the solution to the diffusion equation. In section III we present the expressions for each radiative process under consideration, then in section IV we present the data used in the analysis. In section V we present the results of our analysis, in section VI we comment on X-ray diffuse emission constraints, and, finally, we conclude in section VII.

II Astrophysical Model of Andromeda

In order to calculate the secondary emission from the production of primary or secondary cosmic-ray electrons, we must first model the relevant astrophysical components of Andromeda. In particular, we require a description of the magnetic field model, which determines the synchrotron emission produced, as well as the inter-stellar radiation field (ISRF) that provides the target photon bath for the inverse Compton (IC) scattering. Also relevant to this analysis is a model of diffusion and radiative energy losses, since the relativistic electrons diffuse out and escape the system, while also losing energy through synchrotron emission, IC emission, Coulomb interactions, and bremsstrahlung. In the following section, we briefly describe the model adopted in this work.

II.1 Magnetic Field

The study of radio synchrotron emission in M31 requires some knowledge of the magnetic field in our region of interest, namely within the inner few kpc. In this region, the field structure is turbulent and complex, however estimates of the field strengths as determined by Faraday rotation measures of polarized radio emission are typically around $15\pm 3\mu$G for $r=0.2-0.4$ kpc, and $19\pm 3\mu$G for $r=0.8-1.0$ kpc Hoernes et al. (1998); Gießübel and Beck (2014). At larger radii in the disk of the galaxy, the magnetic field falls off to values of roughly $5\pm 1\mu$G Fletcher et al. (2004). In our previous study of Andromeda where we considered a dark matter origin of the multi-wavelength emission we selected a spatially dependent magnetic field based on these values, approximating some degree of spherical symmetry in the central region, while noting that that model would not be applicable at larger radii where a multi-dimensional field model would be more appropriate. In this analysis however, we treat the magnetic field strength as one of the free parameters in our fit, and so we adopt a simplified constant magnetic field where

[TABLE]

While this is helpful in that it reduces the parameters we need to fit in order to define our model, it comes at the cost of accuracy in capturing the complexity of the field or any spatial dependence. We can consider the constant field as an average over space, with a consequence of this being that we would expect the values of $B_{\mu}$ to be lower than the quoted central values, and at a roughly similar level to that of the disk.

II.2 Inter-stellar Radiation Field

Our inter-stellar radiation field model (ISRF) contains two components: (i) a CMB photon component and (ii) a starlight (SL) component. We have chosen to neglect the infrared (IR) component in the ISRF and instead focus only on the CMB and starlight components in order to simplify the ISRF modelling and fitting procedure, implicitly making the assumption that the starlight component will be the more significant contribution to the IC emission. For the CMB, we simply have a black-body spectrum at $T=2.73$ K and spatial homogeneity. Thus we have,

[TABLE]

For the starlight component, we approximate the spectrum as a black-body with temperature $T=3500$ K, a choice motivated by previous analysis of the ISRF in the Milky Way demonstrating this as a good approximation for starlight spectra. We additionally include a spatial dependence based on the starlight luminosity profile of M31 Courteau et al. (2011) incorporating a bulge component of the form:

[TABLE]

and a disk component:

[TABLE]

Combining these spatial components with the black-body spectral profile yields a starlight photon number density:

[TABLE]

The parameters $r_{b},\>r_{d},\>n,\>b_{n}$ are taken from Courteau et al. (2011) and the factor of $1/135$ in the disk component was chosen to recreate the bulge to disk luminosity ratio in Courteau et al. (2011). The factor $N_{SL}$ is a dimensionless normalization constant that is to be determined in the later sections as a free parameter in our fit. To get a sense for what value this parameter should be, we can consider the stellar luminosity of the inner region of M31. The stellar luminosity within a 1 kpc radius of M31 has previously been reported as $L=10^{9.9}L_{\odot}$ Groves et al. (2012); Draine et al. (2014). We can roughly estimate the luminosity as

[TABLE]

where the bar refers to a spatial average over the volume. Taking the radius to be $\sim 1$ kpc, a stellar luminosity of $L=10^{9.9}L_{\odot}$ corresponds to $\bar{u}_{SL}\approx 5$ eV cm*-3*, or $N_{SL}\approx 5\times 10^{-12}$. A reasonable, albeit somewhat large, range of values for the stellar energy density in the inner regions of galaxies is $\bar{u}_{SL}\approx 1-10$ eV cm*-3* Profumo and Ullio (2010); Porter et al. (2008), which roughly corresponds to a normalization fit range of $N_{SL}\in(10^{-12},10^{-11})$.

II.3 Solution to the Diffusion Equation

After being injected into the system, the CRe undergo both radiative losses and diffusion. Diffusion is particularly important on shorter distance scales, such as the few kpc scales considered in this work, and we have demonstrated in our previous M31 paper (see Fig. 6 from McDaniel et al. (2018)) that it significantly impacts the expected fluxes. The diffusion and radiative energy loss mechanisms of the CRe are accounted for in the diffusion equation:

[TABLE]

where we neglect convection and reacceleration effects which can be safely ignored for energies greater than a few GeV Delahaye et al. (2009, 2010). Particularly in a quiescent galaxy such as M31 with its low star-formation rate Ford et al. (2013); Rahmani et al. (2016), the effects of convection are expected to be less prominent than in galaxies with higher star-formation activity McCormick et al. (2013); Lacki et al. (2010); Veilleux et al. (2005) such as starbursts or even the Milky Way, and thus we treat diffusion as the dominant escape term. In equation 7, $\partial n_{e}/\partial E$ is the electron/positron equilibrium spectrum with units of GeV*-1* cm*-3*, $D(E,\mathbf{r})$ is the diffusion coefficient, $b(E,\mathbf{r})$ is the energy loss term, and $Q(E,\mathbf{r})$ is the CRe source term that we specify in later sections and has units of GeV*-1* s*-1* cm*-3*. In the energy loss term, we include contributions from synchrotron, IC, Coulomb, and bremsstrahlung processes, with the full expression given by:

[TABLE]

The $b^{0}$ coefficients in this expression have units GeV s*-1* with values $b_{syn}^{0}\simeq 0.0254$, $b_{IC}^{0}\simeq 0.76$, $b_{brem}^{0}\simeq 1.51$, and $b_{Coul}^{0}\simeq 6.13$ Longair (2011); Colafrancesco et al. (2006). The photon energy density for the CMB is $u_{CMB}=0.25$ eV cm*-3* and for the starlight photons can be computed from equation (5) to be $u_{SL}(r)=h\nu_{0}^{2}n_{SL}(\nu_{0},r)$, where $\nu_{0}$ is taken to be the peak frequency. Finally, $\bar{n}_{e}$ in equation (8) refers to the average thermal electron density and is taken to be $\bar{n}_{e}\approx 0.01$ cm*-3* Gaensler et al. (2008); Berkhuijsen et al. (2006); Berkhuijsen and Müller (2008); Caprioli (2011); Caprioli et al. (2010).

For the diffusion coefficient we assume a homogeneous power law of the form:

[TABLE]

with $\delta=1/3$ and $D_{0}=3\times 10^{28}$ cm2 s*-1* Strong et al. (2007); Vladimirov et al. (2012); Baltz and Edsjö (1999); Webber et al. (1992). The choices of these parameters are motivated by assuming that M31 has roughly similar diffusion properties to the Milky Way, with these values being determined by measurements of the stable (e.g. B/C) or unstable (e.g. Be10/Be9) secondary to primary ratios, and also supported by studies of the far-infrared - radio correlation in M31 and other galaxies that infer similar values Berkhuijsen et al. (2013); Murphy et al. (2006). Equation (7) can be solved analytically using the Green’s function method (see e.g. Colafrancesco et al. (2006); Ginzburg and Syrovatskii (1964)) and in the steady state case where the left-hand side of equation (7) is set to zero the appropriate Green’s function with free-escape boundary conditions is given by:

[TABLE]

where $r_{h}\approx 5$ kpc is the diffusion zone radius and the locations of the image charges used to implement the free-escape boundary condition are $r_{n}=\left(-1\right)^{n}r+2nr_{h}$. The value $\Delta v$ is defined as $\Delta v=v(E)-v(E^{\prime})$ with

[TABLE]

where we have approximated a spatially independent form of the energy loss term by taking a spatial average of $u_{SL}(r)$ and $B(r)$ in equation (8). In the above expression, $E^{\prime}$ represents the energy of the electron at the source, while $E$ is the interaction energy. The quantity $\sqrt{\Delta v}$ has units of distance, and represents the diffusion length scale of the particles. The final form of the electron equilibrium spectrum is then given by,

[TABLE]

Here we use the full spatially dependent form of the energy loss expression, rather than the homogeneous form used in equation (11).

III Multi-wavelength Emission

Once we have obtained the electron equilibrium spectrum $dn_{e^{\pm}}/dE$ by solving the diffusion equation, we can then proceed to calculate the emissivity $j_{i}$, by integrating the electron spectrum with the power for the given radiative process, namely the synchrotron radiation and IC scattering for our purposes. This gives

[TABLE]

where the factor of two accounts for electrons and positrons and $P_{i}$ is the power of a radiative process $i$ which we calculate in the following sections. From here, the flux density is given by the integral of the emissivity over volume,

[TABLE]

where $d$ is the distance to M31, taken to be $d=780$ kpc Stanek and Garnavich (1998). In this work, we make use of the publicly available RX-DMFIT tool McDaniel et al. (2017) to solve the differential diffusion equation and then to perform the various secondary emission calculations. Models used in this analysis can be obtained from the authors.

III.1 Synchrotron Power

In the presence of ambient magnetic fields, the relativistic CRe undergo synchrotron radiation, producing radio emission. The synchrotron power for a frequency $\nu$ averaged over all directions is Storm et al. (2016); Longair (2011)

[TABLE]

where $r_{0}=e^{2}/(m_{e}c^{2})$ is the classical electron radius, $\theta$ is the pitch angle, and $\nu_{0}=eB/(2\pi m_{e}c)$ is the non-relativistic gyrofrequency. The $x$ and $F$ terms are defined as,

[TABLE]

[TABLE]

where $K_{5/3}$ is the Bessel function of order 5/3.

III.2 Inverse Compton Power

With the photon number density $n(\epsilon,r)=n_{CMB}(\epsilon)+n_{SL}(\epsilon,r)$, and the IC scattering cross-section $\sigma\left(E_{\gamma},\epsilon,E\right)$, the IC power is

[TABLE]

where $\epsilon$ is the energy of the target photons, $E$ is the energy of the relativistic electrons and positrons, and $E_{\gamma}$ is the energy of the photons after scattering (note that $E_{\gamma}=h\nu$ for observing frequency $\nu$ in equation (13)). The scattering cross-section, $\sigma\left(E_{\gamma},\epsilon,E\right)$, is given by the Klein-Nishina formula:

[TABLE]

where $\sigma_{T}$ is the Thomson cross-section and $G(q,\Gamma)$ is given by Blumenthal and Gould (1970):

[TABLE]

where,

[TABLE]

The kinematics of inverse Compton scattering set the range of $q$ to be $1/\left(4\gamma^{2}\right)\leq q\leq 1$ Colafrancesco et al. (2006); Blumenthal and Gould (1970); Rybicki and Lightman (1979).

III.3 Gamma-ray Flux

In addition to gamma rays produced from IC scattering, we also consider gamma rays resulting from the decay of neutral pions produced in cosmic-ray proton collisions. When the pions are produced they decay rapidly within a time span of $\sim 10^{-16}$ s. The gamma rays do not experience diffusion or radiative loss effects, and thus we do not need to consider equation (7). Instead, for a $\pi^{0}$ gamma-ray source injection $Q_{\gamma}$ (in units of GeV*-1* cm*-3* s*-1*), the flux is simply given by integrating over the volume of the source Colafrancesco et al. (2006); Regis and Ullio (2008); McDaniel et al. (2018):

[TABLE]

IV Gamma-ray and Radio Data

The gamma-ray data points are taken from the analysis performed in Ackermann et. al. (2017) Ackermann et al. (2017), where they used 88 months of PASS 8 Fermi data collected between August 4, 2008, and December 1, 2015. Reconstructed events within an energy range of 0.1-100 GeV were considered as well as reconstructed directions within a $14^{\circ}\times 14^{\circ}$ region centered at $(\alpha,\delta)=(10^{\circ}.6847,41^{\circ}.2687)$. SOURCE class events were used excluding those with zenith angle greater than $90^{\circ}$ or rocking angle greater than $52^{\circ}$. The resulting detected emission found in this study was concentrated within the inner 5 kpc, motivating this as the choice of region of interest in our calculations of the gamma-ray emission resulting from IC scattering and pion decay.

Radio observations of M31 have predominantly focused on regions of large radii out to about $\sim 16$ kpc with a particular emphasis on the star-forming 10 kpc “ring” Berkhuijsen et al. (2003); Gießübel et al. (2013); Beck et al. (1980, 1998); Beck and Graeve (1982), or alternatively on the central regions within $r\sim 1$ kpc Walterbos and Graeve (1985); Berkhuijsen et al. (1983). For our purposes, the available data in the 1 kpc region are most useful, as they allows us to focus our analysis on the inner region and make better comparisons between the 1 kpc ROI for radio emission and the 5 kpc ROI for gamma-ray emission. In each of the spectral energy distributions in the following sections, the synchrotron emission is calculated with $r=1$ kpc while the IC and gamma-ray emission are calculated with $r=5$ kpc.

V Results

In our analysis, the multi-wavelength emission in M31 is assumed to be due to the presence of cosmic-rays. We consider here two cosmic-rays production mechanisms: the first is primary production of CRe following a power law with exponential cut-off source injection, which then radiate synchrotron and IC emission. The second source of cosmic-rays we consider is primary production of cosmic-ray protons obeying a power law. The hadronic inelastic interactions of the CRp produce neutral pions that decay into gamma rays, as well as charged pions that decay into secondary CRe which then produce synchrotron and IC emission. Finally, we consider the scenario in which both of these sources provide comparable contributions to the overall cosmic-ray abundance in what we refer to as our “multi-component” model. We then examine to what extent each of these three scenarios can be responsible for the multi-wavelength emission in M31. Of the three cases mentioned, the multi-component model appears the most convincing, while the primary-only and secondary-only models do not easily reproduce the emission in M31 within the range of realistic parameter space.

V.1 Emission from primary cosmic ray electrons

We now define the source term of equation (7) by considering the case in which the cosmic-ray population is dominated by primary electron production obeying a power law with an exponential cutoff:

[TABLE]

In this section, $\alpha_{CR},\>N_{CR_{e}},\>E_{cut}$ in equation (23) along with $B_{\mu}$ and $N_{SL}$ are taken to be free parameters that we adjust to fit the observed radio and gamma-ray spectra. Previous studies of cosmic-ray origins can provide some guidance as to reasonable values for these parameters. For example, values of $\alpha_{e}\sim 2.0-2.3$ have been found to be consistent with production of cosmic rays in supernovae (SNe), as well as suggesting $E_{cut}$ values on the order of a few TeV Delahaye et al. (2010); Di Mauro et al. (2014); Manconi et al. (2017); Fang et al. (2017); Di Bernardo et al. (2013). The normalization $N_{CR_{e}}$ however is poorly constrained, and in section V.1.1 we compare the fit values of $N_{CR_{e}}$ with the corresponding SNe power output.

In table 1 we list the results for the best-fitting model and plot the SED in figure 1. Throughout this analysis we fit the parameters by minimizing the standard $\chi^{2}$ metric using a Nelder-Mead simplex algorithm Nelder and Mead (1965); Press (2007) to two decimal tolerance and employ a penalty function to enforce the fit range constraints. Additionally, the expansion, contraction, and shrink parameters are determined using a dimension dependent implementation Gao and Han (2010). For the source term parameters $\alpha_{e}$ and $E_{cut}$, we see a general agreement with expectations as described above, albeit with a cutoff energy somewhat lower than the TeV level. The starlight component is also suppressedto the lower end of the allowed range, with a normalization factor $N_{SL}=1.02\times 10^{-12}$ which is a factor of $\sim 5$ lower than the value derived in section II.2. We find a magnetic field value of $B_{\mu}=1.7\mu$G which is relatively small in comparison to those discussed in section II.1. While the actual structure of the field would involve higher central values with some spatially dependent fall-off, the average strength of the field over the space can be expected to take a smaller value. However, even with this in mind the magnetic field value is particularly low and likely not representative of the field strength within the inner regions of M31, especially the inner $\sim 1$ kpc where the synchrotron emission is calculated. Thus, we instead seek a configuration that allows for a higher magnetic field value.

One way in which we can potentially achieve a higher magnetic field is to take into consideration the case where the radio emission is due to synchrotron radiation from cosmic-ray electrons, but the IC emission is not sufficient to recreate the Fermi observations, and remain agnostic as to the source of the gamma-ray emission. To do this, we increase the strength of the magnetic field and change $N_{CR_{e}}$ to reproduce the radio emission. In figure 2 we show these fluxes for a few values of the magnetic field, and list the normalization factors in table 2. In this approach, we are essentially assuming that for reasonable magnetic field values the radio synchrotron emission in M31 can be produced predominantly by primary cosmic-ray electrons, while the source of the gamma ray emission remains unaccounted for. In later sections we use this approach in conjunction with comic-ray secondaries to account for the full spectrum of emission.

V.1.1 Cosmic-ray Electron Power

In order to place our value for the source term normalization $N_{CR_{e}}$ into a physical context, we can compare the total power output injected into cosmic-ray electrons with physical processes such as supernovae explosions. The power injected into the CRe for the source term of equation (23) is given by,

[TABLE]

where $V$ is the diffusion volume. Meanwhile, the power injected into cosmic-ray electrons from supernovae is given by the expression,

[TABLE]

where $\Gamma$ is the supernova (SN) rate, $E_{SN}$ is the total energy released in the SN explosion and $\eta_{e}$ is the efficiency of the SN energy transferred to the CRe. The SNe rate can be estimated from the observed star-formation rate (SFR), which in the case of M31 has been measured to be $\sim 0.2-0.4$ M*⊙* yr*-1* Ford et al. (2013); Rahmani et al. (2016). Adopting a value of $SFR=0.25$ M*⊙* yr*-1* Ford et al. (2013), the rate is then given by Horiuchi et al. (2011); Botticella et al. (2012)

[TABLE]

We use the Salpeter initial mass function (IMF) Salpeter (1955) defined over the main-sequence mass range of $0.1-100\>M_{\odot}$ wherein $\psi(M)\propto M^{-2.35}$ and $\psi(M)dM$ gives the number of stars in the mass range $M+dM$. We take $M_{min}=8M_{\odot}$ and $M_{max}=40M_{\odot}$ in line with canonical CC SNe parameters Horiuchi et al. (2011). This yields a SNe rate of $0.17$ per century, and the total energy output of for one supernova explosion is $E_{SN}\sim 10^{51}$ erg. While the efficiency at which energy is imparted to electrons during SN explosions is not well constrained, several estimates suggest values of $\eta_{e}=10^{-5}-10^{-3}$ Tatischeff (2009); Delahaye et al. (2010). Putting these together, we obtain a lower limit on the power injected into CRe in SNe explosions to be $P_{SN,\>e}\approx 5.1\times 10^{35}$ erg s*-1*, and an upper limit of $P_{SN,\>e}\approx 5.1\times 10^{37}$ erg s*-1*. In figure 3 we show the power injected into CRe implied by our best-fit model while increasing the magnetic field and normalizing to the radio data (as in figure 2). We compare this with the estimated range of SNe power output for CRe and see that the necessary normalization to fit the radio data produces a power requirement that is substantially greater than the estimated SNe power budget for the lower magnetic fields, including at the best-fit value when also fitting the gamma-ray data at $B_{\mu}=1.7\mu$G. Although the SNe power calculations involve a great deal of uncertainty, it is unlikely that the uncertainty is so great that it can be reconciled with the power output implied by our parameter model. Potential other cosmic-ray acceleration mechanisms such as PWNe could also contribute to the power total, however we can briefly demonstrate that this contribution is not enough to overcome the difference. For the case of pulsars, the relevant quantity is the spin down luminosity which can be expressed as

[TABLE]

where $\eta$ is the injection efficiency, $W_{0}$ is the pulsar energy output, $t_{0}$ is the typical pulsar decay timescale and $t_{p}$ is the pulsar lifetime Delahaye et al. (2010); Malyshev et al. (2009); Della Torre et al. (2015); Linden and Profumo (2013). If we take as fiducial values, $\eta=0.1$, $W_{0}=10^{50}$ erg, $t_{0}=1$ kyr and assume $t_{p}\approx t_{0}$ we obtain a power contribution from pulsars of

[TABLE]

While this value suggests that PWNe can contribute a significant amount of the CRe power, the estimate here is not sufficient to account for the necessary CRe power of our best-fit models, and thus does not have a significant impact on the results shown in figure 3.

V.2 Emissions from cosmic rays of hadronic origin

We next consider the scenario in which the dominant contribution to the cosmic-ray population is in the form of primary cosmic-ray protons. Inelastic interactions between the CRp and the interstellar medium produce neutral and charged pions. The neutral pions decay into gamma-rays, while the charged pions decay into muons and neutrinos, which in turn decay into neutrinos and secondary CRe. This can be summarized as:

[TABLE]

For the most common astrophysical model of the CRp distribution we assume a simple power law:

[TABLE]

with $N_{CR_{p}}$ in units of GeV*-1* cm*-3*. The resulting source terms for the gamma rays and cosmic-ray electrons have been previously calculated for this choice of CRp source distribution Pfrommer and Enßlin (2004); Mannheim and Schlickeiser (1994); Schlickeiser (2002). Following Pfrommer and Enßlin (2004) for the gamma-ray source term from $\pi^{0}$ decay yields the expression:

[TABLE]

with, $\alpha_{\gamma}=4/3\left(\alpha_{p}-0.5\right)$, and the source for $e^{\pm}$ from the charged pion decay is given by

[TABLE]

as described in Dolag and Enßlin (2000). Here, $n_{N}(r)$ is the nucleon number density, which we take to be proportional to the thermal electron number density with $n_{N}(r)=\frac{1}{1-\frac{1}{2}X_{H}}n_{e}(r)$ where $X_{H}=0.24$ is the the primordial $\prescript{4}{}{\text{He}}$ mass fraction. The thermal electron density $n_{e}(r)$ can be modeled as a beta-fit of the form

[TABLE]

with $\beta=0.49$ and $r_{c}=54^{\prime\prime}$ Liu et al. (2010) and assuming $n_{e,0}\sim 0.1$ cm*-3* Caprioli (2011); Caprioli et al. (2010). The neutral pion mass is $m_{\pi^{0}}=135$ GeV, $\xi$ gives the pion multiplicity taken to be $\xi=2$ for $\pi^{0}$, and $\sigma_{pp}=32$ mbarn is the proton collision cross-section. The shape parameter $\delta$ is given by $\delta=0.14\alpha^{-1.6}_{\gamma}+0.44$. For this case, when fitting to both the radio and gamma-ray data using the same free parameters as the previous section (but with $E_{cut}$ excluded and $N_{CR_{e}},\alpha_{e}$ replaced with $N_{CR_{p}},\alpha_{p}$) we are unable to find a reasonable fit, due to the significant difference between the index required to fit the synchrotron emission to the radio simultaneously as the $\pi^{0}$ decay to the gamma-ray emission. Additionally, if we ignore the contribution from $\pi^{0}$ gamma rays and assume that synchrotron and IC emission from secondary electrons are dominantly responsible for the observed radio and gamma-ray emission we similarly do not find a good fit to the data.

Instead, we determine $\alpha_{p}$ and $N_{CR_{p}}$ by only fitting the $\pi^{0}$ gamma rays to the Fermi data, while leaving the other parameters to be determined separately. With the only contribution to the fit being from the $\pi^{0}$ gamma ray contribution, we are find a best fit with $\alpha_{p}=2.66$ and $N_{CR_{p}}=8.89\times 10^{-8}$ GeV*-1* cm*-3*, also listed in table 3. The gamma-ray spectrum is shown in figure 4 along with a few other values of $\alpha_{p}$, normalized appropriately. In the selection of the remaining parameters that need to be determined (i.e. $B_{\mu}$ and $N_{SL}$) we are mainly constrained by the requirement that we are consistent with the field values in section II.1 while not overproducing the radio emission, and $N_{SL}$ does not result in IC emission that significantly impacts the spectrum from pion decay gamma rays in the Fermi data energy range. In figure 5 we show the result of this procedure with various values for the magnetic fields, and for simplicity a single starlight normalization $N_{SL}=5\times 10^{-12}$ in accordance with the discussion of section II.2. In this figure, two things are evident: first, the gamma-ray emission provides a good fit to the Fermi data; and second, the spectral index required for this fit results in a significant mismatch to the radio data regardless of normalization or field strength. The index of the CRp distribution obtained is in agreement with other studies that suggest $\alpha_{p}\sim 2.5-2.75$ Blasi and Amato (2012); Huang (2003); Paglione et al. (1996).

V.2.1 Cosmic-ray Proton Power

As we did in section V.1.1, we can again compare the power injected into CRp as implied by our fit parameters to the energy budget of SNe produced CRp. The power injection from SNe to CRp is of the same form as the CRe;

[TABLE]

where the only difference is in the value of the power injection efficiency, $\eta_{p}$. While Tatischeff (2009) inferred a value of $\eta_{p}\sim 10^{-5}-10^{-4}$, others have adopted higher values of $\eta_{p}\sim 10^{-3}$ Cassam-Chenaï et al. (2007). Additionally, gamma ray observations suggest that up to $3-30\%$ of the SN kinetic energy can be imparted into the cosmic-ray protons Ackermann et al. (2013); Thompson et al. (2012). We therefore have quite a large range of possible values, finding $P_{SN,\>p}\approx 5.1\times 10^{35}$ erg s*-1* for our lower bound and $P_{SN,\>p}\approx 1.53\times 10^{40}$ erg s*-1* as an upper bound.

To calculate the implied CRp power from our models, we take into account the diffusive properties of the CRp source distribution $n_{CR_{p}}$. Noting that for the heavier cosmic-ray protons the radiative energy losses of equation (8) are unimportant, we can consider only the propagation of the CRp by diffusion. The steady-state distribution of cosmic-ray protons has a characteristic diffusion timescale of $t_{D}(E)\approx r_{h}^{2}/D(E)$ Taylor et al. (2014); Pavlidou and Fields (2001); Gaggero and Valli (2018); Longair (2011); Blasi and Amato (2019), which gives us an injection source term $n_{CR_{p}}(E)/t_{D}(E)$. We then have for the power injected into CRp:

[TABLE]

In figure 6 we show the contours of the implied power injected into CRp. In this case the power determined by the fit parameters still exceeds the estimated SNe power injection. This discrepancy between the SNe estimates and our calculated power is not as extreme as in the primary $P_{CR_{e}}$ scenario for lower magnetic fields, though for higher field values, the CRe power is just over an order of magnitude greater than the upper SNe power, as opposed to the almost two order of magnitude difference for the CRp seen here.

V.3 Multi-component cosmic ray source model

In the previous sections we were working under the assumption that the cosmic-ray source was dominated by either primary production of CRe or hadronic production of secondary CRe and $\pi^{0}$ gamma rays. However, another possible scenario would be where both of these cosmic-ray production mechanisms are incorporated. We can therefore consider a multi-component model that includes contributions of both the primary source as well as the hadronically produced sources. For the gamma-ray source term, the only contribution is from the decay of pions produced in inelastic hadronic collisions as described by (31). The electron source term for the multi-component model is the sum of the source terms in equations (23) and (32):

[TABLE]

with, $\alpha_{\gamma}=4/3\left(\alpha_{p}-0.5\right)$. The best-fit results are listed in table 4 along with a selection of parameter sets with fixed magnetic fields or fixed $N_{SL}$. The SED for the best fit is shown in figure 7.

For the injection indices, we obtain values of $\alpha_{e}=2.04$ and $\alpha_{p}=2.75$, which are both within the ranges discussed in sections V.1 and V.2, while the normalization factors $N_{CRe}$ and $N_{CRp}$ do not deviate significantly from the values found in those sections. The cutoff energy $E_{cut}=1658$ GeV is in line with the $\sim$ TeV level used in previous cosmic-ray studies Delahaye et al. (2010); Di Mauro et al. (2014); Manconi et al. (2017); Fang et al. (2017); Di Bernardo et al. (2013), but higher than in the primary-only case. The magnetic field is also higher here than in the primary only case and is in good agreement with M31 magnetic field estimates. The similarity between the parameters of the multi-component model and the primary-only or secondary-only models is reflected in that for the multi-component model each of the two components (primary and secondary) have separate regimes of dominance. That is to say, the radio is predominantly due to the primary CRe whereas the gamma-rays are mainly due to the neutral pion decay gamma rays. This resolves the discrepancy in the model with purely hadronically produced CRe between the spectrum of radio data and the predicted synchrotron emission. In addition to the best-fit model, we also list in table 4 models in which we hold the magnetic field fixed and fit for the remaining free parameters. We do this as well where we instead hold $N_{SL}$ fixed to the value discussed in II.2 and fit the remaining parameters. In either case, the parameter held fixed is denoted in table 4 by the ‘$a$’ superscript. The spectra for each model in 4 are plotted in figure 8. With different field strengths we are still able to find good fits to the data, with only very slight changes to the $\chi_{min}^{2}$. This suggests that in the multi-component model there is no issue with a suppressed magnetic field as in the primary-only case. Again we see that the starlight normalization is highly suppressed and the IC emission is heavily dominated by the CMB component. Since this appears to be a fairly extreme scenario for the central region of the galaxy, we try to achieve a more reasonable value by holding the normalization fixed at $N_{SL}=5\times 10^{-12}$ as derived in section II.2 and fit the remaining parameters. We are still able to achieve a good fit, however it requires a relatively higher magnetic field of $B_{\mu}=12.4\mu$G in order to suppress the stellar IC component, as well as a low injection index of $\alpha_{e}=1.57$.

We once again compare the power injection into CRe and CRp implied by the parameters of our fit with the estimated SN injected power. Noting that the source term parameters for the CRe and CRp do not deviate significantly from the values found in section V.1 and V.2, similar results in this comparison can be expected here. In fact, that is essentially what we see in figure 9, wherein we show the implied CR power from our models for the various magnetic field values compared with the SN power injection estimates of sections V.1.1 and V.2.1. We see that the implied CRe power injection decreases for models with higher magnetic field (cf. figure 3), while the CRp injection remains relatively constant with some slight increase due to suppression of the primary CRe induced IC emission. However, neither are within their respective ranges for the SN source power. Although there is a discrepancy between the implied power injection of our cosmic-ray parameter sets and the estimated supernovae contribution, the great deal of uncertainty in the SNe power estimates makes it difficult to make concrete statements on the viability of these models on this basis alone.

VI Diffuse X-ray Emission in M31

While the focus of this analysis has been on the radio and gamma-ray emission, X-ray emission in M31 provides another potential avenue to study. However, several observations of the X-ray emission in the bulge of M31 have detected the presence of an unresolved diffuse component, using data from ROSAT Supper et al. (2001, 1997), XMM-Newton Shirey et al. (2001), and CHANDRA Li and Wang (2007). In each of these studies, diffuse X-ray flux in the inner $\sim$ 1 kpc of M31 is observed at a flux level of falling roughly between $\sim 3-5\times 10^{-12}$ erg cm*-2* s*-1*, and can likely be attributed to the presence of thermal hot gas and unresolved X-ray point sources. We note that this observed X-ray emission within 1 kpc has a higher flux than the X-rays produced in any of our cosmic-ray models from the previous sections, even despite the computed X-ray emission being within a 5 kpc radius. We thus conclude that for these cosmic ray models and our astrophysical setup the X-ray emission in M31 does not provide particularly useful information due to the bright diffuse emission in the bulge of M31 being considerably brighter than what we would obtain in our models.

VII Conclusion

We have examined the possibility of a cosmic-ray origin for the multi-wavelength emission in the Andromeda galaxy, specifically addressing the origin of the recently detected extended gamma-ray emission. We considered three models for the production of cosmic rays. First, we considered a primary injection of CRe obeying a power law with an exponential cutoff, then considered production of secondary CRe and gamma rays produced from interactions of a power law distribution of primary CRp. Finally, we looked at a multi-component model that incorporates both of these cosmic-ray sources. We then fit the synchrotron and IC fluxes arising due to the presence of the primary and secondary CRe, as well as the gamma ray emission from neutral pion decay, to available radio data and a recent Fermi gamma-ray detection in M31.

For the primary CRe scenario, we find best fit parameters for the injection spectrum $\alpha_{e}=2.14$ and cutoff energy $E_{cut}=514$ GeV. The injection index is consistent with expected values for CRe sources such as SNR. The cutoff energy is slightly lower than expected, however not wholly inconsistent with expected values on the order of TeV. The magnetic field value of $B_{\mu}=1.7\>\mu$G and the starlight normalization are both suppressed in the fit. We also considered higher magnetic fields and renormalized the synchrotron emission to match the radio data. This suppresses the IC gamma-ray emission, requiring that we account for the Fermi data separately which was done in the multi-component model. We then compared the power injection into CRe implied by our model with the expected range of power injection due to SNe. We saw that even by increasing the magnetic field in order to lower the normalization constant $N_{CRe}$, the power injection implied by our models was well above the expected output from astrophysical sources such as PWNe and SNe.

In the case where we considered contributions from only secondary cosmic rays of hadronic origin, we were unable to find a good fit to both the radio and gamma-ray data simultaneously. Rather, we assumed that the gamma rays were purely from the neutral pion decay and found a CRp distribution index of $\alpha_{p}=2.66$, consistent with previous results for $\pi^{0}$ gamma-ray studies, along with a CRp distribution coefficient of $N_{CR_{p}}=8.89\times 10^{-8}$ GeV*-1* cm*-3*. With this arrangement we then manually selected the magnetic field and starlight energy density, and found that for a variety of field strengths the calculated flux remains below the radio data, and even for a higher selected value of $N_{SL}=5\times 10^{-12}$ there was no conflict between the IC emission and the gamma-ray data. We again compared the power injection into CRp from SNe with the implied power output of our models, and found that the CRp injection is also greater than the estimated SNe output.

Finally, we consider a combined “multi-component model” that incorporates the contributions from both the primary CRe as well as the secondary CRe of hadronic origin. Although here the power budget concerns remained due to minimal variation in the best fit normalization constants, this scenario gives the best overall fit to the data, while still providing similar parameter values as in the primary-only and secondary-only cases. We found the best-fit $\alpha_{e}=2.04$ and $\alpha_{p}=2.75$, both similar to the values discussed in sections V.1 and V.2 respectively, while the best-fit magnetic field was found to be $B_{\mu}=4.8\>\mu$G and $E_{cut}$ was $1658$ GeV. Additionally, the multi-component model offers a large degree of flexibility in the parameter choices, as evidenced by good fits for a range of multiple magnetic fields values and $E_{cut}$ on the order of a few TeV, as well as for higher $N_{SL}$ values in accordance with the observed stellar luminosity in the central region of M31. In our final power comparison we saw similar results as in the primary-only and secondary-only scenarios. That is, both the implied CRe and CRp power in our models were greater than the estimated power output from astrophysical sources, and this held at a wide range of magnetic fields values. This suggests that although the spectra can be fit well with a multi-component model, the input power needed for the cosmic-ray sources is consistently more than an order of magnitude above what is expected from supernova as galactic cosmic-ray accelerators. Furthermore, as mentioned in the introduction and discussed more thoroughly in the original detection paper Ackermann et al. (2017) along with earlier and subsequent Fermi M31 studies Pshirkov et al. (2016); Eckner et al. (2018); Fragione et al. (2019), the gamma-ray emission does not appear to correlate with star-formation or gas-rich regions. CRp produced at larger radii that then diffuse into the emission region may contribute to the observed signal, although this does not address the lack of gas for the CRp to interact with in the interior regions of the galaxy. Another possibility is that the CRp are remnants from a previous period of higher star-formation. However, the stellar population of the bulge is dominated by stars with ages $\gtrsim 4-12$ Gyr Olsen et al. (2006); Saglia et al. (2010); Dong et al. (2018), compared to a CRp escape time of $\sim 10-100$ Myr, which suggests that the majority of CRp would have likely left the system in the time since this higher star-formation activity in M31. This morphological point along with the power discrepancy combine to disfavour a purely CR explanation, particularly one that relies on $\pi^{0}$ gamma-rays from CRp to explain the observed gamma-ray emission.

Acknowledgements.

This study is based on work supported by the National Science Foundation under Grant No. 1517545. A.M. is supported by a Department of Education GAANN fellowship. S.P. is partly supported by the US Department of Energy, grant number DE-SC0010107.