PROPERTIES OF NEARBY STARBURST GALAXIES BASED ON THEIR DIFFUSE GAMMA-RAY EMISSION

The successful operation of the Fermi Gamma-ray Space Telescope provides our first opportunity to infer the potential origins of the diffuse high energy emission from normal galaxies beyond the Milky Way (e.g., Butt 2010). Recent detections of diffuse gamma-ray emission from Local Group, starburst, and Seyfert 2 galaxies demonstrate the system's sensitivity (Nolan et al. 2012; Abdo et al. 2010a; Lenain et al. 2011). With further integration time beyond these detection studies, we can construct gamma-ray spectra in order to constrain models of the cosmic-ray (CR) interactions that give rise to the emission in these sources.

For starburst galaxies in particular, the normalization of the gamma-ray spectrum helps to constrain, among other properties, the volume-averaged gas density affecting the CR distributions. With this and the radio normalization, the strength of the so-called calorimeter model of CR cooling may be determined for these sources. At issue is the discrepancy in massive starbursts between the minimum energy magnetic field strength, estimated from the radio flux, and the equipartition field strength, estimated from gas surface density and assuming hydrostatic balance (Thompson et al. 2006). Given their large surface densities and the linearity of the FIR–radio correlation over decades in flux, Thompson et al. (2006) argue for relatively high field strengths in starbursts. These values can be an order of magnitude above the minimum energy field estimates, and thus significantly alter our understanding of the energy budgets in starburst galaxies.

In this work, the gamma-ray, radio, and FIR spectra of the archetypal starburst galaxies NGC 253 and M82 are quantitatively compared to self-consistent models of the various energy loss and emission mechanisms for CRs in normal galaxies. A broad parameter space is explored, including density and magnetic field strength, to constrain the physical conditions of these starburst nuclei that best agree with their full spectral energy distributions (SEDs).

The data were obtained using the Fermi Large Area Telescope (LAT; Atwood et al. 2009). The LAT is a pair-conversion detector sensitive to photons with energies between roughly 20 MeV and 300 GeV. The publicly available Fermi Science Tools (Abdo et al. 2009) were used, and the threads reviewed on the Fermi Science Center Web site⁴ for the binned likelihood analysis were followed.

Almost 3.5 yr of Pass 7 LAT data were downloaded from 10° regions of interest around NGC 253 and M82. To avoid gamma rays from Earth's limb, a zenith cut of 100° was invoked in the data filtering step. The rock angle was constrained to be less than 52° and various other data flags and filters were set as suggested by the LAT team. With the large field of view of the LAT, the background must be modeled in order to ensure photons only associated with the target source are counted. The 2FGL catalog was used as an initial background, and then gttsmap was used to search for unmodeled sources that showed up with high significance. Two additional sources were identified within 5° of NGC 253, and three were found near M82. These new sources were incorporated into the background model with the likelihood routine gtlike. Thirty logarithmically spaced energy bins per were used for the counts cube and binned exposure map. The spatial binning for both was set to 01 per pixel. An additional 10° was added to the exposure map radius to account for the larger point-spread function (PSF) at low energies. The P7SOURCE_V6 instrument response functions were used. The final energy binning was settled upon to maximize the signal-to-noise and detection significance, based on the test statistic (TS), in each bin.

It should be noted that no source with a significant TS was detected above 100 MeV at the location of M81. However, the 061 separation between M82 and M81 is smaller than or comparable to the PSF of the LAT below 1 GeV, so they are currently difficult to resolve. For this work, it was assumed that none of the flux from the position of M82 originates from M81. The FIR and radio luminosities of M82 are more than an order of magnitude higher than those of M81, so the gamma-ray contribution from M81 is expected to be negligible.

