Diffuse non-thermal emission in the disks of the Magellanic Clouds

The Magellanic Clouds, two dwarf galaxy companions to the Milky Way, are among the Fermi Large Area Telescope (LAT) brightest gamma-ray sources. Aiming at a comprehensive modeling of the non-thermal electromagnetic and neutrino emission in both Clouds, we self-consistently model the radio and gamma-ray spectral energy distribution from their disks based on recently published Murchison Widefield Array and Fermi/LAT data. All relevant radiative processes involving relativistic and thermal electrons (synchrotron, Compton scattering, and bremsstrahlung) and relativistic protons (neutral-pion decay following interaction with thermal protons) are considered, using exact emission formulae. Our joint spectral analyses indicate that radio emission in the Clouds has both primary and secondary electron synchrotron and thermal bremsstrahlung origin, whereas gamma rays originate mostly from neutral-pion decay with some contributions from relativistic bremsstrahlung and Compton scattering off starlight. The proton spectra in both galaxies are modeled as power laws in energy with similar spectral indices, ~2.4, and energy densities, ~1 eV/cm3. The predicted 0.1-10 GeV neutrino flux is too low for detection by current and upcoming experiments. Our analyses confirm earlier suggestions of a largely hadronic origin of the gamma-ray emission in both Magellanic Clouds.


Introduction
The spectral energy distributions (SEDs) of cosmic rays (CRs) outside their galactic accelerators are important for determining basic properties of CR populations and for assessing the impact of the particle interactions in the magnetized plasma in galactic disks and halos. Knowledge of these distributions is generally limited as it is usually based only on spectral radio observations. When, in addition, non-thermal (NT) X/γ-ray observations are available, reasonably detailed spectral modeling of the CR electron (CRe) and proton (CRp) distributions in starforming environments can be very useful. Sampling the SED over these spectral regions yields important insight into the emitting CRe and possibly also CRp, whose interactions with the ambient plasma may dominate (via π 0 decay) high-energy emission ( ∼ > 100 MeV).
The nearby Magellanic Clouds (MCs), two irregular dwarf satellite galaxies of the Milky Way, may exemplify the level of spectral modeling currently feasible in a joint analysis of radio and γ-ray measurements. The Large MC (LMC) is located at a distance d = 50 kpc (Foreman et al. 2015); for the purpose of this study, its structure is modeled as a cylinder with radius r = 3.5 kpc and height h = 0.4 kpc (Ackermann et al. 2016). Similarly, for the Small MC (SMC) these quantities are d = 60 kpc, r = 1.6 kpc, and h = 4 kpc (Abdo et al. 2010b).
Both MCs have been detected with the Fermi Large Area Telescope (LAT) as extended sources (Abdo et al. 2010a,b).
Based on recent data (Lopez et al. 2018), there is statistically significant emission along the SMC "bar" and "wing" where there is active star formation. A spectro-spatial analysis of the LAT data suggests its integrated γ-ray emission to be mostly hadronic (≥90%). Similar analyses of the LMC with multi-source models (Foreman et al. 2015;Ackermann et al. 2016;Tang et al. 2017) attempted to extract the spectral content of these sources in a statistically viable way (including details on the quality of the fits). The sources were found to be point-like (pulsars, supernova remnants, plerions, unidentified background active galactic nuclei) and extended. The extended sources included the largescale disk denoted as 1 E0 ∼ G1, and smaller-scale components denoted as E1+E3 ∼ G2 (the region west of 30 Doradus); E2 ∼ G3 (the LMC northern region); and E4 ∼ G4 (the LMC western region) -these smaller-scale components possibly encompass different sources with different spectra. These spectral and/or spatial analyses revealed a galaxy-scale component dominating the emission, with a spectral shape suggesting a lepto-hadronic (Foreman et al. 2015) or hadronic (Ackermann et al. 2016;Tang et al. 2017) origin. 2 In this paper we focus on the emission in the MC disks in an attempt to determine the mean values of their magnetic fields and CRe and CRp energy densities. In spite of their detailed LAT data analyses, the above-mentioned studies were not based on self-consistent modeling of the broadband radio-γ SED of the two galaxies. SED modeling is important because a firm proof of the (mostly) pionic nature of the γ-ray emission should be based on a quantitative estimate of the CRe spectrum. Clearly, this spectrum is primarily determined from measurements of synchrotron radio emission, which was not accounted for in the earlier works. In an attempt to clarify, and possibly remove, some uncertainties inherent in previous modeling work, here we reassess key aspects of (average) conditions in the MC large-scale disks. Employing a one-zone model for the extended disk emission, we self-consistently carry out detailed calculations of the emission by CRe and CRp. The radio and γ-ray data used in our analyses are taken from the following papers: Lopez et al. (2018;SMC), Tang et al. (2017;LMC), Ackermann et al. (2016;LMC), SMC, LMC).
In Section 2 we briefly review the observations of extended NT emission from the MC disks. In section 3 we review the IRoptical radiation fields permeating the MCs. In Section 4 we describe calculations of the MC disk SEDs and perform fits to the data. Prospects of neutrino detection are discussed in Section 5, followed by a conclusion in Section 6.

