MULTIMESSENGER MODEL FOR THE STARBURST GALAXY M82

M82 is a near-starburst galaxy. It can be seen nearly edge-on (77°; Mayya et al. 2005) and has a gas content mostly concentrated in the inner 2 kpc. This galaxy presents a high luminosity both in the far-infrared (FIR) and X-ray domain (10⁴⁴ erg s⁻¹ and 10⁴⁰ erg s⁻¹, respectively; e.g., Ranalli et al. 2008). As part of the M81 group, M82 shows hints of an encounter with some of its members 1 Gyr ago, specifically with the dwarf galaxy NGC 3077 and the former M81, materialized in an intergalactic gas bridge of 20 kpc. Up to 10 kpc above the plane of M82, galactic superwinds can be detected. Following H i streamers, the external part of the disk (greater than 5 kpc) presents a warped form. On the other hand, the inner part (300 pc) harbors a starburst. Around this region, a molecular ring of 400 pc radius and an NIR bar of ∼1 kpc length can be detected (Telesco et al. 1991). From the tips of this bar, two symmetrical arms emerge (Mayya et al. 2005). This led to a change in the morphological classification. Previously, an irregular shape was assumed due to optical appearance: bright star-forming knots interspersed by dusty filaments (O'Connell & Mangano 1978); now it seems more likely that M82 is an SBc galaxy. Freedman et al. (1994) established a distance for M81 of 3.63 ± 0.34 Mpc, thanks to the discovery of new Cepheids by the Hubble Space Telescope (HST) in this galaxy. A few years later, Sakai & Madore (1999) found a distance of 3.9 (±0.3)_random(±0.3)_systematic Mpc for M82, based on the detection of the red giant branch stars using HST photometry. Both distances are consistent within errors; the latter is the one used in this work. The starburst region is located in the inner part of the galaxy (radius of ∼300 pc, height of ∼200 pc; e.g., see Völk et al. 1996 and references therein; Mayya et al. 2006); meanwhile the rest of the emission extends to a thin disk up to 7 kpc (see, e.g., Persic et al. 2008 and references therein). A significant amount of molecular material has been found in the central region, where most of the starburst activity (see the discussion on the supernova (SN) explosion rate below) is located. The content of the gas mass in the whole galaxy ranges from 2.9 × 10⁹ M_☉ (H i) by Crutcher et al. (1978) and 2.9 × 10⁸ M_☉ (H₂) by Young & Scoville (1984) to the more recent results found by Casasola et al. (2004) of 7 × 10⁸ M_☉ (H i) and 1.8 × 10⁹ M_☉ (H₂), where different assumptions for distance and the normalization of luminosity were made. Weiss et al. (2001) report around 2 × 10⁸ M_☉ (H₂) for the starburst region, using separate methods for the determination of CO and H₂ densities (see this paper for further details). This value essentially agrees with those estimated from 450 μm dust continuum measurements (Smith et al. 1991) and from CO(2 → 1) intensities (Wild et al. 1992), and is, therefore, the one used in the determination of the uniform density for the model of this paper (∼180 cm⁻³).

With the advent of Fermi, the upgrade of MAGIC to a stereoscopic observatory, and with VERITAS online, it is important to have the most detailed model of M82 to be able to feedback from observations our knowledge of the cosmic-ray population, and the physical environment of nearby starbursts. In addition, with the recent ongoing discussion about the starbursts' contribution to the neutrino background, it is important to have consistent and detailed models of neutrino production in specific galaxies, where the assumptions made in previous studies could be studied in detail and observational detectability with, e.g., km³ observatories could be assessed.

The aim of this study is to obtain multifrequency/multimessenger predictions for the photon and neutrino emission coming from the central region of M82, where the inner starburst is located. We seek a consistent model, where the different components of the emission (from radio to TeV photons and neutrinos) can all be tracked to one and the same original cosmic-ray population, and result as a consequence of all electromagnetic and hadronic channels from the primary and subsequently produced secondary particles. With this model at hand, we explore a range of uncertainties in the parameters, which may correspondingly lift up or lower the high-energy end of the spectrum. In order to perform such a study we have updated the code ${\cal Q}$-diffuse (Torres 2004) in several ways. ${\cal Q}$-diffuse solves the diffusion-loss equation for electrons and protons and finds the steady state distribution for these particles subject to a complete set of losses in the interstellar medium (ISM). Subsequently, it computes secondaries from hadronic interactions (neutral and charged pions) and Coulomb processes (electrons), and gives account of the radiation or decay products that these particles produce. Secondary particles (photons, muons, neutrinos, electrons, and positrons) that are in turn produced by pion decay are calculated too. Additional pieces of the code compute the dust emissivity, and the IR–FIR photon density, which is consistently used both as a target for inverse Compton scattering and to model the radiation at lower frequencies. For radio photons, we compute, using the steady distribution of electrons, synchrotron radiation, free–free emission, and absorption. Finally, opacities to γγ and γZ processes are computed, as well as absorbed γ-ray fluxes, using the radiation transport equation. Our implementation of ${\cal Q}$-diffuse includes several upgrades. Some are strictly technical, allowing for a more versatile and automatic input–output interaction; some are physical: we have developed neutrino-production subroutines for all neutrino channels in the decay of positively and negatively charged pions and make use of a recent parameterization of pp interactions by Kelner et al. (2006). The flow of the ${\cal Q}$-diffuse code is shown in Figure 1. Several loops can be distinguished there; they are used in order to determine a self-consistent set of parameters such that all predictions for the different bands of the electromagnetic spectrum are consistent with data.

Zoom In Zoom Out Reset image size

Previous studies of diffuse high-energy emission, and of electron and positron production, with different levels of detail and aims go back to the early years of γ-ray astronomy. A summary of these first efforts is summarized by Fazio (1967), Ginzburg & Syrovatskii (1964), and Ramaty & Lingenfelter (1966), then followed by Maraschi et al. (1968), and Stecker (1977), among many others. Secondary particle computations have a similarly long history: see, e.g., Stecker (1969) and Orth & Buffington (1976). Relatively more recent efforts include Drury et al. (1994), Moskalenko & Strong (1998), Strong & Moskalenko (1998), Markoff et al. (1999), and Fatuzzo & Melia (2003). Here, the general ideas used by Paglione et al. (1996) and Blom et al. (1999), when modeling nearby starburst galaxies, are followed. These, in turn, closely track Brown & Marscher's (1977) and Marscher & Brown's (1978), regarding their studies of close molecular clouds.

We emphasize here that this approach lacks detail of modeling of single sources: we deal not with modeling isolated SN explosions or SN remnants, nor with the issues of them exploding within the wind bubbles of their progenitors and the concurrent losses in these environments, but rather on average properties of the whole starburst region, seen as a cosmic-ray (CR) injector. To mix these two aims would require an effort that is beyond both current observational capabilities of what is known in M82 and what is possible for this work: essentially, to make quantitative estimates in this way one would need to understand the features of each and all the SN explosions in M82, model them in detail individually, and from that build (i.e., from stacking of a source-by-source modeling) an averaged scenario. This has not been achieved even for our own Galaxy. Similarly, our approach to the modeling of the superwind is also simple in order to encompass it with the rest of the components of the model (see Everett et al. 2008 for a more detailed analysis of the properties of the Galactic superwind; see also Gallagher & Smith 2005).

The assumption of a uniform distribution of accelerators in the inner region of the starburst justifies to have (and to compute it via solving the full diffusion-loss equation) the existence of a steady state distribution therein. The limited size and compactness of the starburst region is compatible with a steady state electron distribution in that central region. The steady state particle distribution is computed as the result of an injection distribution being subject to losses and secondary production in the ISM. At sufficiently high energies, the injection proton emissivity is assumed to have the following form as a function of proton kinetic energies: Q_inj(E_p, kin) = K(E_p, kin/GeV)^−pexp(−E_p, kin/E_p, cut), where p is a power-law index, K is a normalization constant, and E_p, cut is an energy cutoff in the accelerated particles injected (assumed as 100 TeV, with a range explored below). Units are [Q]= GeV⁻¹ cm⁻³ s⁻¹. To get a numerical prior for the normalization K is assumed to come from the total power transferred by SNe into kinetic energy of CRs within a given volume; thus, the SN rate is essential to fix how high is the cosmic-ray sea $\int _{E_{\rm p,\,kin,\, min}}^{E_{\rm p,\, kin,\, max}} Q_{\rm inj}(E_{\rm p,\, kin}) E_{\rm p,\, kin} dE_{\rm p,\, kin} \equiv { \eta {\cal P} {\cal R} }/ {V}$. ${\cal R}$ is the rate of SN explosions, V being its volume, and η the transferred fraction of the SN explosion power (${\cal P} \sim 10^{51}$ erg) into CRs. The average rate of power transfer is assumed to be 10%. Our assumption of the p value is not a priori; it is rather an a posteriori choice made by what it seems to be a slightly better fit to all multiwavelength data and acts as an average description of the injected cosmic-ray sea in the star-forming region, where multiple nearby shocks could contribute.

As noted by Domingo-Santamaría & Torres (2005), at low energies the distribution of cosmic rays is probably flatter, e.g., it would be given by Equation (6) of Bell (1978), correspondingly normalized. We have numerically verified that to neglect this difference at low energy does not produce any important change in the computation of secondaries, and especially on γ-rays at the energies of interest. Here, we assume, e.g., as in Berezhko et al. (2006), that electrons are injected into the acceleration process at the shock fronts. We choose the electron injection rate such that the electron–proton ratio is a constant to be determined from the synchrotron observations. The diffusive transport equation takes care of the changes produced after injection onto this distribution.

The diffusion-loss equation is given by (see, e.g., Longair 1994, p. 279; Ginzburg & Syrovatskii 1964, p. 296) −D▽²N(E) + N(E)/τ(E) − d/dE[b(E)N(E)] − Q(E) = −∂N(E)/∂t. In this equation, D is the scalar diffusion coefficient, Q(E) represents the source term appropriate to the production of particles with energy E, τ(E) stands for the confinement timescale, N(E) is the distribution of particles with energies in the range E and E + dE per unit volume, and b(E) = −(dE/dt) is the rate of loss of energy. The functions b(E), τ(E), and Q(E) depend on the kind of particles considered. In the steady state, ∂N(E)/∂t = 0, and when the spatial dependence is considered to be irrelevant D▽²N(E) = 0. Then the diffusion-loss equation can be numerically solved by using the corresponding Green function G, such that for any given source function, or emissivity, Q(E), the solution is $N(E) = \int _E^{E_{\rm max}} d{E^\prime } Q({E^\prime }) G(E,{E^\prime }).$ The confinement timescale will take into account that particles can be diffusing away, can be carried away by the collective effect of stellar winds and SNe, or can be affected by pion production (for CRs). Pion losses (which are catastrophic, since the inelasticity of the collision is about 50%) produce a loss timescale τ⁻¹_pp = (dE/dt)^pion/E (see, e.g., Mannheim & Schlickeiser 1994). Thus, in general, for energies higher than the pion production threshold τ⁻¹(E) = τ⁻¹_D + τ_c⁻¹ + τ⁻¹_pp. τ_c, the convective timescale, is ∼R/V, where V is the collective wind velocity. Strickland et al. (1997) found a value of 600 km s⁻¹ for the wind velocity, and we use this value below. However, we studied the response of the model to an arbitrary doubling or halving of this value and found that the uncertainty so introduced is smaller than those produced by other parameters. Finally, assuming an homogeneous distribution of accelerators in the central hundreds of parsecs of M82, the characteristic escape time is that from the homogeneous diffusion model (Berezinskii et al. 1990, p. 50–52 and 78) τ_D = R²/(2D(E)) = τ₀/(β(E/GeV)^μ), where β is the velocity of the particle in units of c, R is the spatial extent of the region from where particles diffuse away, and D(E) is the energy-dependent diffusion coefficient, whose dependence is assumed ∝E^μ, with μ ∼ 0.5, and τ₀ is the characteristic diffusive escape time at ∼1 GeV. Results that follow use τ₀ = 10 Myr, but a range was explored to judge uncertainties. For electrons, the total rate of energy loss considered is given by the sum of that involving ionization, inverse Compton scattering, adiabatic, bremsstrahlung, and synchrotron radiation. The Klein–Nishina cross section is used.

Zoom In Zoom Out Reset image size

A natural comparison can be made between this approach and previous studies of the central starburst region of M82, in particular Akyüz et al. (1991), Völk et al. (1996), Paglione et al. (1996), and Persic et al. (2008). The earlier work of Akyüz et al. focused on γ-ray emission at the GeV regime (integral flux above 100 MeV) and basically provides a simple order of magnitude estimation with data available at that time. Völk et al. also did not provide a multifrequency model, but rather consider only neutral pion decay from protons to obtain an order of magnitude estimation of the flux at high and very high gamma-rays. This approach was also taken (with minor differences) by Pavlidou & Fields (2001) and Torres et al. (2004). These computations (having a different aim than the one here, which is providing a consistent model all along the electromagnetic spectrum) do not include a computation of secondaries and their emission at low energies. The work by Paglione et al. (1996) is more complex, and it was discussed in detail, in particular when compared with results by ${\cal Q}$-diffuse on NGC 253, by Torres (2004) and Domingo-Santamaría & Torres (2005), where we refer the reader. There are several different physical and methodological considerations embedded in their model than the ones assumed here, in addition to not providing high-energy (TeV) predictions. None of this work considers neutrino emission. The most recent study on M82, from Persic et al. (2008), also differs from ours in some key parameters and, most importantly, in the method and assumptions for the modeling. Among the physical parameters for which values are given by Persic et al., the distance (3.6 Mpc, which is actually the one to M81), the dust emissivity index (σ = 1), the proton to electron primary ratio (N_p/N_e = 200), the magnetic field (180 μG), and the slope of the primary injection spectrum (2.4) are all different from the ones used or derived in the present study. The method by which these (latter) values are fixed by Persic et al. differs from ours. For instance, Persic et al. obtained the magnetic field by assuming equipartition with primary particles, while we find the magnetic field by a multifrequency analysis that takes into account also secondary particles produced: as we show in Figure 2, left panel, and was found also in the cases of other highly dense environments, the steady population of electrons that would result after solving the diffusion loss equation injecting only primary and only secondary (i.e., those electrons and positrons resulting from knock-on and pion decay in the inner region of M82) particles are comparable, so that it is not a safe assumption to consider that the primary population dominates in computing the equipartition field. In the numerical treatment by Persic et al. (2008), the path is thus inverse to ours. They seemingly fix the normalization of the proton spectrum starting from an assumed N_p/N_e factor from equipartition, with the normalization for electrons N_e in turn obtained from the radio spectrum. In our case, we start from protons, with the SN explosions injecting primaries which are left to interact by all possible processes producing electrons, which in turn radiate synchrotron. It is this (not at all negligible) secondary electrons contribution that is taken into account when computing B and N_p/N_e in order for the data to be consistently matched at all frequencies (see Figure 1 for further clarification). Neither equipartition is a priori assumed, nor is the primary electron population a priori fixed from the radio spectrum in our approach. Finally, neutrino emission is scaled from gamma-ray fluxes by Persic et al., while it is computed by use of correspondingly parameterized cross sections up to tertiary particles interactions in this work.

The steady population of protons and electrons of the starburst galaxy M82 which results from our model is shown in Figure 2, right, where a power-law index p = 2.1 for the injection of primary relativistic hadrons is used and, to fix numerics, an exponential cutoff is set to 100 TeV (no significant difference is observed even at energies as high as 10 TeV if it is changed to half this value). A range of values is studied for the normalization factor for the hadronic injection spectrum too, as it depends on other sensitive parameters, subject to uncertainties. For instance, the rate of SN explosions per year was assumed to be 0.3 SN yr⁻¹ in earlier works, but recently it is more commonly referred to as ∼0.1 SN yr⁻¹. In the present work, both a high and a low value for the SN rate were studied, revealing a range of uncertainties as can be seen in Figure 4, where we have also considered other exponential cutoffs.

Indeed, the discussion about the SNe explosion rate is ongoing. The highest value of 0.3 SN yr⁻¹ is used by, e.g., Völk et al. (1996) and other authors thereafter, and it is based on Kronberg et al. (1985) who compile different estimations, starting from 0.1 SN yr⁻¹ by Kronberg & Wilkinson (1975), based on the total nonthermal emission, going up to 0.16 yr⁻¹ (based only on IR excess), then 0.2–0.3 by Kronberg & Sramek (1985) based on direct monitoring of variability of discrete sources, and still going up to 0.3 SN yr⁻¹ (Rieke et al. 1980), based on estimating the number of new radio sources. But critics of Kronberg & Sramek (1985) and in general to high values of the SN remnant rate have arisen based mainly on the fact that the rate of detection of new radio sources does not correspond to them (see, e.g., de Grijs et al. 2001; McLeod et al. 1993; Bartel et al. 1987). Lower values for the rates (0.07–0.08 SN yr⁻¹) have also recently been favored by Fenech et al. (2008) who made an 8-day deep MERLIN radio imaging of supernova remnants (SNRs) in the starburst and use two different methods (on one side, they assumed that the SNRs in M82 that are more radio luminous than Cas A are younger than it; on the other, assuming that the SNRs are in free expansion and using the measured N(<D) − D plot) to compute the explosion rate.

The primary electron population is taken to be proportional to the injection proton spectrum, scaled by the ratio between protons and electrons, N_p/N_e. The value is set later on, together with others, e.g., the magnetic field, in a recursive application of our model, as a fit to radio data. Together with the primary electron population, all secondary electrons produced are taken into account when computing the steady state. As a comparison, in Figure 2, right, we also plot the steady particle populations present in NGC 253 (from Domingo-Santamaría & Torres 2005). Both curves present a similar behavior and level as those corresponding to M82, emphasizing the similarity of the initial characteristics of these two starbursts as possible high-energy sources. With the steady electron population determined, we proceed to compute a multifrequency spectrum from radio up, and compare with experimental data.

Figure 3 shows the best match from our model: a multifrequency spectrum is overplotted to previous radio and IR data (see caption for details). At the lowest frequencies, radio data are no longer reduced to the central starburst region. Due to the angular resolution of the observing instrument, emission at this frequency comes from the whole galaxy and cannot be separated from that coming only from the inner region, so our model is just a component of the total radiation below a few 10⁸ Hz. Along the part of the spectrum in which synchrotron emission dominates, several parameters can be determined through repeated iterations, e.g., magnetic field, EM, slope of proton injection, and N_p/N_e ratio. In the range of uncertainty mentioned above, the lowest value for the magnetic field, 120 μG, corresponds to the highest steady distribution (i.e., with the highest SNe explosion rate assumed, and a value N_p/N_e = 30). Inversely, a low steady distribution implies the need for greater losses by synchrotron in order to match the data, so leading the magnetic field up to 290 μG. The general fit to observations is quite good in each case. These values of magnetic field are in agreement with previous results for M82 (Völk et al. 1989) and are also similar to the ones found for the disk of Arp 220 (Torres 2004) or NGC 253 (Domingo-Santamaría & Torres 2005), as well as compatible with measurements in molecular clouds (Crutcher 1988, 1994, 1999). On the other hand, the slope of the proton injection p is determined also here by matching the corresponding slope for radio emission. The range of uncertainties does not seem to affect it significantly (see below).

Zoom In Zoom Out Reset image size

Tabatabaei et al. (2007) confirm, in a study of M33, that nonthermal emission is spatially associated with star-forming regions and propose that the higher the (synchrotron photon) energy the more important the contribution of localized structures (SNRs) may be. Although this may be true, the spatial correlation by itself is not really a proof of it, since star-forming regions also have an enhanced cosmic-ray sea, and thus an enhanced production of secondary electrons which emit nonthermally. In fact, Bressan et al. (2002) have shown that the contribution of radio SNRs to the nonthermal emission cannot be dominant in galaxies, finding that in our Galaxy in particular it cannot be more than about 6%. In any case, this is an uncertainty in determining the magnetic field in our model that must be emphasized: we are considering that most of the radio emission comes from cosmic-ray electrons that have been either injected by SNRs (primaries) or other accelerators or produced in cosmic-ray interactions (secondaries from knock-on, charged pion decays). If there is a contribution to the total nonthermal radio signal made by unresolved SNRs, this would diminish the estimation of the average magnetic field in the starburst. This distinction cannot, however, be made with the data now at hand on M82.

Regarding the FIR emission, starburst data generally require a dust emissivity law, ν^σB(E, T), where typically σ = 1.5 and B(E, T) is the Planck function. This graybody peaks at ∼45 K and has a luminosity of 4 × 10¹⁰ L_☉ (Telesco & Harper 1980, corrected by distance). However, at higher frequencies, simple graybody emission cannot explain observations, and an excess appears at NIR wavelengths. Lindner (2003) proposed a dust cloud model in which an envelope of dust grains surrounds a luminous source. Each concentric shell decreases its density with increasing radius and has its own temperature and associated flux. As a simpler but still accurate approximation, we show here that the addition of just one secondary graybody is enough to fit the data well for our aims, the temperature is warmer (T≃ 200 K) than the main one, and has a smaller luminosity of 7 × 10⁹ L_☉. Small changes in the warm graybody luminosity or maximum of the peak are possible with results similar to the ones presented here, that we have explored. The fit to the IR data is quite good, as can be seen in Figure 3.

The model discussed in this work yields inverse Compton fluxes of just a few percent or less than the upper limits at X-ray energies (see, e.g., Persic et al. 2008 and references therein). The diffuse emission is overwhelmed by emission from compact objects.

Figure 4 shows the γ-ray differential flux that results from this model coming from the central region of M82, together with the sensitivity curves of possible instruments observing it. The goodness of this prediction is that it is consistent with the whole of the rest of the electromagnetic spectrum: the particles generating this emission, both hadrons and leptons, are the same ones generating the radiation seen at lower energies. The model confirms that there is not an expectation for detection by EGRET (consistent with the EGRET upper limit by Torres et al. 2004, and the stacking analysis of EGRET data by Cillis et al. 2005) but predicts M82 as a source for Fermi. The total flux estimations at high and very high energies are as follows: for E > 100 MeV, 2.6 ×  10⁻⁸ (8.3 ×  10⁻⁹) cm⁻² s⁻¹; for E > 100 GeV, 8.8 ×  10⁻¹² (2.8 × 10⁻¹²) cm⁻² s⁻¹, with the values parentheses representing the results obtained with the lower energy cutoff and lower SN explosion rate. Separate contributions are plotted from each γ-ray channel: neutral pion decay, bremsstrahlung and inverse Compton (against CMB, FIR and NIR photon densities). Finally, Figure 4 also presents neutrino fluxes coming from the inner part of the starburst galaxy. The separate contribution of each neutrino process can be seen there, together with the total flux. As in the previous case, values in parentheses represent the results obtained with the lower energy cutoff and lower SN explosion rate. The total flux estimations at high and very high energies are as follows: for E > 100 GeV, 1.2 × 10⁻¹¹ (3.9 × 10⁻¹²) cm⁻² s⁻¹; for E > 1 TeV 3.8 × 10⁻¹³ (9.9 × 10⁻¹⁴) cm⁻² s⁻¹.

Zoom In Zoom Out Reset image size

In the light of our results, only if the highest end of the predictions happens to be a realistic representation of the galaxy, the MAGIC telescope could detect it above 300 GeV with 5σ in 50 hr. Although the time required by MAGIC alone would be unrealistic in order to obtain a detection for the lowest end of the predictions, with the upcoming arrival of MAGIC II the time to expend on this source could be affordable. The actual estimations for MAGIC II sensitivity are a factor of 2–3 better than for MAGIC I, so the flux coming from the starburst galaxy would be detected within 50 hr. Similar estimations would apply for the VERITAS array. But indeed, the instrument of choice for this source would be the forthcoming CTA (now in Design Study). Presenting an instrument acceptance extending both to the lower and higher energy ends of that of the previous telescopes, with a 1 TeV sensitivity 1 order of magnitude or more better than MAGIC I, M82 should, if this model is a realistic representation of it, be a bright source for CTA. With the latter, studies on the possible cutoff in the proton spectra could be made, although M82 would still look like a pointlike source.

We have explored the dependence of these results on the initial injection slope p (enforcing it to be different) finding that, within a reasonable range, this is not relevant. In Figure 5, we show this by analyzing the predictions of our model for a −2.3 average injection slope. We show both the multiwavelength predictions in radio-IR and the gamma-ray emission. We see that no relevant differences (within uncertainties) are found in the parameters that give rise to these curves.

Zoom In Zoom Out Reset image size

We have mentioned that care should be exercised with equipartition estimates, since after the work by Beck & Krause (2005), it was shown that the usual formula to compute the equipartition field (or the minimum-energy estimate of total magnetic field strengths from radio synchrotron intensities) may be of limited practical use. It is based on the ratio $\cal K$ of the total energies of cosmic-ray protons and electrons (which is generally computed without a full model of the system). In addition to other nontrivial technical problems mentioned by Beck & Krause, if energy losses of electrons are important, the number density between these particles increases with particle energy and the equipartition field may be underestimated significantly; the correct value can be computed only by constructing a model of gas density and cosmic-ray propagation. Beck & Krause (2005) already emphasize that starburst galaxies and regions of high star formation rate in the central regions and massive spiral arms of galaxies have high gas densities and mostly flat radio spectra, and nonthermal bremsstrahlung and other losses are important. Using the steady population of relativistic particles that we have found as a result of our model, we determine the energy content in such a component, and from there we impose equipartition in order to see what would be the magnetic field determined in this way. It results in 150 μG, which is close to (but not the same as) the estimation by Weaver et al. (2002). Further analysis of the details of the equipartition issue in starbursts will be presented elsewhere.

Loeb & Waxman (2006) have suggested that hadronic processes in starbursts can produce a large enough background of diffuse high-energy neutrinos to be observable with ICECUBE. They arrive at this conclusion (criticized by Stecker 2007) by making order-of-magnitude estimations under the assumptions: (1) that protons lose essentially all of their energy to pion production, and (2) that a lower limit to the energy loss rate of the protons can then be obtained from the synchrotron radio flux by assuming that all of the electrons (and positrons) which are radiating are from pion decay. Even when photodissociation in starbursts is a minor effect (Anchordoqui et al. 2008) assumption (1) is arguable, as Stecker emphasized, because particles are subject to diffusion and convection by winds in addition to pion losses, i.e., the galaxy is not a complete calorimeter, especially at the highest energies. Assumption (2) is arguable since synchrotron radiating electrons are not only the secondaries, but the primaries too. Figure 6 quantifies this for M82 and compares the contribution to the synchrotron radiation of the inner starburst if only primary or only secondary electrons are considered. From about E_e − m_e ∼ 10⁻¹ to ∼100 GeV, the secondary population of electrons slightly dominates (see Figure 2, left), with the primary population contributing significantly. It would then be an error to normalize the whole radio emission to the secondaries alone (or primaries alone) and then use this to fix the energy loss rate of protons, from which to estimate the neutrino emission. As far as we are aware, this paper presents the first full neutrino emission computation from a starburst galaxy which is self-consistent with the emission at all electromagnetic frequencies as part of the same model. Given the detail needed for such a description, a generalization of our M82 results to all starbursts is not recommended. We remark, however, that if this is done, and if we support it by basic assumptions on the starbursts' number density within the horizon, the neutrino diffuse emission would be below the Waxman–Bahcall limit.

Zoom In Zoom Out Reset image size

The sky-average sensitivity of the ICECUBE installation (IC-22, with 22 strings) at 90%C.L. to a generic E⁻² flux of ν_μ is 1.3(2.0) × 10⁻¹¹ TeV⁻¹ cm⁻² s⁻¹ (E/TeV)⁻² depending on the analysis (Bazo Alba et al. 2008). This is not enough to detect M82 directly (see Figure 4). In fact, using our estimation of the neutrino flux from M82 and following the comments made by Anchordoqui et al. (2004), we find than <2 events per year would be expected in the full ICECUBE facility, and a more definitive assessment of its sensitivity to such a signal will need to await further refinement of angular and energy resolutions, via improved knowledge of the detector response.

We acknowledge support by grants AYA 2006-00530, CSIC-PIE 200750I029, and AYA2008-01181-E/ESP. The work of E. de Cea del Pozo has been made under the auspice of an FPI Fellowship, grant BES-2007-15131. We acknowledge Luis Anchordoqui for discussions.