For NGC 253, the integrated flux between 100 MeV and 100 GeV was (1.19 ± 0.15) × 10⁻⁸ photons cm⁻² s⁻¹, with a photon index of −2.24 ± 0.06, and TS of 113. For M82, a total flux of (1.44 ± 0.22) × 10⁻⁸ photons cm⁻² s⁻¹ was found, with a photon index of −2.25 ± 0.08, and TS of 180. These results are consistent with initial detections (Abdo et al. 2010b) and the second LAT catalog (Nolan et al. 2012). They also extend to meet the TeV observations (Acero et al. 2009; Acciari et al. 2009). Light curves show no detectable variability in eight-week bins. The differential spectra of both galaxies are shown in Figures 1 and 2, and are summarized in Table 1. A linear fit yields χ² values which indicate that the spectra are inconsistent with a simple power law. In particular, there may be an apparent turnover in the spectra below a GeV indicative of the expected CR proton interactions in these sources (Blom et al. 1999). The following section presents models of gamma-ray and radio emission from the expected CR populations in starburst nuclei.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The likely origin of gamma rays from starburst galaxies is the interaction of supernova-accelerated CRs with a large mass of dense, primarily molecular, gas. The high temperatures of dense molecular clouds in starbursts (e.g., Paglione et al. 2004; Sakamoto et al. 2011) indicate significant penetration and ionization by CRs (Suchkov et al. 1993). The CR electrons themselves emit synchrotron radiation, so radio spectra of starburst regions serve as a significant constraint as mentioned before. The imbedded, dusty star-forming regions should radiate strongly in the FIR, and if the CR electrons lose most of their energy within them, this would be a natural explanation for the FIR–radio correlation and calorimetry. However, the gamma-ray spectra are critical for helping to reduce degeneracies and ambiguities in the physical solutions.

To quantify the starburst contribution to the gamma-ray and radio emission from starburst galaxies, the steady-state CR distributions were modeled. The models are significantly updated versions of earlier studies (Paglione et al. 1996), similar in approach to more recent work (e.g., Torres 2004). The steady-state CR diffusion-loss equation (Longair 1981) is now integrated over energy, where the original work solved a time integration, requiring a slightly different Green's function. For electrons, as before, energy losses due to ionization, non-thermal bremsstrahlung, synchrotron radiation, inverse Compton (IC) scattering, and escape are accounted for (Paglione et al. 1996; Torres 2004, their Appendix). The full Klein–Nishina cross section is now included (Schlickeiser 2002, p. 88) for IC scattering since relativistic effects become important at the high energies to which Fermi and ground-based arrays are now sensitive. Bremsstrahlung losses are now included in the continuous CR energy loss rate rather than approximating it with an energy-independent CR loss timescale. Also, the photon energy density for IC is now constrained using the starburst templates of Siebenmorgen & Krügel (2007) fit to the IR spectra. Finally, positrons (see below) are now included.

The CR protons are affected by ionization and escape losses, plus pion production above a total energy of 1.22 GeV. The cross-section parameterizations of Kamae et al. (2006) are now employed for all pion species (π^{±, 0}). Secondary CR electrons are created via knock-on or Coulomb interactions, and the decay of negatively charged pions. Secondary positrons are created through the decay of positive pions. The distributions of these secondary leptons from charged pions are also determined using the parameterizations of Kamae et al. (2006).

Power-law primary CR electron and proton injection spectra of the form Q(E) = KE^−s, in particles cm⁻³ s⁻¹ MeV⁻¹, are assumed, where E is the particle total energy. The previous model used a power law in kinetic energy, similar to other studies (Torres 2004), which has a conveniently simple overall normalization. The total energy is chosen in this work not only to be consistent with other theoretical studies (Lacki et al. 2010; Bell 1978), but as a deliberate contrast from the previous work to test the effect of a different injection spectrum. This choice mostly affects the proton spectrum at low energies, and thus the pion spectra, which, via secondary production, could alter both the radio and gamma-ray emission, particularly near the pion bump.