Observations of extended emission
The MCs have been extensively observed over a wide range of radio-microwave, (soft) X-ray, and γ-ray bands. As mentioned above, point-source emission and extended small-and largescale emission have been detected from both galaxies. In this paper we focus on the extended large-scale disk emission of each galaxy because this emission traces the mean galactic properties of NT particles and magnetic fields. The spectral datasets used in our analysis are public (either tabulated or plotted) and are fully specified in Table 1 (SMC) and Table 2 (LMC). In this section we briefly review the observations most relevant to NT emission leaving out details that are discussed in the cited papers.
• Radio. In a multifrequency radio continuum study of the MCs, For et al. (2018) presented closely-sampled 76-227 MHz Murchison Widefield Array data, supplemented by previous radio measurements at lower (for the LMC) and higher (for both MCs) frequencies. Measured fluxes include emission from MC (and background) point sources, which were estimated to contribute 11% and 23% of the measured emission of the SMC and the LMC, respectively. Based on statistical analyses of a set of four fitting spectral models (power-law [1PL], curved PL, and double PL [2PL]) with either free and fixed high-ν index), they found that for the SMC the best-fitting model is a 1PL (in the 76-8550 MHz band) with α(85.5 MHz -8.55 GHz) = 0.82 ± 0.03 (consistent with Haynes et al. 1991). Instead, for the LMC the preferred model is a 2PL (19.7-8550 MHz) with a low-ν (19.7-408 MHz) index α 0 = 0.66 ± 0.08 (consistent with Klein et al. 1989) and a (fixed) high-ν index α 1 = 0.1 suggesting that synchrotron radiation and thermal free-free (ff) emission dominate at low and high frequencies, respectively.
• X-rays. The observed diffuse X-ray emission in the MCs has been determined to be of thermal origin (Wang et al. 1991: Einstein Observatory;Points et al. 2001: Röntgen Satellit (ROSAT) Position Sensitive Proportional Counters (PSPC); Nishiuchi et al. 1999Nishiuchi et al. , 2001: Advanced Satellite for Cosmology and Astrophysics (ASCA)); as such, it is not directly relevant to our SED analysis except for the estimated thermal plasma density, which is required to calculate the pionic emission (see below).
• Gamma rays. Both MCs were detected at >100 MeV γ rays as galaxy-scale extended sources whose emission is dominated by a diffuse component. This component has been interpreted as mainly originating from CRp interacting with the interstellar gas.
The SMC was first detected with 17 months of Fermi/LAT data (Abdo et al. 2010b) as a ∼3 o -size source, in which emission was not strongly correlated with prominent sites of star formation. More recently Lopez et al. (2018), based on 105 months of LAT Pass8 data, produced maps of the extended >2 GeV emission (no signal at >13 GeV) and found statistically significant emission along the galaxy-scale quietly star-forming bar and wing. Within a set of single-component spectral models (i.e., PL; broken PL, with frozen spectral shape and free normalization, representing pulsars; and exponentially cutoff PL, representing pionic emission) the exponentially cutoff PL provides the best fit to the total γ-ray spectrum, with only a marginal improvement of the fit if the broken PL is added. The Lopez et al. analysis suggests that although pulsars may contribute ≤10% at >100 MeV, the extended emission is mainly (≥90%) of pionic origin, with a γ-ray emissivity ∼ > 5 times smaller than the local (Galactic) emissivity.
The LMC was marginally (>4.5 σ) detected with the Compton Gamma Ray Observatory (CGRO) Energetic Gamma Ray Telescope (EGRET) (Sreekumar et al. 1992) and then confirmed (33 σ) with 11 months of Fermi/LAT data (Abdo et al. 2010a). More recent studies have focused on the spatial and spectral modeling of the γ-ray surface brightness distribution; these are briefly reviewed below. The first LAT-based spectro-spatial analysis of the γ-ray surface brightness distribution was performed by Foreman et al. (2015) using 5.5 years of data, in the photon energy range 0.2−20 GeV. They modeled the CR distribution and γ-ray production based on observed maps of the LMC interstellar medium, star formation, radiation fields, and radio emission. The γ-ray spectrum was described by means of analytical fitting functions for the π 0 -decay, bremsstrahlung, and inverse-Compton yields. These processes were estimated to account for, respectively, 50%, 44%, and 6% of the total 0.2−20 GeV emission. In particular, they inferred the CRp spectral index to be q p = 2.4 ± 0.2 and the equipartition (with CRp) magnetic field to be B eq = 2.8µG. A subsequent spectro-spatial analysis (Ackermann et al. 2016; A+16) of six years of LAT (P7REP) data in the 0.2-100 GeV band reported extended and point-like emissions. The dominant extended emission is in the form of a disk-scale component, denoted as E0 in their paper, and additional degree-scale emissions from several regions of enhanced star formation (including 30 Doradus). If pionic, the E0 emission component implies a population of CRp with ∼1/3 the local Galactic density, whereas the superposed small-scale emissions imply local enhancements of the CRp density by factors of at least 2-6. The spectrum of the E0 component (which does not include emission from the degree-scale emissions from several regions) shows some curvature (see their Fig.7-top): with no reference to the underlying emission mechanism, A+16 fitted the data with a tabulated function derived from the local (Galactic) gas emissivity spectrum (as a result of its interactions with CRp), an exponentially cutoff PL, a broken PL, and a log-parabola, and concluded that the bestfitting analytical model was the log-parabola. They noted, however, that the similarity between the log-parabola and tabulatedfunction models (see their Fig.9) suggests a pionic nature of the Article number, page 2 of 11 large-scale disk emission. More recently Tang et al. (2017;T+17) analyzed eight years of LAT Pass8 data and deduced a deeper, more spectrally extended 0.08-80 GeV spectrum of the large-scale disk component (identified as G1; this too did not include emission from other regions) largely overlaps with E0 of A+16. The G1 shape, now extending down to 60 MeV, was determined to be best described as a π 0decay hump from a CRp spectrum harder than that in our Galaxy, although they could not rule out CRe relativistic bremsstrahlung completely. In this analysis we use the T+17 G1 data as our reference set for the diffuse large-scale LMC disk emission, and the earlier A+16 E0 data as an auxiliary set for a consistency check.

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 cosmic microwave background (CMB) and the extragalactic background light (EBL). The CMB is a pure Planckian described by a temperature T CMB = 2.735 (1 + z) K and energy density u CMB = 0.25 (1 + z) 4 eV cm −3 . The EBL originates from direct and dustreprocessed starlight integrated over the star formation history the Universe. It shows two peaks, corresponding respectively to the cosmic infrared background (CIB, at ∼100 µm) that originates from dust-reprocessed starlight, and the cosmic optical background (COB, at ∼1 µm) that originates from direct starlight (e.g., Cooray 2016). The two peaks are described as diluted Planckians, characterized 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 aT 4 , where a is the Stefan-Boltzmann constant). The dilution factors of the cosmic fields are C CMB = 1, C CIB = 10 −5.629 , and C COB = 10 −13.726 . A recent updated EBL model, based on accurate galaxy counts in several spectral bands, is due to Franceschini & Rodighiero (2017); locally (z = 0) it can be numerically approximated as a combination of diluted Planckians, with A 1 = 10 −5.629 , T 1 = 29 K; A 2 = 10 −8.522 , T 2 = 96.7 K; A 3 = 10 −10.249 , T 3 = 223 K; A 4 = 10 −12.027 , T 4 = 580 K; A 5 = 10 −13.726 , T 5 = 2900 K; A 6 = 10 −15.027 , T 6 = 4350 K; A 7 = 10 −16.404 , T 7 = 5800 K; A 8 = 10 −17.027 , T 8 = 11600 K. Local radiation fields in the MCs arise from their intrinsic stellar populations and (given its close proximity) also from the Milky Way; we refer to this as the galactic foregound light (GFL). Similar to the EBL by shape and origin, the GFL is dominated by two thermal humps, IR and optical. The corresponding energy densities (in eV cm −3 ) are (i) LMC: u IR = 0.12 and u opt = 0.20, estimated from Foreman et al. 2015 (u IR+opt = 0.32) and the IR-opt SED (from NED 3 ); (ii) SMC: u IR = 0.026, u opt = 0.062, scaling down the LMC values taking (from the SED, cf. NED) the optical and IR peaks to be factors of 0.215 and 0.144, respectively, of the corresponding LMC values.

