WINDS, CLUMPS, AND INTERACTING COSMIC RAYS IN M82

M82 is a nearby (D = 3.9 Mpc; Sakai & Madore 1999) and well-studied example of a central starburst occurring in a moderate mass disk galaxy (e.g., O'Connell & Mangano 1978). The galaxy is nearly edge on, so that much of the starburst zone is optically obscured by the dusty ISM within the disk. This orientation, however, allows the galactic wind to be clearly seen in wavelengths extending from the X-rays to the radio. M82 is a member of the M81 galaxy group and is connected to M81 and its surroundings by tidal debris that is readily visible in H i 21 cm line maps (Yun et al. 1994; Chynoweth et al. 2008). This structure is a result of a close passage about 0.2 Gyr ago between M82 and M81 that is the probable trigger of the M82 starburst (Yun et al. 1994). The presence of extraplanar H i surrounding M82 is likely to also be a result of this interaction (Yun et al. 1993).

The stellar body of M82 is a high surface brightness disk with a moderate radial scale length that is only mildly distorted in its outer regions (Ichikawa et al. 1995). The kinematics of the central region are dominated by the presence of a bar that extends over much of the ∼400 pc diameter starburst region (Wills et al. 1997; Greve et al. 2002) and is surrounded by a molecular gas ring or tightly wound arms. This region contains the majority of the molecular mass of ∼3 × 10⁸ M_☉ (Naylor et al. 2010), which is comparable to the mass of H i surrounding the galaxy (Yun et al. 1993). The SFR within the starburst region of ∼10 M_☉ yr⁻¹ (e.g., Gao & Solomon 2004; Förster-Schreiber et al. 2003a) exceeds that of the entire Milky Way. M82 thus should exhaust its internal gas supply in ∼0.1 Gyr, which is significantly less than the timescales of typical galaxies.

M82's ISM consists of compact, dense molecular clouds and dense warm, ionized gas, both embedded in a hot, low density medium in approximate overall pressure equilibrium (Westmoquette et al. 2009). Pressures in the starburst region are extremely high, on the order of 0.5–1.0 × 10⁷ K cm⁻³ (O'Connell & Mangano 1978; Smith et al. 2006; Westmoquette et al. 2007). Observations of CO lines suggest molecular gas densities of 10³–10⁴ cm⁻³ with kinetic temperatures of  40 K (Wild et al. 1992; Mao et al. 2000), whereas triatomic molecules such as HCO and HNC naturally suggest higher densities between 10^4.0 and 10^4.8 cm⁻³ with temperatures between 50 and 500 K (e.g., Mühle et al. 2007; Fuente et al. 2008; Naylor et al. 2010). CO measurements also constrain the mass of the molecular gas to 3 ± 1 × 10⁸ M_☉ (Naylor et al. 2010; Wild et al. 1992), much of which may be in the form of fragmented clouds. Observations in the near-infrared give warm, ionized gas densities on the order of ∼10–600 cm⁻³ with a gas mass of ∼8 × 10⁶ M_☉ (Förster-Schreiber et al. 2001; Westmoquette et al. 2009).

Due its intense star formation, M82 produces its bolometric luminosity of ∼4 × 10¹⁰L_☉ in only a small region. Assuming a cylindrical geometry with a radius of 200 pc and scale height of 100 pc, we find an internal radiation field energy density of U_rad ≈ 10³ eV cm⁻³ in the starburst zone. Thus both the thermal pressure and energy density are substantially enhanced in M82 relative to these quantities in the ISM near the Sun. If the magnetic field energy density also is increased by factors of ∼10³ over that in the solar neighborhood, then we would expect to find magnetic fields in M82 with average strengths of B ∼ 150 μG.