This is strictly a one-zone model only simulating the expected nuclear emission from starburst galaxies. No diffusion to, or emission from, the extended disk is considered. There are a variety of reasons for this choice besides basic simplicity. Given that these sources are so near the flux and angular resolution detection thresholds already, the much lower emissivities from the disk are insignificant to Fermi. We originally modeled the disk of NGC 253 (Paglione et al. 1996) and indeed found its gamma-ray emissivity to be negligible compared to the central starburst, in the same way as its radio emission is much reduced (Ulvestad 2000). Recently, Torres et al. (2012) estimated more realistically the steep decline (by an order of magnitude) of the CR proton density outside the inner starburst in NGC 253. The diffuse gamma-ray and radio emission from the Milky Way Galactic center also outshine the disk (Crocker et al. 2011). Our nearest neighbor, the LMC, one of the very few extragalactic sources with resolved gamma-ray emission, similarly only shows emission from its most massive site of star formation (Butt 2010).

Using a power-law injection spectrum in total energy also requires a new normalization. As usual, the normalization K is determined by the supernova rate, Ψ, starburst volume, V, and the efficiency of power transfer between the supernovae and CRs, η, such that

$$\begin{equation} \int ^\infty _{m_pc^2} Q(E) E_{{\rm kin}}\ dE = \eta P \Psi /V \. \end{equation} \tag{ 1 }$$

The supernova energy, P = 10⁵¹ erg. A minimum CR kinetic energy, E_kin ≳ 2m_pv²_s ∼ a few MeV, is required for acceleration in a shock front given shock velocities of 10⁴ km s⁻¹ (Bell 1978), and is much smaller than the proton rest mass energy. The relative normalization between the proton and electron CR distributions at high energies then depends on the injection spectral index s (Bell 1978; Lacki et al. 2010; Persic et al. 2008; Schlickeiser 2002, p. 472), and as such is no longer an independent free parameter:

$$\begin{equation} N_p/N_e = (m_p/m_e)^{(s-1)/2}\ \. \end{equation} \tag{ 2 }$$

Typical values of s = 2.0–2.4 yield proton-to-electron ratios of 50–200, which are consistent with Milky Way observations and effectively constrain the choice of s. Note that this is a simplification which differs from the treatment of Lacki et al. (2010; see also discussion in Persic et al. 2008) in order to reduce the number of free parameters. This normalization mostly impacts the secondary contribution to the CR electron+positron flux, and introduces variations of at most factors of ∼2 from other treatments.

A convection loss timescale τ_c = 1 Myr is assumed, which is consistent with estimates of starburst winds and CR transport time (or length) scales (Heesen et al. 2009; Murphy et al. 2008). So the total energy-dependent escape timescale is the combined diffusion and convection terms, or

$$\begin{equation} \tau _{{\rm esc}}^{-1}(E) = \tau _0^{-1} \beta E_{_{\!\rm GeV}}^{0.5} + \tau _c^{-1}\ \. \end{equation} \tag{ 3 }$$

In sum, the model variables are s, the volume-averaged gas (H₂ + H) density n, magnetic field strength B, diffusion timescale τ₀, and photon field energy density U_ph (from the integrated starburst SED template).

3.1. Model Results: Steady-state CR Distributions

Figure 3 shows the results from a few models broadly indicating the effects of varying density and magnetic field strength on the steady-state CR electron+positron distributions. Higher densities increase the ionization losses at low energies, leading to an overall flatter, or more curved, spectrum. High magnetic fields increase synchrotron radiation losses, subsequently decreasing the CR electron distribution, particularly at high energies. As a result, high energy emission from IC scattering and non-thermal bremsstrahlung are reduced. Losses from IC scattering are less important; the overall results are quite insensitive to changes in U_ph or the FIR spectrum template. Increased losses due to a shorter convection timescale (τ_c) are also relatively small, primarily lowering and slightly steepening the CR proton spectrum around TeV energies. Longer diffusion timescales (higher τ₀; see Torres et al. 2010) result in slightly harder and higher CR spectra above 10 GeV or so. This effect is also comparatively small, reaching factors of just up to a few above a TeV, and negligibly changing the radio spectrum.