SED models
Our main objective in this study is to determine CR electron and proton spectra in the MC disks by spectral modeling of NT emission in all the relevant energy ranges accessible to the observations. The particle, gas density, magnetic, and radiation field distributions clearly vary significantly across the disk, obviating the need for a spectro-spatial treatment. A modeling approach based on a solution to the diffusion-advection equation has been applied in the study of the Galaxy and several nearby galaxies. However, even for just a diffusion-based treatment to be feasible and meaningful, the spatial profiles of key quantities, such as the particle acceleration sources, gas density, magnetic field, and generally also the diffusion coefficient have to be specified. Given the generally very limited observational basis for reliably determining these profiles, a parameter-intensive spectrospatial modeling approach is rarely warranted. An example is the very approximate diffusion-based approach adopted in modeling NT emission in the disks and halos of the star-forming galaxies NGC4631 and NGC4666 (Rephaeli & Sadeh 2019) for which reasonably detailed radio spectra and spatial profiles are avail-Article number, page 3 of 11 A&A proofs: manuscript no. pap_publisher

Fig. 2.
CRe spectra in the Magellanic Clouds. The curves are the secondary spectra (dashed), primary spectra (dotted), and the primaries' local PL representations (solid). In each panel the normalization of the primary CRe spectrum to its PL counterpart is chosen such that they nearly overlap. The primary CRe injection indices, q i , are 2.28 (SMC) and 2.23 (LMC). The two LMC CRe spectra panels refer, respectively, to the T+17 and A+16 Fermi/LAT datasets, as secondary spectra are derived from the pionic fits to the γ-ray data.
able; even so, the results of these analyses are not conclusive given the substantial uncertainty in the values of key parameters.
Radio emission in the MC disks is electron synchrotron in a disordered magnetic field whose mean value B is taken to be spatially uniform, and thermal ff from a warm ionized plasma. A significant CRp component could yield additional radio emission by secondary e ± produced by π ± decays, and γ-ray emission from π 0 decay. In addition, γ-ray emission is produced by CRe scattering off photons of local and background radiation fields. The calculations of the emissivities from all these processes are well known and standard.
We assume the CR spectral distributions to be timeindependent and locally isotropic with spectral PL form. The primary CRe spectral density (per cm 3 and per unit of the electron Lorentz factor γ) is, N e (γ) = N e,0 γ −q e for γ min < γ < γ max with γ min >> 1. As discussed later, this spectrum proves to be a good approximation to the actual primary steady-state spectrum in the relevant electron energy range. The secondary CRe spectrum can be analytically approximated as a smoothed 2PL, exponentially cut off at high energies (see below). The assumed