M82 is well known for its highly visible galactic wind, which has been observed from the X-rays to the radio (e.g., O'Connell & Mangano 1978; Seaquist & Odegard 1991; Ohyama et al. 2002; Stevens et al. 2003; Engelbracht et al. 2006; Strickland & Heckman 2007). At most wavelengths the outflow has an approximately bipolar structure that emerges nearly perpendicular to M82's stellar disk. Optical observations of emission lines from ionized gas entrained in the wind indicate outflow speeds of ∼500–600 km s⁻¹ (Shopbell & Bland-Hawthorn 1998, and references therein). Fits to the X-rays, however, indicate higher terminal wind speeds of up to ∼1400–2200 km s⁻¹ (Strickland & Heckman 2009). Traditionally the M82 wind has been interpreted as a hot outflow powered by shock heating from supernovae (Chevalier & Clegg 1985), but cosmic rays also could be a factor in driving the wind (Breitschwerdt et al. 1993). The variation in speeds and excitation levels within the M82 wind indicate that it has a complex structure, and so may also be subject to a variety of driving mechanisms.

An advantage of M82 for cosmic ray studies is that the supernova rate can be directly determined through radio studies of the numbers and expansion rates of supernova remnants. Although values as high as ν_SN = 0.2–0.3 SN yr⁻¹ have been used for the supernova rate in M82 (de Cea del Pozo et al. 2009a), we adopt the more conservative estimate of ν_SN = 0.07 yr⁻¹ from Fenech et al. (2008). This is likely to be a lower bound on the true rate, and is consistent with other estimates of the M82 SFR (Fenech et al. 2010).

Supernovae are the assumed accelerators of cosmic rays. Of the energy resulting from the explosion, only a small fraction is transferred to cosmic rays. In the Milky Way, the typical value for this efficiency factor, η, is 10% (Blandford & Eichler 1987). We discuss models for η = 0.04–0.2 in Section 5. The majority of the energy transferred to cosmic rays goes into protons and only a small fraction goes into electrons. The ratio of protons to electrons is here assumed to be N_p/N_e ∼ 50 (Blandford & Eichler 1987). The spectrum of cosmic rays accelerated by supernovae is generally assumed to be a power law with spectral index p ∼ 2.1–2.2.

Of the available radio data sets, a collection of single-dish observations by Klein et al. (1988) has the largest range of frequencies (from 10⁷ to 10¹⁴ Hz). Thermal dust emission becomes important starting at ∼10¹¹ Hz. The lowest frequency data will likely overestimate the flux of the starburst due to a large beam size. Because of this we only consider data from 10⁸ to 10¹¹ Hz in the Klein et al. data set. More recently, Williams & Bower (2010) observed the starburst region from 1 to 7 GHz with the Allen Telescope Array. These interferometer measurements have significantly smaller fluxes than the single-dish observations: the interferometer fluxes are ∼12% smaller (see Section 4).

Starburst galaxies were anticipated to be γ-ray sources because of their high star-formation and supernova rates. Recently, γ-rays above 700 GeV were detected from M82 with the ground-based Cherenkov telescope VERITAS. Measurements were made between ∼0.9 and ∼5 TeV with a flux upper limit (99% confidence level) at ∼6.6 TeV (Acciari et al. 2009). M82 has also been detected in γ-rays with the Fermi γ-Ray Space Telescope at energies from 200 MeV to 300 GeV (Abdo et al. 2010).

The cosmic ray transport equation in a homogeneous medium (Longair 2011), with diffusion omitted and advective losses added, is

$$\begin{equation} \frac{\partial N(E,t)}{\partial t} = - \frac{\partial }{\partial E} \left[ \frac{{\it dE}}{{\it dt}} N(E,t) \right] + Q(E,t)-\frac{N(E,t)}{\tau _{{\rm adv}}}, \end{equation} \tag{ 1 }$$

where N(E, t)dE is the number density of particles with energy between E and E + dE at time t, dE/dt < 0 is the rate at which a particle's energy changes due to radiation and collisions, Q(E, t)dE is the rate at which cosmic rays with energies between E and E + dE are injected per unit volume, and τ_adv, which we assume is independent of E, is the rate at which particles are advected out of the region.

We assume Q is independent of t and seek time independent solutions of Equation (1). Although the time independent version of Equation (1) is a linear ODE, with solutions that can be written down in closed form, the solutions are not very useful for our parameter study because the complicated form of dE/dt makes them difficult to evaluate. Therefore, we make the approximation

$$\begin{equation} N(E) \approx Q(E) \tau (E), \end{equation} \tag{ 2 }$$

where

$$\begin{equation} \tau (E)^{-1} \equiv \tau _{{\rm adv}}^{-1} + \tau _{{\rm loss}}^{-1} \end{equation} \tag{ 3 }$$

is the total energy loss rate and

$$\begin{equation} \tau _{{\rm loss}} \equiv - \frac{E}{{\it dE}/{\it dt}} \end{equation} \tag{ 4 }$$

is the energy loss rate due to radiative and collisional processes. Equation (2) becomes exact for τ_adv/τ_loss ≪ 1, and is accurate to better than 10%–20% in cases where Q is a power law in E with index near 2, as is thought to be the case for acceleration of cosmic rays by strong shocks, the leading theory of cosmic ray origin.

As supernovae are the assumed drivers of cosmic ray acceleration, the source function for cosmic rays must be related to the total energy input from supernovae:

$$\begin{equation} \int _{E_{{\rm min}}}^{E_{{\rm max}}} Q(E) E {\it dE} = \frac{\eta \nu _{{\rm SN}} E_{51}}{V}, \end{equation} \tag{ 5 }$$

where ν_SN is the supernova rate, V is the volume of the starburst region, η is the fraction of the supernova energy transferred to cosmic rays, and E₅₁ = 1 is 10⁵¹ erg, the typical energy from a supernova explosion. Then, assuming the source function is a power law of the form Q(E)∝E^−p, the particle spectrum becomes

$$\begin{equation} N(E) = \frac{(p-2)}{E_{{\rm min}}^{-p+2}} \frac{\eta \nu _{{\rm SN}} E_{51}}{V} E^{-p} \tau (E). \end{equation} \tag{ 6 }$$

Energy losses for protons include ionization, Coulomb interactions, and pion production. For electrons, losses include ionization, bremsstrahlung, inverse Compton scattering, and synchrotron radiation. Cosmic ray lifetimes for values representative of our best fit models are plotted in Figure 1, and energy loss rates are plotted in Figure 2. Energy loss rates are explained in the Appendix.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Only a small fraction of the supernova power goes into cosmic-ray electrons, as compared to cosmic-ray protons. We assume the source functions for protons and electrons can be related by (Torres 2004)

$$\begin{equation} Q_{e}(E) = \frac{N_{e}}{N_{p}} Q_{p}(E). \end{equation} \tag{ 7 }$$

We assume the ratio of protons to electrons to be N_p/N_e = 50.

In this section, we discuss pions and the particles resulting from their decay. Above ∼1 GeV, the dominant energy loss for cosmic-ray protons is pion production, and so in this regime, the source function for pions can be calculated from the proton spectrum. Assuming a delta function approximation, which is appropriate for smooth power-law distributions of cosmic rays, the source function for pions is given by (Kelner et al. 2006)

$$\begin{equation} q_{\pi }(E_{\pi }) = c n \int \delta (E_{\pi } - K_{\pi }T_{p}) \sigma _{pp}(E_{p}) N_{p}(E_{p}) {\it dE}_{p}. \end{equation} \tag{ 8 }$$

Integrating over proton energy, this becomes

$$\begin{equation} q_{\pi } (E_{\pi }) = \tilde{n} \frac{c n}{K_{\pi }} \sigma _{pp} \left(m_{p}c^{2} + \frac{E_{\pi }}{K_{\pi }} \right) N_{p} \left(m_{p}c^{2} + \frac{E_{\pi }}{K_{\pi }} \right), \end{equation} \tag{ 9 }$$

where K_π ≈ 0.17 is the mean fraction of the kinetic energy of the proton, T_p, transferred to the pion per collision.

While the cross sections for pion production from proton–proton interactions have been widely studied, there is no clear agreement as to which parameterization should be used. Of those available, we considered the model of Kamae et al. (2005) which includes the effects of the diffraction dissociation process and scaling violations, and the model of Kelner et al. (2006), whose energy weighted total cross section is based on the numerical simulations of proton–proton interactions by the sibyll code. Comparing the two cross sections, the models disagree at energies between 1 and 10 GeV. However, the γ-ray source functions calculated from the two models are nearly identical. The differences in cross section add less than a factor of ∼2 in uncertainty (Domingo-Santamaría & Torres 2005; Lacki et al. 2013). As the Kelner et al. (2006) and Kamae et al. (2005) cross sections give such similar results, we chose to use Kelner et al.'s model for simplicity.

Kelner et al. (2006) give the total inelastic cross section as

$$\begin{eqnarray} \sigma _{{\rm inel}} (E_{p}) = &(34.3 + 1.88L + 0.25L^{2}) \nonumber \\ &\times \left[ 1 - \left(\frac{E_{{\rm th}}}{E_{p}} \right)^{4} \right]^{2} {\rm mb}, \end{eqnarray} \tag{ 10 }$$

where L = ln (E_p/TeV). This gives us the cross section for neutral pion production. For charged pions, we can multiply the inelastic cross section by the multiplicity of the charged pions. To obtain the multiplicity, we take the ratio of the inclusive cross sections for the production of charged pions to the inclusive cross section for the production of neutral pions. We use the inclusive cross sections modeled in Dermer (1986) for this purpose.

Charged pions are relatively short lived and quickly decay into muons and then into electrons, positrons, and neutrinos, π⁺ → μ⁺ + ν_μ and $\mu ^{+} \rightarrow e^{+} + \nu _{e} + \overline{\nu }_{\mu }$. Because the muon moves with nearly the same speed as the pion, in the laboratory frame, their source functions are equivalent, q_μ(γ_μ) ≃ q_π(γ_π). Then, the secondary electron and positron source functions are given by

$$\begin{equation} q_{e}(\gamma _{e}) = \int _{1}^{\gamma _{e}^{\prime }({\rm max})} d\gamma _{e}^{\prime } \frac{P(\gamma _{e}^{\prime })}{2\sqrt{\gamma _{e}^{\prime 2}-1}} \int _{\gamma _{\mu }^{-}}^{\gamma _{\mu }^{+}} d\gamma _{\mu } \frac{q_{\mu }(\gamma _{\mu })}{\sqrt{\gamma _{\mu }^{2}-1}}, \end{equation} \tag{ 11 }$$

where $\gamma _{\mu }^{\pm } = \gamma _{e}\gamma _{e}^{\prime } \pm \sqrt{\gamma _{e}^{2}-1} \sqrt{\gamma _{e}^{\prime 2}-1}$, $\gamma _{e}^{\prime }({\rm max}) = 104$ and the electron/positron distribution in the muon's rest frame is (see Schlickeiser 2002)

$$\begin{equation*} P(\gamma _{e}^{\prime }) = \frac{2 \gamma _{e}^{\prime 2}}{\gamma _{e}^{\prime 3}({\rm max})} \left(3 - \frac{2 \gamma _{e}^{\prime }}{\gamma _{e}^{\prime }({\rm max})} \right). \end{equation*}$$

Primary and secondary electron and positron spectra can be seen in Figure 3.

Zoom In Zoom Out Reset image size

Neutral pions have extremely short lifetimes and quickly decay into γ-rays, π⁰ → γ + γ. The source function for γ-rays is given by (Schlickeiser 2002; Pohl 1994)

$$\begin{equation} q_{\gamma ,\pi }(E_{\gamma }) = 2 \int _{E_{{\rm min}}}^{\infty } \frac{q_{\pi ^{0}}(E_{\pi })}{\sqrt{E_{\pi }^{2} - m_{\pi }^{2}c^{4}}} {\it dE}_{\pi }, \end{equation} \tag{ 12 }$$

where E_min = E_γ + (m_πc²)²/(4E_γ).

In addition to production by the decay of neutral pions, γ-rays are also produced by inverse Compton scattering and by relativistic, non-thermal bremsstrahlung from cosmic-ray electrons. In the case of production by bremsstrahlung, the source function for the γ-rays is given by (Stecker 1971)

$$\begin{equation} q_{\gamma ,{\rm brem}}(E_{\gamma }) = c n \sigma _{{\rm brem}} E_{\gamma }^{-1} \int _{E_{\gamma }}^{\infty } N_{e}(E_{e}) {\it dE}_{e}, \end{equation} \tag{ 13 }$$

where σ_brem = 3.38 × 10⁻²⁶ cm² and N_e(E_e) represents the combined cosmic-ray electron/positron spectrum.

For inverse Compton scattering, the γ-ray source function is given by (Rybicki & Lightman 1979)

$$\begin{equation} q_{\gamma ,{\rm IC}}(E_{\gamma }) = \frac{3 c \sigma _{T}}{16 \pi } \int _{0}^{\infty } d\epsilon \frac{v(\epsilon)}{\epsilon } \int _{\gamma _{{\rm min}}}^{\infty } d\gamma \frac{n_{e}(\gamma)}{\gamma ^{2}} F(q, \Gamma) \end{equation} \tag{ 14 }$$

with (Schlickeiser 2002)

$$\begin{equation*} \gamma _{{\rm min}} = \frac{E_{\gamma }}{(2 m_{e} c^{2})} \left[ 1 + \left(1 + \frac{m_{e}^{2}c^{4}}{\epsilon E_{\gamma }} \right)^{1/2} \right], \end{equation*}$$

where E_γ is the energy of the resulting γ-ray, is the energy of the incident photon, γ is the energy of the electron. Here, the electron spectrum, n_e(γ), is in units of cm⁻³. The function F(q, Γ) is part of the Klein–Nishina cross section and is given by (Blumenthal & Gould 1970)

$$\begin{equation*} F(q, \Gamma) = 2q \ln (q) + (1 + q - 2q^{2}) + \frac{\Gamma ^{2}q^{2} (1 - q)}{2(1 + \Gamma q)}, \end{equation*}$$

where

$$\begin{equation*} \Gamma = \frac{4 \epsilon \gamma }{(m_{e}c^{2})}\quad {\rm and }\quad q = \frac{E_{\gamma }}{\Gamma (\gamma m_{e}c^{2} - E_{\gamma })}. \end{equation*}$$

For the blackbody spectrum, v (), we use an isotropic, diluted, modified blackbody spectrum (Persic et al. 2008, and references therein)

$$\begin{equation} v(\epsilon) = \frac{C_{{\rm dil}}}{\pi ^{2} \hbar ^{3} c^{3}} \frac{\epsilon ^{2}}{e^{\epsilon / k T_{d}} - 1} \left(\frac{\epsilon }{\epsilon _{0}} \right)^{\sigma = 1}, \end{equation} \tag{ 15 }$$

where C_dil is a spatial dilution factor (given by the normalization U_rad = ∫v()d) and ₀ corresponds to ν = 2 × 10¹² Hz.

We expect that the majority of the radio flux comes from non-thermal synchrotron emission and thermal free–free emission. To calculate the radio synchrotron spectrum, we utilize the emission coefficient given by Longair (2011),

$$\begin{equation} j_{\nu }^{{\rm synch}} = \left(- \frac{{\it dE}}{{\it dt}} \right) N(E) \frac{{\it dE}}{d\nu }, \end{equation} \tag{ 16 }$$

where N(E) is the combined electron and positron spectrum and

$$\begin{equation} {-} \frac{{\it dE}}{{\it dt}} = \frac{4}{3} \sigma _{T} c \left(\frac{E}{m_{e}c^{2}} \right)^{2} \frac{B^{2}}{2 \mu _{0}} \end{equation} \tag{ 17 }$$

gives the synchrotron losses as a function of energy. If we make the simplifying assumption that an electron of energy E radiates away its energy at the critical frequency, ν_c, such that ν ≈ ν_c ≈ γ²ν_g ≈ γ² × eB/(2πm_e), then

$$\begin{equation} j_{\nu }^{{\rm synch}} = \frac{4}{3} \sigma _{T} m_{e} c^{3} \left(\frac{\nu }{\nu _{g}^{3}} \right)^{1/2} \frac{B^{2}}{2 \mu _{0}} N(E). \end{equation} \tag{ 18 }$$

Because of the warm, ionized gas in the ISM, we expect some portion of the radio spectrum to come from thermal bremsstrahlung or free–free emission and absorption. The free–free emission coefficient is given by (Rybicki & Lightman 1979)

$$\begin{equation} j_{\nu }^{{\rm ff}} = 6.8 \times 10^{-38} T^{-1/2} Z^{2} n_{e} n_{i} e^{-h\nu / k_{B}T} \overline{g}_{{\rm ff}} \end{equation} \tag{ 19 }$$

with units of erg s⁻¹ cm⁻³ Hz⁻¹. Additionally, we must consider the inverse process, free–free absorption. The absorption coefficient is given by (Rybicki & Lightman 1979)

$$\begin{equation} \kappa _{\nu }^{{\rm ff}} = 0.018 T^{-3/2} Z^{2} n_{e} n_{i} \nu ^{-2} \overline{g}_{{\rm ff}}\; {\rm cm}^{-1}. \end{equation} \tag{ 20 }$$

In the "small-angle, classical region," the mean Gaunt factor is (Novikov & Thorne 1973)

$$\begin{equation} \overline{g}_{{\rm ff}} = \frac{\sqrt{3}}{\pi } \ln \left[ \frac{1}{4 \zeta ^{5/2}} \left(\frac{k_{B}T}{h \nu } \right) \left(\frac{k_{B}T}{Z^{2} Ry} \right)^{1/2} \right], \end{equation} \tag{ 21 }$$

where Ry = 13.6 eV is the Rydberg constant and ζ = e^γ ≈ 1.781 with γ being Euler's constant.

For both emission and absorption, the radiative transfer equation is (Rybicki & Lightman 1979)

$$\begin{equation} \frac{dI_{\nu }}{d\tau _{\nu }} = -I_{\nu } + S_{\nu }, \end{equation} \tag{ 22 }$$

where S_ν ≡ j_ν/κ_ν is the source function. We discuss solutions to Equation (22) in Section 4.3.1.

We assume that the cosmic rays sample the mean density in the disk. This assumption is not trivial, since we expect that most of the cosmic rays are accelerated in the hot, low density medium, which is flowing out of the disk, while most of the mass is in the form of dense, molecular gas which has a small filling factor. Cosmic rays will sample the mean density if (1) most of the magnetic field lines they are accelerated on connect to a "fair sample" of the interstellar gas, and (2) cosmic rays can propagate a sufficient distance along the field lines, or diffuse sufficiently across the field lines, to encounter dense gas before they—and the field lines—are borne out of the disk by the wind (see Figure 4). The timescale for the latter is τ_adv, the escape time introduced in Equation (5).

Zoom In Zoom Out Reset image size

By a "fair sample," we mean that the ratio of interstellar gas mass to magnetic flux is uniform within the starburst region. Since the covering factor of the molecular clouds is below unity, this suggests that the field lines have a substantial random component, as is the case in the Milky Way. In the absence of any information to the contrary, we assume that the field lines are indeed uniformly loaded, so that the average gas density on any magnetic flux tube within the disk is the mean interstellar density, f_moln_mol + f_ionn_ion + (1 − f_mol − f_ion)n_hot ≈ f_moln_mol. If there are N_mc molecular clouds in the disk and the mean dimension of the disk is L ∼ V^1/3 where V = 2πR²H, then the mean distance l_c to a cloud is of order $L N_{{\rm mc}}^{-1/3}$. If the supernova remnants where the cosmic rays are accelerated and the molecular clouds are spatially uncorrelated, then l_c is the mean minimum distance that cosmic rays must travel to encounter dense gas.

In order to estimate the time it takes cosmic ray particles to travel a distance l_c, we assume that the cosmic rays are self-confined (see Kulsrud 2005 for a pedagogical treatment): any cosmic ray density gradient produced by a point source of cosmic rays, or by the global distribution of cosmic ray sources within the starburst region, is associated with an anisotropy of the cosmic rays in velocity space. If the mean cosmic ray velocity, or streaming speed, exceeds the Alfvén speed v_A, then the cosmic rays destabilize short wavelength Alfvén waves. The waves grow large enough to scatter the cosmic rays, reducing their drift speed to v_A. In this picture, the propagation time is l_c/v_A. The cosmic rays will encounter dense clouds if

$$\begin{equation} \frac{v_A \tau _{{\rm adv}}}{l_c} \sim \frac{v_A}{v_{{\rm adv}}} \left(\frac{H}{R} \right)^{2/3} \frac{N_{{\rm mc}}^{1/3}}{(2 \pi)^{1/3}} > 1, \end{equation} \tag{ 23 }$$

where in the last step we have identified τ_adv with an advection time H/v_adv.

In Section 4.4 we will show that both the radio and γ-ray spectra are very well fit by a simple family of models in which the magnetic field strength B and advection velocity v_adv are related by the empirical formula

$$\begin{equation} v_{{\rm adv}} = 4.0 B_{\mu } - 600, \end{equation} \tag{ 24 }$$

where B_μ is the magnetic field in μG and v_adv is in km s⁻¹. Equation (24) holds for 150 < B_μ < 350. Substituting this relation into Equation (23), using H/R = 0.5 (see Table 1) and solving for N_mc, we find that cosmic rays will encounter dense clouds before being advected out of the Galaxy if

$$\begin{equation} N_{{\rm mc}} > 8 \pi \left(2.00 - \frac{300}{B_{\mu }} \right)^3 n_{{\rm hot}}^{3/2}. \end{equation} \tag{ 25 }$$

According to Equation (25), for the value of N_mc needed to satisfy the criterion varies from 0 to 7 as B_μ ranges from 150 to 350 with n_hot = 0.33 cm⁻³ (see Table 1). With the molecular gas mass of the M82 starburst region estimated to be about 3 × 10⁸ M_☉, and giant molecular cloud masses in the Milky Way estimated at 10⁴–10⁶ M_☉, this criterion would seem to be easily satisfied. Moreover, the drift speed actually exceeds v_A if the waves are strongly damped, making Equation (23) even easier to satisfy. A more careful assessment based on a Monte Carlo simulation supports this conclusion (E. Boettcher et al. 2013, in preparation). This effect may be important for the higher energy cosmic rays. Once the cosmic rays reach the weakly ionized clouds, the waves which scatter them are strongly damped, and the cosmic rays interact primarily through physical collisions (Everett & Zweibel 2011).

Though we incorporate the multiphase nature of the ISM into our model for the radio spectrum (see Section 4.3.1), this is not necessary in the case of the γ-rays because of our assumption that they sample the mean density (Section 4.1) which is almost entirely due to molecular clouds. The ionized gas mass is well-constrained by observations to be M_ion = 8 × 10⁶ M_☉. In Section 4.4 we describe models with molecular masses of M_mol = 2–4 × 10⁸ M_☉ corresponding to mean densities of 〈n〉 = 280–550 cm⁻³.

Pion decay is the main mechanism for the production of γ-rays. When advective losses dominate over energy losses, the γ-ray flux is proportional to the gas density. However, in the calorimeter limit, where cosmic-ray proton energy losses dominate over advective losses, the γ-ray flux is essentially independent of density and so the number of γ-rays depends only on the number of protons. Other parameters that affect the γ-ray spectrum include the supernova rate (ν_SN = 0.07–0.1 SN yr⁻¹) and the distance (d = 3.5–3.9 Mpc), both of which are observationally constrained.

In addition to the supernova rate and distance, wind advection speed (v_adv) affects the γ-ray spectrum. For our cosmic ray spectra, we use a combined lifetime (see Equation (3)) which includes energy losses and advection. The advection time is given by τ_adv = H/v_adv where H is the height of the starburst region and v_adv is the speed of particles in the wind in the starburst region. Of the parameters that affect the γ-ray spectrum, the least constrained is the wind speed. For the γ-ray spectrum, the main constraint on the wind speed is the Fermi observations. The proton lifetime peaks between 1 to 10 GeV, and it is this portion of the proton spectrum that contributes to the energy range for the γ-ray spectrum observed by Fermi (0.3–30 GeV). So, if the wind speed becomes too large, the maximum proton lifetime decreases and the γ-ray spectrum will no longer agree with the Fermi observations. If this is the case, then the galaxy is no longer calorimetric.

Assuming values typical of the Milky Way for parameters related to cosmic ray acceleration by supernovae, we find that the γ-ray flux from neutral pion decay fits the TeV energy γ-ray data and all but the two lowest-energy Fermi data points as seen in Figure 5 and discussed in Section 4.4. As the γ-ray flux from neutral pion decay does not account for all of the observed flux below ∼1 GeV, there must be another emission mechanism responsible for the total γ-ray flux at these energies. Two possible mechanisms are inverse Compton emission and bremsstrahlung. For γ-ray emission from the inverse Compton effect, we found that the γ-ray flux in the starburst region was on the order of 2 × 10⁻¹⁰ GeV cm⁻² s⁻¹ at ∼1 GeV assuming an approximate radiation field energy density of 1000 eV cm⁻³ (see Table 1). While our model currently includes inverse Compton emission from only primary electrons, we found that it to be a significant contribution from the total γ-ray flux at both GeV and TeV energies. Inclusion of emission from secondary electrons and positrons will increase the total γ-ray flux by a factor of three at minimum. We will provide a more accurate model for inverse Compton emission in future work.

Zoom In Zoom Out Reset image size

We also calculated the γ-ray flux from non-thermal, relativistic bremsstrahlung. As bremsstrahlung depends on the electron spectrum, the flux is dependent upon the magnetic field strength and wind speed. We find that bremsstrahlung makes a significant contribution to the γ-ray flux below ∼10 GeV (see Figure 5). Thus, to accurately model the observed γ-ray flux at lower energies, we must include both neutral pion decay and bremsstrahlung in our work.

In the case of the radio spectrum, we have several more parameters that affect our fit than for the γ-ray spectrum. These parameters are the ionized gas density (n_ion), the magnetic field strength (B), and the wind speed (v_adv). Free–free emission and absorption both scale as the square of the gas density. Thus, the lower the gas density, the lower the frequency at which the radio spectrum turns over and the lower the contribution of free–free emission to the radio spectrum at higher frequencies. Additionally, synchrotron emission is scaled by the square of the magnetic field strength (B), and the wind speed (v_adv) affects the cosmic ray "dwell time" in the starburst. To quantitatively distinguish between the effects of our parameters, we use χ² tests to calculate the goodness-of-fit.

4.3.1. Multiphase Solutions for Radiative Transfer

The emergent radio spectrum is due to free–free and synchrotron emission, modified by free–free absorption as discussed in Section 3.3. Fitting a model to the observed radio spectral energy distribution therefore depends on knowing the degree of free–free absorption. This is most pronounced at lower frequencies and these observations are essential for obtaining a physically realistic fit to the radio data for M82. We tested several models and found that the free–free covering factor was always near unity. A simpler model consisting of a cylindrical wall of warm, ionized gas surrounding a lower density hot medium provides a reasonable representation of the free–free radio absorption in M82.

The multiphase nature of the ISM is incorporated by considering the different properties of the warm and hot ionized gas. The warm, ionized gas is assumed to be magnetically threaded and produces both synchrotron and free–free emission ($j_{{\rm WIM}} = j_{\nu }^{{\rm synch}} + j_{\nu }^{{\rm ff,ion}}$) as well as free–free absorption ($\kappa _{{\rm WIM}} = \kappa _{\nu }^{{\rm ff,ion}}$). Synchrotron self-absorption is found to be negligible in M82, so we need only consider free–free absorption. The hot, low density gas also produces synchrotron, and due to its lower density, negligible amounts of free–free emission and absorption ($j_{{\rm hot}} = j_{\nu }^{{\rm synch}} + j_{\nu }^{{\rm ff,hot}} \approx j_{\nu }^{{\rm synch}}$).

Adopting our assumed cylindrical geometry with a shell of warm, ionized gas, we take the shell width to be l/2 where l is the radius of the cylinder. Because the galaxy is nearly edge-on, we see emission from both the front and the back of the cylinder. Emission from the front of the cylinder is

$$\begin{equation} I_{\nu , {\rm front}} = \frac{j_{{\rm WIM}}}{\kappa _{{\rm WIM}}} (1 - e^{- \kappa _{{\rm WIM}} l/2} ). \end{equation} \tag{ 26 }$$

Emission from the rear of the cylinder can be absorbed within both the front and back warm ionized gas shells. Thus,

$$\begin{equation} I_{\nu , {\rm back}} = \frac{j_{{\rm WIM}}}{\kappa _{{\rm WIM}}} (1 - e^{- \kappa _{{\rm WIM}} l/2} ) e^{- \kappa _{{\rm WIM}} l/2}. \end{equation} \tag{ 27 }$$

Then, the radiative intensity emergent from the hot, diffuse gas is given by

$$\begin{equation} I_{\nu , {\rm hot}} = 2 j_{{\rm hot}} r_{i} e^{- \kappa _{{\rm WIM}} l/2}, \end{equation} \tag{ 28 }$$

where r_i is the radius of the central volume of hot, diffuse gas (r_starburst = r_i + l/2). Then, the total emergent intensity from the starburst region is

$$\begin{equation} I_{\nu } = \frac{j_{{\rm WIM}}}{\kappa _{{\rm WIM}}} (1 - e^{- \kappa _{{\rm WIM}} l} ) + 2 j_{{\rm hot}} r_{i} e^{- \kappa _{{\rm WIM}} l/2}. \end{equation} \tag{ 29 }$$

This model ignores the presence of a radio halo. Observations of the radio halo shows it has a steeper spectral index than that of the main body and thus can become dominant at low frequencies (Seaquist & Odegard 1991; Adebahr et al. 2012). We therefore constructed a model that assumes all of the 74 MHz flux from M82 arises in an unabsorbed synchrotron halo. This assumption is consistent with the parameters of our cylindrical model, where any 74 MHz emission from the main starburst zone is heavily absorbed. We assume a standard −0.7 spectral index for the halo. The radio halo's spectral energy distribution then is simply

$$\begin{equation} I_{\nu , {\rm halo}} = I_{0, {\rm halo}} \nu ^{-0.7}, \end{equation} \tag{ 30 }$$

where I_{0, halo} is determined by our lowest frequency radio data point. The final expression for intensity is

$$\begin{equation} I_{\nu } = \frac{j_{{\rm WIM}}}{\kappa _{{\rm WIM}}} (1 - e^{- \kappa _{{\rm WIM}} l} ) + 2 j_{{\rm hot}} r_{i} e^{- \kappa _{{\rm WIM}} l/2} + I_{\nu , {\rm halo}}. \end{equation} \tag{ 31 }$$

Figure 6 shows the radio spectrum predicted by the simple cylindrical model, which shows the contributions from the model components. Emission from the hot, diffuse gas, Equation (28), is represented by the dotted blue line; emission from the warm, ionized gas, Equations (26) and (27), is represented by the dot-dashed green line; and emission from the halo, Equation (30), is represented by the dashed red line. Note that the free–free high frequency emission is underestimated since we cannot use a single value for the mean electron density (see Section 4.4). Our best-fit model predicts a flux from the starburst region at 333 MHz in agreement with direct measurements from Adebahr et al. (2012).

Zoom In Zoom Out Reset image size

4.3.2. Radio Observations

4.4.1. χ² Tests

One of the primary goals of developing a one-zone model is to be able to find a set of parameters which produce an acceptable fit to both the radio and γ-ray spectra. To find the best-fit combined solutions, we use χ² tests on both spectra. We begin by using observations and previous studies to determine the range of values to test for each parameter. Previous studies have calculated the equipartition magnetic field strength to be ∼150 μG (de Cea del Pozo et al. 2009b). Thus, we test values from 50 to 350 μG. For the wind speed, model fits to X-ray observations of the wind give wind speeds on the order of 1500 km s⁻¹ (Strickland & Heckman 2009), while optical observations give wind speeds on the order of 500 km s⁻¹ (Shopbell & Bland-Hawthorn 1998). The two measurements do not inherently contradict each other as they measure different materials in the wind. Thus, we test values from 0 to 2000 km s⁻¹. We have fixed the total mass of the molecular and ionized gases, and in doing so, we have also fixed the proton density. Observations constrain the ionized gas densities to be 10 to 600 cm⁻³ with kinetic temperatures of ∼8000 K (Förster-Schreiber et al. 2001). We test values from 50 to 500 cm⁻³ (see Table 2). Additionally, we test two different cosmic ray spectral indices (p = 2.1, 2.2) and three different molecular gas masses (M_mol = 2 × 10⁸, 3 × 10⁸, 4 × 10⁸ M_☉). We took advantage of the computational simplicity of our model to run 49,000+ cases.

Table 2. Fitted Model Parameters

Physical Parameters	Tested Range	Reference
Magnetic field strength (B)	50–350 μG	1, 2
Advection (wind) speed (vadv)	0–2000 km s−1	3, 4
Ionized gas density (nion)	50–500 cm−3	5

References. (1) de Cea del Pozo et al. 2009b; (2) Thompson et al. 2006; (3) Shopbell & Bland-Hawthorn 1998; (4) Strickland & Heckman 2009; (5) Förster-Schreiber et al. 2001.

Download table as:  ASCIITypeset image

In determining the best fit to the radio spectrum, we performed χ² tests using the available radio observations, discussed above. Fits to the 11 radio data points were used to constrain the ionized gas density (n_ion), the advection speed (v_adv), and the magnetic field strength (B) (see Figure 7). For the γ-rays, fits to the eight data points (excluding upper limits) were used to place additional constraints on the advection speed and the magnetic field strength.

Zoom In Zoom Out Reset image size

In looking for the joint best-fit model, we found that for the radio spectrum there is a degeneracy between magnetic field strength and wind speed such that as wind speed increases, an increase in magnetic field strength allows a similarly good fit (see Figure 7). We explain the reasons for this degeneracy below. For the γ-ray spectrum we also find a valley of solutions (several minima as opposed to a single absolute minimum). If we restrict the possible γ-ray models to those deemed to have acceptable radio spectra, we find that the χ² tests for the γ-ray spectrum limits the tested wind speeds to only a few.

The best-fit solution has a χ² value of χ² = 22.6 for the radio spectrum and χ² = 9.6 for the γ-ray spectrum. It occurs for an ionized gas density of 100 cm⁻³, a magnetic field strength of B = 275 μG, and a wind speed of v_adv = 500 km s⁻¹ (Table 3). The models within 3σ have magnetic field strengths of B = 225–350 μG, advection speeds of v_adv = 300–700 km s⁻¹, and ionized gas densities of n_ion = 50–250 cm⁻³. Figure 8 shows all radio spectra for models within 3σ of the minimum for a single spectral index and molecular gas mass. Of the 16,000+ models tested, we found 666 radio models within 3σ while only 89 had both radio and γ-ray models within 1σ.

Zoom In Zoom Out Reset image size

Table 3. Best-fit Model Parameters

Physical Parameters	Best-fit Value
Magnetic field strength (B)	275 μG
Advection (wind) speed (vadv)	500 km s−1
Ionized gas density (nion)	100 cm−3
Spectral index (p)	2.1
Molecular gas mass (Mmol)	4 × 108 M☉

Note. Results for $\chi ^{2}_{{\rm radio}}$ = 22.6, $\chi ^{2}_{\gamma }$ = 9.6.

Download table as:  ASCIITypeset image

The best-fit γ-ray model can be seen in Figure 5 and the corresponding radio model can be seen in Figure 6. The radio spectrum agrees to within the errors for 10 of 11 data points, while the γ-ray spectrum agrees with 6 of 8 data points. It can be clearly seen in Figure 6 that the radio model does not agree with the 92 GHz data point. The flux at 92 GHz has been previously shown to be mostly due to free–free emission (Williams & Bower 2010). Because we have limited the ionized gas to a single density, we achieve the turnover at low frequencies due to free–free absorption but not the required flux at high frequencies due to free–free emission. To achieve both the observed free–free absorption and emission would require another level of detail to the model. Thus, in maintaining the simplicity of the model, we are not able to produce the required free–free flux at high frequencies.

4.4.2. Wind Speed and Magnetic Field Strength

4.4.3. Gas Mass and Spectral Index

One of the main questions we seek to answer in this paper is whether or not M82 is a calorimeter. We define a calorimeter such that cosmic-ray proton or electron energy losses dominate over advective or diffusive losses. In the case of the cosmic-ray electrons, it is clear that the galaxy is an electron calorimeter as energy losses dominate over advective losses for all wind speeds tested (see Figure 1).

However, it is more difficult to decide exactly when a galaxy is no longer a proton calorimeter. In the case of the cosmic-ray protons, the advection timescale is about equal to the peak of the energy loss lifetime (see Figure 1) meaning that the galaxy is a 50% calorimeter. How calorimetric a galaxy is depends on not only on the wind speed, which determines the advection timescale, and the average density, which determines the energy loss timescale, but also on the cosmic ray acceleration efficiency. Throughout the paper, we have assumed cosmic ray acceleration efficiency (η) of 10%. To investigate the impact of acceleration efficiency, we tested two alternative values: a decreased efficiency of η = 4% and an increased efficiency of η = 20% (that could also be achieved with an efficiency of 10% and an supernova rate of 0.14 yr⁻¹).

For the 4% case, the resulting radio models consistently required higher magnetic field strengths due to a lack of electrons because of the decrease in the primary electron population and because of a decrease in the secondary electron/positron population because of the decrease in primary protons. The best models required wind speeds on the order of 0–100 km s⁻¹ (see Figure 9). Wind speeds of this magnitude give advection timescales that are longer than the proton energy loss timescales, resulting in an 80%–100% calorimeter. The best fit had a χ² value nearly a factor of two larger than the best fit for an efficiency of η = 10% for both the γ-ray and radio spectra. For all three acceleration efficiencies, we found 317 of the 49,000+ models to be within 1σ for both the radio and γ-ray spectra. Of those 317 models, only 45 had an acceleration efficiency of 4%. Based on these results and the extremely low wind speeds, we can effectively rule out the 4% case.

Zoom In Zoom Out Reset image size

Conversely, for the 20% efficiency case, we found that significantly higher wind speeds were required to achieve good fits because of an increase in the overall electron/positron population. The best models required wind speeds on the order of ∼1500 km s⁻¹ (see Figure 9). While models for an efficiency of 10% agree with measurements of wind speeds in the optical (Shopbell & Bland-Hawthorn 1998), wind speeds on the order of 1500 km s⁻¹ agree with X-ray measurements by Strickland & Heckman (2009).

Comparing the resulting advection timescales to the proton energy loss timescales gives a 35% calorimeter. We find that for this higher acceleration efficiency, the minimum χ² value for the radio spectrum is smaller than in the 10% case, while the χ² value for the γ-ray spectrum is just slightly higher than before (χ² = 10.8 versus χ² = 9.6). Thus, while we rule out the 4% acceleration efficiency, we consider the 20% case to be plausible.

Acceleration efficiency directly affects our initial cosmic ray population. However, it is not the only factor. For example, we use the lower bound of 0.07 SN yr⁻¹ for our supernova rate. If we scaled the supernova rate up and the acceleration efficiency down in equal measure, our initial population of cosmic rays would be unchanged. While we have selected one particular set of parameters for our initial cosmic ray population, it is simple to obtain the same population from a slightly different set of initial parameters and as such, we would again produce a degenerate relationship between magnetic field strength and wind speed. Thus, we would produce the same fits as presented above.

It remains a curious question as to why starburst galaxies lie on the same FIR–radio correlation as the normal, quiescent galaxies. The general explanation for the FIR–radio correlation in disk galaxies is the electron calorimeter theory by Voelk (1989). This paper postulates that both FIR and radio synchrotron emission are proportional to the supernova rate as FIR emission is essentially starlight reradiated by dust grains and radio synchrotron emission comes from energetic electrons accelerated by supernovae (Voelk 1989). As such, it is required that cosmic-ray electrons radiate away their energy within the galaxy to achieve the observed radio synchrotron emission. However, in dense environments such as starburst galaxies, there is a second source of electrons/positrons from proton–proton collisions.

The overall importance of secondary electrons/positrons is dependent upon the ratio of primary protons to primary electrons (N_p/N_e). The significance of secondary electron/positrons decreases as the ratio of primary protons to electrons decreases. In our model, we leave this ratio as an independent free parameter fixed at 50 based upon supernovae studies in the Milky Way. Schlickeiser (2002) makes a case for this parameter being model dependent such that the ratio depends on the masses of protons and electrons and the selected spectral index, N_p/N_e = (m_p/m_e)^{(p − 1)/2}. Our preferred spectral index is p = 2.1, which gives a ratio of N_p/N_e ∼ 60. As our free parameter selection is smaller than this value, the importance of secondary electrons/positrons would only increase by making the ratio model dependent.

Thus, while we find that M82 is an electron calorimeter, the overall picture is complex. It is true that all cosmic ray electrons are mainly confined to the starburst core. However, because the galaxy is not also a proton calorimeter, there is an inherent loss of any secondary electrons and positrons that would have resulted from proton–proton collisions in the core. Thus, any electrons produced in the starburst region (primary or secondary) are confined to the core but the total possible population of electrons/positrons is not equivalent to the confined population. This is one of the mechanisms responsible for keeping the galaxy on the FIR–radio correlation. Without some loss of protons to a advection in a galactic wind, there would be a larger population of electrons/positrons resulting in a higher radio flux.

Our approach to modeling cosmic rays in starburst galaxies builds on that taken in previous models (de Cea del Pozo et al. 2009b; Paglione & Abrahams 2012). However, we place an emphasis in our models on varying and constraining a variety of possible wind (advection) speeds rather than assuming a single, static wind speed (de Cea del Pozo et al. 2009b; Paglione & Abrahams 2012). Additionally, whereas some previous models have had more free parameters (Lacki et al. 2010), our approach reduces the varied parameters to a select few. All other parameters are chosen to be in agreement with observations. Thus, we compare large numbers of models and strive to produce fits to the radio and γ-ray spectra that are more intuitive and observationally informed.

Paglione & Abrahams (2012) favor higher mean gas densities, ranging from 100 to 1000 cm⁻³, and strong magnetic fields, ranging from 200 to 1000 μG. Their best solution resulted in a density of n = 600 cm⁻³ and a magnetic field strength of B = 450 μG. In contrast, we favor smaller magnetic field strengths (225 to 350 μG). However, though we test smaller densities (280 to 415 cm⁻³), our best fits have a comparable density of 550 cm⁻³ (see Table 1). De Cea del Pozo et al. (2009b) favor models with larger supernova rates (0.1–0.3 SN yr⁻¹) and use a significantly smaller average ISM density (180 cm⁻³). They produce models with magnetic fields strengths from 120 to 290 μG. This range is comparable to our own results.

Like our model, both de Cea del Pozo et al. (2009b) and Paglione & Abrahams (2012) incorporate free–free absorption into their radio models with spectra that turn over at ∼700 MHz and ∼200 MHz, respectively. Neither incorporates a halo synchrotron emission component to allow for comparison with low frequency radio data. Additionally, Paglione & Abrahams (2012) assume a constant convection loss timescale of 1 Myr which is sufficiently longer than the energy loss timescale for cosmic rays (at their preferred average density of 600 cm⁻³) that they find the galaxy to be calorimetric for both protons and electrons.

In addition to creating a model for cosmic rays in starburst galaxies, we address the question of calorimetry in starburst galaxies. Calorimetry is also address by the recent paper (Ackermann et al. 2012). A semi-analytical formula which depends on the SFR and supernova acceleration efficiency of a galaxy is used to determine the calorimetric efficiency estimate, while we compare a fully calculated energy loss lifetime against the advection timescale for a variety of wind speeds (see Sections 3.1 and 5.1). Ackermann et al. (2012) find that galaxies with SFR ∼10 M_☉ yr⁻¹ have calorimetric efficiencies of 30% to 50%, a similar result to our own findings.