Zoom In Zoom Out Reset image size

Altering the injection slope s affects the steady-state spectral slopes accordingly, but can also substantially change the contribution of secondary particles because of its effect on the proton-to-electron ratio. For example, for flatter slopes, secondary particles contribute more to the lepton distribution (Equation (2)), causing a more pronounced bump around 0.1 GeV in the spectrum. Note that this secondary bump is very subtle compared to the primary spectrum, so adjusting it by factors of a few changes the total lepton spectrum very little. Thus, the results are not very sensitive to whether a power-law injection spectrum in total or kinetic energy is chosen, or to the specific treatment of N_p/N_e.

The CR proton spectrum depends primarily on the density. As with the electrons, higher densities slightly flatten the proton spectrum at low energies from ionization losses. However, high energy losses are also increased via pion production, so the overall effect is that the proton distribution decreases, and slightly flattens, with density. The effect of using an injection spectrum that is a power law in total, rather than kinetic, energy mostly appears in the resultant pion distributions (and thus the secondary electron+positron distributions). The pion spectrum is marginally lowered (by factors of ∼2) at the 0.1 GeV bump, but is unaffected at high energies. This small difference seems to have a negligible effect on the resultant gamma-ray and radio emission spectra as the secondary contribution is slight, and the high energy pion distributions are unaltered.

3.2. Model Results: Gamma-Ray and Radio Spectra

Radio and gamma-ray spectra were generated from the steady-state CR distributions. Gamma rays are created by IC scattering and non-thermal bremsstralung from primary and secondary CR electrons+positrons, plus the decay of neutral pions from CR protons interacting with ambient gas. Radio spectra are generated by synchrotron radiation from these CR leptons, thermal free–free absorption and emission, plus warm dust emission at millimeter wavelengths. We neglect TeV absorption and pair production for now as their expected impact is relatively small and difficult to discern with current sensitivities (Torres 2004).

The gamma-ray spectrum of a star-forming galaxy above a GeV is dominated by pion decay from CR protons interacting with the interstellar medium (Figure 4). At lower gamma-ray energies, bremsstrahlung and/or IC scattering from CR electrons is important. The IC contribution is very sensitive to the high energy tail of the CR lepton distribution, and contributes most when there is a sufficient density of very high energy CRs, that is, when the magnetic field—and associated energy losses—is low. In fact, sensitive MeV data would be very useful to constrain the TeV CR lepton distributions in these galaxies. Higher densities enhance not only pion decay emission, but also non-thermal bremsstrahlung, since the secondary CR electron and positron distributions increase with density given that they too arise from CR protons (via charged pion decay). Therefore, the gamma-ray spectrum normalization is most sensitive to density, with some dependence on magnetic field strength below a GeV.

Zoom In Zoom Out Reset image size

The radio spectrum normalization is very sensitive to the magnetic field. Higher gas densities tend to reduce the radio signal somewhat by flattening the CR electron distribution. The flattening, or curvature, in the radio spectrum (Williams & Bower 2010; Thompson et al. 2006) is a result of ionization reactions affecting mostly lower energy CR electrons. Higher densities also enhance the bump feature in the secondary distributions, adding to the curvature. To a much lesser extent, the far-infrared (FIR) emission from dust also affects the radio and gamma-ray spectra by providing a deep reservoir of seed photons for IC scattering. So a high FIR luminosity increases CR leptonic losses, consequently diminishing slightly the non-thermal bremsstrahlung and synchrotron emission, while enhancing the IC gamma-ray flux. The FIR–radio correlation hints at the interrelatedness of all of these processes in star-forming galaxies. Neither the radio nor gamma-ray emission is very sensitive to the diffusion timescale. Together the radio and gamma-ray normalizations effectively constrain n and B.

3.3. Model Fits to NGC 253 and M82