SMC
The emission spectrum from π 0 -decay has more constraining power than a generic PL; thus, our modeling procedure begins with fitting a pionic emission profile to the γ-ray data with free normalization, slope, and high-energy cutoff. We then use the deduced (essentially, fully determined) secondary CRe spectrum together with an assumed primary CRe spectrum to calculate the combined synchrotron emission using an observationally estimated value of B. The fit to the radio data also includes, at high frequencies, a thermal bremsstrahlung component computed with previously determined values of the gas density and temperature. Finally, the full Compton X-ray and γ-ray yield is determined.
The π 0 -decay γ-ray flux is computed using the emissivity in Eq. (15) of Persic & Rephaeli (2019a), where the total (thermal) gas proton density is n g = n HI + 2 n H 2 with the aver-Article number, page 4 of 11 M. Persic and Y. Rephaeli: Relativistic particles in the Magellanic Clouds Fermi/LAT spectrum (Fig.1). With these values, the estimated CRp energy density The closely related π ± -decay secondary CRe spectrum, N se (γ), has no free parameters once the CRp spectrum is determined 4 . We use N fit se (γ), with parameter values reported in Table  4, to compute the corresponding leptonic yields. However, the uncertainty in the total (HI and H 2 ) gas density clearly affects the energy loss rate b(γ) and the resulting spectral normalization and shape of N se (γ), hence also of N e (γ), and their yields. To compute the synchrotron emission we assume B = 3.5µG, based on estimates of the ordered (1.7 µG) and random (3 µG) fields in the SMC (Mao et al. 2008).
Once the secondary CRe synchrotron yield has been computed, its low-frequency residuals from the data are modeled using a CRe spectrum, N e (γ), with N e0 = 3.7 10 −8 cm −3 , q e = 2.23, γ max = 10 4 . This spectrum is the 1PL approximation, for 100 ∼ < γ ∼ < γ max , of the actual steady-state primary CRe spectrum. Neglecting diffusion and advection losses, this primary spec- is the CRe spectral injection rate. Over the mentioned electron energy range, we find that q i = 2.28 provides a decent match between the shapes of the curved spectrum and the 1PL spectrum; this value is typical for γ-ray emission from Galactic CR accelerators (e.g., supernova remnants, the Crab nebula). The relative normalization between the two spectra can be based on the measured flux density at some radio frequency or on imposing the same CRe energy density on the 4 Denoting the total cross-section for inelastic pp collisions σ pp (E p ) (see Eq.[79] in Kelner et al. 2006), the secondary CRe injection spectrum is Q se (γ) = (8/3) m e (c/k π ± ) n g N p0 [m p + (4m e /κ π ± ) γ] −qp σ pp (m p + E π ± /k π ± ) for γ thr < γ < γ max se , with γ thr se = m π ± /(4m e ) = 68.5 and γ max se =m max π ± /(4m e ) ≃ (E max p − 3/2m p )/(4m e ). The corresponding steady- is the radiative energy loss term. In a magnetized medium consisting of ionized, neutral, and molecular gas, it is b(γ) = 2 j=0 b j (γ) where b j (γ) are loss terms appropriate to each gas phase (see Eqs. [4]-[6] of Rephaeli & Persic 2015). N se (γ) is analytically approximated by N fit se (γ) = N se,0 γ −q 1 (1 + γ/γ b1 ) q 1 −q 2 exp[−(γ/γ b2 ) η ], where q 1 , q 2 and γ b1 , γ b2 are the low-and high-energy spectral indices and breaks; η gauges the steepness of the high-end cutoff. two spectra over the electron energy range of interest (see discussion in Rephaeli & Persic 2015). In our case the normalization of N pe (γ), and hence k i , can be found by using both N pe (γ) and N se (γ) to compute the total (primary plus secondary) synchrotron yield, which in turn is fitted to the radio synchrotron data. In doing this, the acceptable match (in the relevant energy range) between the curved and PL spectra allows us to use N e (γ) without substantial loss of accuracy. Given this procedure, the uncertainties in the gas density ultimately affect N e (γ) as well.
The spectra of the synchrotron-emitting CRe are shown in Fig. 2. The resulting primary electron component is quite dominant, with ∼70% of the total synchrotron and Compton yields (Fig. 3), but the normalization, N e0 , is not strongly constrained due to its dependence on the assumed magnetic field. If measurements of NT-X-ray and ∼1-30 MeV γ-ray fluxes, which correspond respectively to the CMB and CIB peaks, become available, most of the uncertainty could be removed (e.g., Persic & Rephaeli 2019a). In addition to this modeling uncertainty, the considerable coupling between the CRe spectral parameters implies an appreciable range of the deduced values of N e0 and q e . This range can be roughly bracketed by a flatter slope (q e = 2.2) and lower normalization (N e0 = 2.9 10 −8 cm −3 ) or, alternatively, a steeper slope (q e = 2.3) and higher normalization (N e0 = 6.7 10 −8 cm −3 ), both with the same γ max as the reference model (Table 4).
At high frequencies a relatively flat, ∝ ν −0.1 , component represents diffuse thermal ff emission (Spitzer 1978). This flux may be gauged to the Hα flux if both emissions come from the same emitting volume (HII regions) because in this case the relevant warm-plasma parameters (temperature, density, filling factor) are the same. So the measured (optical) Hα flux may be used to predict the (radio) ff emission. We model the thermal ff emission combining Eqs. (3) and (4a) of Condon (1992) and using T e = 1.1 10 4 • K (Toribio San Cipriano et al. 2017) and F(H α ) = 1.6 10 −8 erg cm −2 s −1 (Kennicutt et al. 1995). The primary synchrotron and thermal ff normalizations are spectrally quite far apart so they do not significantly affect each other. The radio model is shown alongside data in Fig. 3-right.
Although they are subdominant, secondary CRe contribute appreciably to the total CRe population of the SMC. Whereas the primary spectrum is approximately PL, the secondary spectrum is clearly curved (see Fig. 2). The resulting curvature of the total CRe spectrum is reflected in the total synchrotron spectrum (computed using Eq. (9) of Persic & Rephaeli 2019a for N e (γ) and its straightforward generalization for N se (γ)). Therefore, the total radio spectrum, which consists of synchrotron and thermal ff emission at low and high frequencies, is not a smooth 2PL, but shows some extra structure. A third-order polynomial (in log units) outperforms the For et al. (2018) "preferred" 1PL model (∆BIC> 0; Table 3). This polynomial (four free parameters; Fig. 3-left) is matched by the physical model described above, namely a low-frequency component representing synchrotron emission (three free parameters) plus a high-frequency ∝ ν −0.1 PL representing thermal ff emission (one free parameter; Fig. 3-right).
With the CRe spectra essentially determined, we can calculate the Compton and NT-bremsstrahlung yields from CRe scattering off, respectively, CMB-EBL-GFL photons and thermal plasma nuclei using, respectively, Eqs. (2.42) and (3.1) of Blumenthal & Gould (1970). The photon fields are described in section 3; for the plasma nuclei we assume a density n i = 0.033 cm −3 (Mao et al. 2008) and metallicity Z 2 = 1.2 (Sasaki et al. 2002). Although the shape of the γ-ray spectrum is distinctly pionic, NT-bremsstrahlung peaks near the pionic peak (ν ∼ 10 23 Article number, page 6 of 11 M. Persic and Y. Rephaeli: Relativistic particles in the Magellanic Clouds  (4) and Table 4 of For et al. (2018). + Third-order polynomial: log(S ν ) = c 0 + c 1 x + c 2 x 2 + c 3 x 3 , with x = log(ν/GHz). • χ 2 calculations use actual fluxes (Table 3 of  Hz) where it contributes ∼3% of the measured flux; Comptonized starlight (EBL+GFL) emission, which peaks at ν ∼ 10 22 Hz, contributes ∼ < 1% of the flux at the lowest LAT data point. The predicted Comptonized CMB spectral flux is also shown in Fig. 4. We found no need to introduce a further broken PL component in the γ-ray band representing pulsars, in agreement with Lopez et al. (2018). The resulting SED model, overlaid on data, is plotted in Fig. 4. An upper limit on the mean magnetic field, B max , can be estimated by assuming that the measured radio emission is fully produced by secondary CRe. By matching the measured and predicted emission levels we infer a limiting value B lim ≃ 7 µG. This upper limit, which is twice the observationally deduced value, is unrealistically high because the shapes of the measured and predicted synchrotron spectra do not match well, most notably at low frequencies, and because primary CRe obviously contribute to the measured radio flux. A lower limit, B min , can estimated using the fact that if B decreases the secondary synchrotron yield decreases as well. In order to keep the synchrotron yield unchanged, the primary yield must increase by increasing the spectral normalization (N e0 ) and possibly also varying the spectral shape (q e , γ max ) of the primary CRe spectrum. At low energies the primary spectrum is steeper than the secondary, which means that at some point the resulting synchrotron spectrum becomes progressively and systematically steeper than the <230 GHz radio data; this yields B min ≃ 2 µG. Based on our SED analysis, the value B = 3.5µG suggested by Mao et al. (2008) lies comfortably between the estimated bounds.

LMC
Our fitting approach differs from that adopted for the SMC. The radio spectrum is best fit by a 2PL model interpreted as a combination of (individually 1PL) synchrotron and thermal ff components , which indicates that synchrotron emission is generated by a 1PL CRe spectrum. This is unlike the case of the SMC for which the interplay of primary and secondary CRe with different spectra results in a more structured, non-PL form. With the primary CRe spectrum determined from fitting the radio data, the Compton γ-ray yield is calculated. To this a modeled pionic component is added and fit to the γ-ray data. Doing so enables the determination of the secondary CRe spectrum and the resulting yields, which turn out to be very minor in comparison with those of primary CRe; thus, no iteration of the fitting process is required.
The radio flux is dominated by synchrotron emission with α = 0.66 at the lowest measured frequencies, and by thermal ff emission with α = −0.1 at the highest frequencies . Adopting these spectral indices for the two components, we fit the full radio database with T e = 1.3 10 4 • K and F(Hα) = 2 10 −7 erg cm −2 s −1 (Toribio San Cipriano et al. 2017;Kennicutt et al. 1995) for the thermal ff component, and B = 4.3µG (Gaensler et al. 2005) 5 , N e0 ≃ 3.5 10 −7 cm −3 and γ max ≃ 2.4 10 4 for the synchrotron. The cutoff is needed for consistency with the measured flux at ν ∼ > 3 10 23 Hz. We should mention at this point that subsequent results (see below) suggest that the secondary CRe are only ∼10% of the total. So the (quasi-)1PL primaries dominate the CRe population (Fig. 2) and the emerging synchrotron emission. Therefore, the combined synchrotron and thermal ff radio spectrum matches the smooth 2PL profile of the For et al. (2018) fit. The radio model is shown in Fig. 5.
Next, we calculate the X-ray and γ-ray Compton and NT bremsstrahlung yields from (primary) CRe scattering off CMB-EBL-GFL photons and thermal plasma nuclei. To this end, we assume the photon fields described in section 3 and, respectively, a plasma characterized by an average density n i = 0.012 cm −3 (Sasaki et al. 2002;Cox et al. 2006) and metallicity Z 2 = 1.2 (Sasaki et al. 2002). Both emissions fall short of reproducing the γ-ray data.
The pionic emission is computed with n HI = 1 cm −3 (from M HI = 3.8 10 8 M ⊙ ; Staveley-Smith et al. 2003) and n H 2 = 6.6 10 −2 cm −3 (from M H 2 = 5 10 7 M ⊙ , measured from CO emission; Fukui et al. 2008). From our fitting, the matching CRp spectrum has q p = 2.38 and E max p = 80 GeV for the T+17 data ( Fig. 6-top) and q p = 2.6 and E max p = 25 GeV for the A+16 data (Fig. 6-bottom). The corresponding CRp energy densities are, respectively, 1 and 1.5 eV cm −3 , the difference being largely due to the somewhat higher flux level in the A+16 data. In spite of the observational uncertainties in the two datasets, we see that for either dataset the pionic yield largely dominates the γ-ray emission, so the pionic nature of the measured γ-ray spectrum appears well established. In addition to the dominant pionic component a subdominant 7% (5%) leptonic contribution peaks at 500 (400) MeV, similarly to the pionic contribution: Comptonized COB and NT bremsstrahlung emissions contribute, respectively, 6.5% (4%) and 0.5% (1%) to the peak of the observed emission in the T+16 (A+16) database.
The secondary CRe spectrum is derived and fitted analytically, as described in section 4.1; the fitting parameters are listed in Table 4 for both sets of LAT data. As mentioned, in both cases the radiative contribution of the secondary CRe is minor; this means that the corresponding primary PL spectra are essentially   (Table 2) are shown as dots: γ-ray data are from T+17 (top; reference set) and A+16 (bottom; auxiliary set). Model predictions are shown as curves. Emission components are plotted by the following line types: synchrotron, dot-long-dashed; thermal ff, dot-shortdashed; total radio, solid; Comptonized CMB, long-dashed; Comptonized starlight (EBL+FGL), short-dashed; NT bremsstrahlung, dotted; pionic, solid. For each type of NT leptonic yield secondary, primary, and total emissions are denoted as curves of progressively higher flux. In the X-ray and γ-ray range, on top of the emissions resulting from the various radiative processes a solid line indicates the cumulative emission. the same, 6 and (for each dataset) no fitting iteration is required. The primary and secondary CRe spectra are shown in Fig. 2.
Primary CRe largely dominate the CRe population in the LMC, even more so than in the SMC. However, their density is effectively unconstrained because it depends on assuming a magnetic field value (even though deduced from Faraday-rotation measurements, as in our case) rather than using a value derived by modeling the NT-X-ray emission as Comptonized CMB radiation (e.g., Persic & Rephaeli 2019a,b) or the hard X-ray-MeV emission (when it becomes available) as Comptonized IR radiation. Appreciable modeling uncertainty also stems from the inter-dependence among primary CRe parameters which results in a range of acceptable values of the spectral index and normalization, q e = 2.26 with N e0 = 2.22 10 −7 cm −3 and γ max = 2.0 10 4 , as well as q e = 2.40 with N e0 = 6.21 10 −7 cm −3 and γ max = 2.7 10 4 . 7

Neutrino emission
With an apparent dominant π 0 -decay origin of the γ rays produced in the interstellar medium of both MCs, it is clear that π ±decay neutrinos are also produced. Our calculations (following Kelner et al. 2006) of the predicted muon-and electron-neutrino spectra of the two galaxies indicate (Fig. 7) that the broadly peaked (∼0.1-10 GeV) neutrino flux is too low for detection by current and upcoming ν projects. This conclusion is based on the estimated observation time needed to detect the LMC with an experiment with detection sensitivity comparable to the Antarcticabased IceCube+DeepCore Observatory, the most sensitive current or planned ν-detector at neutrino energies 10 < E ν /GeV < 100 (e.g., Bartos et al. 2013). The effective area of the observations (Abbasi et al. 2012) is A eff (E ν ) = 40 E ν 100 GeV 2 cm 2 (for ν µ , two times smaller for ν e ) in this energy range (Bartos et al. 2013). Only in the narrow energy range 10-50 GeV do the IceCube+DeepCore sensitivity and the LMC predicted diffuse spectral ν-flux effectively overlap (cf. the T+17 LAT dataset). In this band the LMC flux can be approximated as dN ν dE ν ∼ 10 −10  dominated by atmospheric neutrinos produced by cosmic rays in the northern hemisphere; its energy spectrum is approximately flat in the relevant energy range, at a level ∼10 2 (dΩ LMC /1.6 sr) Gev −1 yr −1 (Bartos et al. 2013 and references therein). So the net 10-50 GeV background rate is ∼4×10 3 (∆Ω LMC /1.6 sr) yr −1 ∼ 60 yr −1 assuming the LMC angular radius to be 5 deg). Based on this crude estimate, observation of diffuse GeV neutrinos from the LMC disk would imply S /N < 0.2.