The Spitzer IR spectra of NGC 253 and M82 are shown in Figures 5 and 6 with some closely matching starburst SEDs from Siebenmorgen & Krügel (2007). The photon densities from these models are 171 and 274 eV cm⁻³ for NGC 253 and M82, respectively. The assumed supernova rates are 0.1 and 0.08 yr⁻¹ for M82 and NGC 253, respectively, and their assumed distances are 3.9 and 3.5 Mpc.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The radio and gamma-ray spectra were simultaneously fit using hundreds of model runs over a broad parameter space. Given the possibility that the magnetic fields in such strong starbursts may be quite high (Thompson et al. 2006), field strengths up to mG were included. Despite the small effect, we also investigated very short and long diffusion timescales, from 0.1 to 30 Myr since recent gamma-ray studies of supernova–cloud interactions indicate that long timescales reproduce the observations best (Torres et al. 2010). The best-fitting models were determined using the minimum χ². The radio data were selectively chosen to include only the central starburst disks (data from within the inner arcmin, or central 300–500 pc where the molecular and millimeter continuum emission are strongest), not the more extended winds/halos. In starburst galaxies, this extended emission tends to have a much steeper spectrum (Williams & Bower 2010; Heesen et al. 2009; Elmouttie et al. 1997) and is derived from a more evolved CR population rather than one directly interacting with the central starburst and contributing to the pion decay gamma rays. The volumes for these galaxies were assumed to be disks of 150 pc radius and 70 pc thickness. This volume may underestimate their full molecular extent, especially in M82, but this distinction is not so critical here; fluxes are generated from the modeled emissivities by multiplying by the volume (Torres 2004), thus canceling with the 1/V term in the normalization (Equation (1)).

The best-fit models are shown with the data in Figures 7–9. There is a large mismatch in both the noise and number of data points between the radio and gamma rays. Thus, the scatter in the gamma-ray models is understandable since the normalization is dictated mostly by the radio data. To stress the physical implications of the new Fermi data, fits were done artificially increasing the weighting of the gamma-ray data by the ratio of the numbers of radio and gamma-ray points. Very similar solutions result, though over a narrower range of parameter space.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Lacki et al. (2011) compared the first Fermi catalog detections and TeV data to a collection of published model predictions and concluded that they were for the most part broadly consistent with each other (many of these models were constrained merely by non-detections at GeV energies and so considerably overestimated the observations). The scatter in the models in their Figure 1 mirrors the 3σ envelope of fits in Figures 7–9.

Normalizing the radio and gamma-ray models to match the spectra yields power transfer efficiencies η ≈ 6% for both galaxies. All fits within 3σ of the minimum χ² result in efficiencies well below 20%, which is consistent with theoretical predictions (Markiewicz et al. 1990).

Solutions with low χ² values (Figures 10 and 11) do include low magnetic fields around the minimum energy value, B ∼ 150 μG, but with very low gas densities of ∼100 cm⁻³. Based on studies of dust and molecular gas in these galaxies, volume-averaged densities of at least several hundred H₂ +H cm⁻³ are expected (Paglione et al. 2004; Sakamoto et al. 2011; Weiß et al. 2001). Models with higher densities imply higher magnetic fields of several hundred μG, approaching the calorimeter model expectations. Valid solutions also result for short confinement timescales. Torres et al. (2010) found that a diffusion coefficient of D = 10²⁶ cm² s⁻¹ is strongly constrained by Fermi and TeV observations of supernova–cloud interactions. Given τ₀ ≈ R²/2D, and setting R to the starburst radius of 150 pc, yields τ₀ = 8 Myr. We restrict the following discussion just to models where n > 300 cm⁻³ and τ₀ = 10 Myr.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

These results generally favor the high magnetic field calorimeter model for starburst CR populations (Thompson et al. 2006). The implication is that the CR cooling rate, primarily via synchrotron radiation and ionization losses, is faster that the escape rate. Figure 12 shows the loss timescales for electrons using models consistent with the radio and gamma-ray data for NGC 253. Here, the CR lifetime is simply approximated as either the energy loss timescale E/b(E), where b(E) is the CR energy loss rate (Torres 2004; Paglione et al. 1996), or the steady-state lifetime, N(E)/Q(E). Even for lower densities and magnetic fields, the electron loss timescales are shorter than the escape timescale. Even protons appear mostly calorimetric, especially for low energies and high densities. So it would seem that the condition that τ_loss ≪ τ_esc is satisfied for electrons, particularly for higher values of n and B and given that diffusion timescales are likely longer than 1 Myr (recall that N depends weakly on τ₀). Therefore, the minimum energy estimate for magnetic field strength may indeed underestimate the true value in starburst galaxies.

Zoom In Zoom Out Reset image size

Another argument in favor of the high field (high density) solutions comes from a rough examination of the transport and interactions of CRs during their lifetimes within the starburst region. CRs encounter a column of primarily molecular gas which has been measured from CO studies to have beam-averaged values of up to 10²² H₂ cm⁻² (Paglione et al. 2004, 2001; Weiß et al. 2001), where typical single-dish beam sizes encompass the entire central starburst region (300–600 pc). One may approximate this column density as nΔl, where Δl is akin to the mean free path of the CR. This path length is then of the order of 100 pc for low average gas densities (30 cm⁻³) and around a few pc for high values (10³ cm⁻³). The former result implies that the typical CR leaves the starburst region within its lifetime, rather than losing the majority of its energy within it, contradicting the steady-state assumption and calorimetry. But it is hard to understand independent studies of very warm cloud core temperatures in starbursts if CRs interact in such a limited way with molecular cloud cores (explored below). The high-density (and high magnetic field) solutions imply that CRs lose the majority of their energy on spatial scales of the order of the sizes of dense molecular cloud clumps or cores. This result is much more consistent with model expectations and observations from photodissociation and/or CR dominated regions (Matsushita et al. 2010; Papadopoulos 2010), the calorimeter model, the high measured cloud temperatures (Panuzzo et al. 2010), and the size scales of Galactic molecular clumps at such densities (McQuinn et al. 2002). Confirming this expectation are the measured volume densities and peak column densities of individual, ∼30 pc scale cloud complexes in NGC 253 and M82 (10^22–24 H₂ cm⁻² and 10^3–4 cm⁻³; Sakamoto et al. 2011; Weiß et al. 2001). Comparable magnetic field strengths, and volume and column densities, have been inferred across the full ∼40 pc extent of the massive star-forming cloud Sgr B2 in the Galactic center (Crutcher et al. 1996).

Further justification for the high-density and magnetic field solutions may be derived from the expected ionization rates from penetrating CRs in starburst molecular clouds. Ionization rates, ζ = ζ_e + ζ_p, are estimated using the best-fit models for CR electron+positron and proton distributions, and the cross sections and procedures from Padovani et al. (2009). Only the H₂ cross sections are used in this analysis since they are the most significant. The steady-state CR fluxes are extrapolated down to the H₂ ionization energy of 15.6 eV for the integration. Note that this extrapolation overestimates the true ionization rate since other absorption processes such as H₂ excitation and dissociation may occur as gas column densities get large. Padovani et al. (2009) indicate that at low energies (<1 MeV and <100 MeV for electrons and protons, respectively), the CR spectrum becomes rather flat, with a plateau that declines with the encountered column density. Deep within clouds, the CR proton contribution dominates over CR electron ionization given the apparently more gradual decline of ζ_p with column density.