Conclusion
The SMC and the LMC, star-forming Galactic satellite galaxies, are among the brightest sources in the Fermi/LAT γ-ray sky. We self-consistently modeled the radio-γ SED of both galaxies using the latest available radio and LAT data, using exact emissivity formulae for the relevant emission processes. Both SEDs were modeled with the radio data interpreted as a combination of NT electron synchrotron emission and thermal electron bremsstrahlung, and the γ-ray data as a combination of π 0decay emission from CRp interacting with the ambient gas plus leptonic emission. The CRp spectra appear similar in the two galaxies with q p ∼ 2.4, and u p ∼ 1 eV cm −3 . Our spectral findings are qualitatively and quantitatively new for the SMC, and confirm and strengthen previous results for the LMC. In detail, we quantified the Lopez et al. (2018) suggestion of a mostly pionic origin of the LAT-measured γ-ray emission of the SMC, only ∼ < 3% of the peak flux being by NT bremsstrahlung; for the LMC, our analysis suggests that 93%, 6.5%, and 0.5% of the LMC peak γ-ray flux (at 0.5 GeV) are accounted for by pionic, Comptonized starlight, and NT-bremsstrahlung emission, in broad agreement with Foreman et al. (2015). As we have noted, for the SMC disk we used the Lopez et al. (2018) LAT dataset which is only mildly contaminated by an estimated ∼10% contribution from unresolved pulsars, whereas for the LMC disk we used LAT datasets based on maps that do not include emission from local gas and CR inhomogeneities associated with actively star-forming regions (A+16; T+17). Thus, our pionic emission modeling results are unlikely to be appreciably affected by emission from individual sources. As stated in Section 4, our treatment is based on determining the particle steady-state spectral distributions in the disk by fits to the radio and γ-ray data, using previously deduced mean disk values of the gas density, magnetic and radiation fields. An assessment of the combined uncertainties in the predicted particle radiative (and neutrino) yields is based first on the precision level of the observationally determined values of n g and B. Although the deduced CRe normalization is (essentially) independent of n g , its dependence on B is significant, N e0 ∝ B −(α+1) , which is roughly B −1.7 . The deduced CRp normalization is N p0 ∝ 1/n g . Given the substantial uncertainties in both n g and B, it is clear that these also impact the overall normalization of the particle spectra, which can be uncertain by a factor of 2-3. However, the constraining power of the fits to the measurements is mostly in the well-known spectral shapes of the predicted synchrotron and π 0 -decay processes. As explained in the previous section, the spectral fits are quite good, with an average level of statistical uncertainty of ∼20% on the primary CRe spectral parameters. Even though substantial, the overall level of uncertainty is not large enough to question our qualitative conclusion on the pionic nature of the γ-ray emission.
The SMC radio data have large error bars below 230 MHz where q e would be best determined, are unsampled (but for one point) between 230 and 1400 MHz where the effect of the spectral truncation would be observed, and are only sparsely sampled above 1400 MHz. This leaves both q e and γ max poorly determined. The LMC radio spectral data permit determination only of q e as no truncation is apparent; we estimated γ max by estimating its lower limit implied by the radio spectrum and its upper limit implied by the γ spectrum. A larger value would shift the Comptonized starlight (EBL+FGL) blue peak to higher frequencies and to a higher flux level, deteriorating the fit to the γ spectrum at frequencies log( ν) > 24. Our joint analysis of the broadband SED of the MCs using exact emissivity formulae, and the leptonic radio and γ emissions are coupled and are based on independent measures (independent of particle-field energy equipartition considerations) of the magnetic field, so the pionic emission is formally determined by modeling the residuals of the LAT data after the leptonic yield is accounted for. In general, uncertainties on the CRe spectrum affect the determination of the CRp spectrum because the latter is obtained by fitting the γ data once the leptonic yields have been properly accounted for. However, we feel our CRp results here are reasonably solid because of the clear dominance of pionic emission in the γ spectrum.
For the SMC the LAT data are sufficient to define the CRp spectral parameters, q p = 2.4 and E max p = 30 GeV. For the LMC, however, some inconsistency between the T+17 and A+16 LAT datasets implies the derived CRp spectrum to be, respectively, flatter (q p = 2.4) with no clearly discernible high-energy truncation (E max p ∼ < 80 GeV) and steeper (q p = 2.6) with a clear truncation (E max p = 25 GeV) and a ∼50% higher normalization. Both q p values are consistent within the errors with the Foreman et al. (2015) earlier estimate. The reason for the discrepancy between the two datasets is not obvious, as the T+17 and A+16 data extraction regions largely overlap with each other and with the LMC disk, and both have been cleaned for identified sources (point-like or extended). It should be noted that in the T+17 dataset the quality of the four highest-frequency points does not allow a definite estimate of E max p , whereas the four lowest-frequency points enable a clear view of the rising portion of the pionic hump. Ongoing Fermi/LAT measurements would likely improve the quality of the spectrum and clarify this issue. For both galaxies we obtained truncation energies systematically higher for CRp than for CRe, probably because energy Article number, page 9 of 11 A&A proofs: manuscript no. pap_publisher losses are much less efficient for the former than for the latter in the relatively low-density MC interstellar gas.
The secondary CRe spectra are, at best, as well determined as those of the parent CRp. Our model secondary spectra (with no free parameters) are relatively well defined in the MCs. However, the primary CRe spectra, determined by fitting radio data after secondary CRe are accounted for, are less well constrained. First, the measured radio fluxes include emission from individual background and MC-disk point sources that contribute ∼20% of the total emission and are spectrally similar to the overall emission , so there is some intrinsic uncertainty in the predicted extended emission's spectral shape and normalization (N e0 is biased high). Secondly, because of the q e -N e0 degeneracy the uncertainty on N e0 affects the determination of the spectral slope. Thus, primary CRe are rather poorly constrained in both MCs. In spite of this, our assumption of a (truncated) 1PL representing primary CRe is sound. The PL is a fair approximation, over the relevant CRe energy range, of the actual spectrum computed accounting for realistic energy losses in the MC disk environments. The approximation is particularly good for the LMC.
A possible underestimate of the gas content in the MC due to the presence of dark (unseen) neutral gas (i.e., HI gas) optically thick in the 21 cm line and/or H 2 gas with no associated CO emission, could increase the uncertainty on the total (atomic and molecular) gas density. The CO (J=1-0) transition, often used to detect the total H 2 mass, may be hard to detect in low-Z galaxies such as the MC (Rolleston et al. 2002;Requena-Torres et al. 2016), so the relatively bright [C II]λ158 µm line can be used instead: indeed, >70% of the total H 2 mass may not have been traced by CO (1-0) in dwarf galaxies, but is well traced by [C II]λ158 µm (Madden et al. 2020). Another important aspect to investigate is the mixing of dust into the interstellar medium (ISM) and the spatial variations of their properties (gas clumpiness affects emission levels, see above). Comparing the distributions of the IR dust emission and several gas tracers (Hα, 21 cm, CO emission), their different emission processes highlight the distribution of gas under different conditions. In the LMC this analysis reveals that (i) dust emission, sampled in the Multiband Imaging Photometer for Spitzer (MIPS) 70 and 160 µm bands, is well mixed with the large-scale HI 21 cm emission; (ii) Hα from star-forming H II regions is confined to the massive star formation region where also warmer (∼120 K) 24 µm emission is found; and (iii) the Spitzer InfraRed Array Camera (IRAC) 8 µm band that traces polycyclic aromatic hydrocarbons correlates well with the HI gas, but is absent from the massive star-forming H II regions. All in all, the dust emission revealed by the combined IRAC 8 µm and MIPS 24, 70, and 160 µm bands traces all three phases of the ISM gas (Meixner et al. 2006). Underestimating the gas content would affect the particle spectra in the following ways: (i) overestimating the CRp spectral normalization; (ii) underestimating the Coulomb and bremsstrahlung energy losses, b 0 (γ) and b 1 (γ) (Rephaeli & Persic 2015), which would cause the computed steady-state secondary CRe spectrum to be, respectively, higher and steeper; (iii) biasing the derived primary CRe spectrum lower and flatter.
Uncertainties in the 3D structure of the MCs (mainly for the SMC; Abdo et al. 2010b) directly affect only the CRe and CRp spectral normalizations, to match the spectral flux when the emitting volume changes; this holds for the relevant leptonic emissivities, which depend linearly on the CRe normalization. Thus, model SEDs are nearly unaffected by variations of the MC structure parameters. Another source of uncertainty in N e0 stems from assuming a magnetic field (for both galaxies), in our case field values derived by Faraday rotation studies of extragalactic polarized sources seen through the MCs. These field strengths are measured along the line of sight to the background sources; thus, they are line-averaged fields, whereas the field in the synchrotron emissivity expression is a volume average. The two values may or may not be the same; the potential mismatch introduces another systematic uncertainty in the determination of N e0 from fitting the radio spectrum. For example, in principle, if B is biased high, N e0 and all leptonic yields are biased low, so the resulting N p0 is biased high.
The main breakthrough needed in MC SED studies is the measurement of diffuse NT X-ray emission from their disks. As exemplified in our studies of radio-lobe SEDs (Persic & Rephaeli 2019a,b), the NT 1 keV flux, interpreted as Comptonized CMB emission, sets the normalization of the CRe spectrum. In the MC case, given that the excellent pionic fit to the γ-ray emission firmly defines the secondary spectrum, the 1 keV flux would effectively measure the primary normalization. Modeling the radio spectrum as synchrotron radiation would then provide the primary spectral shape and, importantly, the magnetic field. A spectral determination of B would bypass the need to assume particle-field equipartition, and would return a volume-averaged value better representing the mean field than the line averages yielded by Faraday rotation-based measurements. Deeper X-ray observations of both MCs, with better spatial and spectral resolution than currently available will be highly Article number, page 10 of 11 M. Persic and Y. Rephaeli: Relativistic particles in the Magellanic Clouds beneficial to searching the diffuse disk emission for a NT component.