This analysis yields median values of ζ ≈ 6 × 10⁻¹² s⁻¹ for NGC 253 and M82 (Figure 13). This CR ionization rate is a few orders of magnitude above the estimate in M82 based on CR heating and the observed cloud temperatures, ζ ∼ 4 × 10⁻¹⁵ (Suchkov et al. 1993), and factors of 1–10 above the rate for ULIRGs (Papadopoulos 2010). Taking ζ from Suchkov et al. (1993) as the fiducial rate necessary to achieve the observed cloud kinetic temperatures in M82 and NGC 253, and reducing the mean ζ from this work by that magnitude requires a column density of ∼10²² H₂ cm⁻² (Padovani et al. 2009), consistent with the observed column density. Further, Figure 13 indicates trends between ζ and density and magnetic field strength. Solutions for lower n and B require more attenuation of the low energy CR flux to match the fiducial ζ. From Padovani et al. (2009), column densities of up to 10²⁴ H₂ cm⁻² or more are needed to reduce ζ by the required factors of 10⁴ (Figure 14). Such column densities exceed the gas mass of their entire starburst regions, which further supports the higher density solutions, and thus the high B calorimeter model. At very high densities, however, lower columns are needed to reach the fiducial ζ, and $\Delta l = N_{{\rm H}_2}/n$ can yield path lengths much less than a pc, which is unreasonably small. Therefore, the ionization rate analysis disfavors the extremely high- and low-density (and magnetic field) solutions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, theoretical and observed scaling relations between B, and $N_{{\rm H}_2}$ and n complete the argument against very low magnetic fields in these starbursts. While the radio and gamma-ray solutions indicate a degeneracy between n and B (Figures 10 and 11), which also roughly follows the expected scaling relation $B \propto \sqrt{n}$, the magnetic field should also increase proportionately with $N_{{\rm H}_2}$ (Crutcher 1999). This scaling is shown in Figure 14. Again, magnetic fields of several hundred μG are favored.

Self-consistent models of the steady-state CR distributions in starburst galaxies match both the radio and gamma-ray spectra of NGC 253 and M82. The models solve the CR diffusion-loss equation given a power-law spectrum of supernova-accelerated particles, accounting for energy and particle losses due to ionization, synchrotron radiation, IC scattering, bremsstrahlung, secondary particle generation, and escape. The radio normalization depends mostly on the magnetic field, with some sensitivity to density in mostly the curvature of the spectrum. The gamma-ray spectrum, now observable with sufficient sensitivity with Fermi and TeV detections, is mostly normalized by the density, with some dependence on the magnetic field below a GeV. Therefore, with both radio and gamma-ray data, the densities and magnetic fields in a starburst region may be independently and less ambiguously constrained than from radio data alone. In particular, the question of calorimetry (high magnetic fields) in massive starbursts may be addressed.

In NGC 253 and M82, these model results favor the high magnetic field solutions anticipated by Thompson et al. (2006). The higher densities also implied are consistent with molecular gas studies. Further, the ionization rates, approximated from the steady-state CR distributions with high magnetic fields and densities, also match expectations. These results confirm that CR penetration and heating is an important contributor to the warm temperatures observed in starburst galaxies. These starburst nuclei may be populated to a significant filling factor with giant molecular clouds resembling Sgr B2 in magnetic field (∼500 μG), column density (10²³ H₂ cm⁻²), and density (10³ cm⁻³).

The authors thank the Fermi-LAT team for their support, particularly at local workshops and through the Fermi Science Center Web site. We are especially grateful to T. Kamae for his guidance. We thank Patrick Lii for his early contributions to the starburst modeling code, and the CUNY/AMNH REU program that sponsored him (NSF AST-0851198). We thank A. Marscher and A. Carramiñana for helpful discussions. This work was supported in part by the NASA New York Space Grant Consortium based at Cornell University (No. NNX10AI94H). This work was supported in part by grant 63388-00 41 from the Professional Staff Congress of the City University of New York. The authors especially thank the referee for vital suggestions that greatly strengthened this work. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This work is based in part on observations made with the Spitzer Space Telescope, obtained from the NASA/IPAC Infrared Science Archive, both of which are operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with the National Aeronautics and Space Administration.

Facility: Fermi (LAT) - Fermi Gamma-Ray Space Telescope (formerly GLAST)