Cosmic-ray energy densities in star-forming galaxies

The energy density of cosmic ray protons in star forming galaxies can be estimated from π0decay γ-ray emission, synchrotron radio emission, and supernova rates. To galaxies for which these methods can be applied, the three methods yield consistent energy densities ranging from Up ∼ 0.1 − 1 eV cm−3 to Up ∼ 102 − 103 eV cm−3 in galaxies with low to high star-formation rates, respectively.


Introduction
Active star formation (SF) in galaxies leads to acceleration of protons and electrons (CRp, CRe) via the Fermi-I diffusive shock acceleration mechanism in supernova (SN) remnants (Fermi 1954; Ginzburg & Syrovatskii 1964;Bell 1978;Protheroe & Clay 2004).Timescales of starburst (SB) activity in galaxies are comparable to galactic dynamical timescales, τ SB ∼ τ dyn ∼ 10 8 yr.On the other hand, in a SB region the characteristic timescales for protons to gain energy (by acceleration) and to lose it by collisions with interstellar gas particles (leading to heating and pion production) are typically much shorter than τ SB (e.g., Persic & Rephaeli 2010).
The respective lengths of these timescales suggest that in galaxies, during typical episodes of SF, particle distributions are likely to attain asymptotic steady-state form.Additionally, under virial equilibrium, CR and magnetic fields are in an approximate minimum-energy state; this is equivalent to having CR and magnetic fields in (approximate) energy equilibrium (e.g., Longair 1994).
The equipartition assumption enables deduction of the CRp energy density in star forming galaxies (SFGs), U p (the main contribution to the particle energy density), indirectly from the electron energy density which can be estimated from radio synchrotron measurements if the protonto-electron (p/e) energy density ratio is known or can be reliably estimated.
Another way to derive U p is based on measuring the GeV-TeV γ-ray emission, which is largely from CRp interactions with ambient gas that yield neutral (π 0 ) and charged pions; π 0 decays into γ-rays.Thus, γ-ray measurements allow a direct estimate of U p ; this is now possible given the recent detection of γ-ray emission from SFGs.
a e-mail: persic@oats.inaf.itb e-mail: yoelr@wise.tau.ac.ilAlso, since energy-loss timescales are shorter than SF timescales, U p can be estimated from the observed rate of core-collapse SN and the deduced CRp characteristic residence time in the galactic disk, given a realistic estimate of the fraction of SN kinetic energy that is channeled into particle acceleration.
This paper describes measures of galactic CR energy densities based on the methods outlined above.Although not strictly independent, these methods are based on very different observationally deduced quantities -radio emission, γ-ray emission, and the rate of core-collapse SN.We find that the three methods yield consistent U p for a sample of 14 galaxies (and galactic environments) with widely varying levels of SF activity, from quiescent to intense starbursts.These are the only galaxies for which γ-ray data (and/or radio data and SN rates) are available.After reviewing the γ-ray, radio, and SN methods (sect.2,3,4), the respective results are discussed in section 5 and summarized in section 6.

Estimating U p from γ-ray emission
Detection of GeV-TeV emission from several SFGs provided the basis for observationally-based estimates of U p in the central SB region and across the galactic disk.In most SB galaxies, such as the two nearby ones M 82 and NGC 253, the central SB region with a radius of ∼ 200 pc, is identified as the main site of particle acceleration.Here, the injection particle spectrum is assumed to have the nonrelativistic strong-shock index q = 2.A theoretical N p /N e ratio, predicted from charge neutrality of the injected CRs, is likely to hold in this source region, as is the assumption of energy equipartition of particles and magnetic fields.
Adopting the convection-diffusion model for energetic electron and proton propagation and accounting for all the relevant hadronic and leptonic processes, the steady-state energy distributions of these particles, in both the nuclear SB region and across the full disk, can be determined in the context of a detailed numerical treatment (e.g., Paglione  ).The procedure is similar when SF is not undergoing a burst confined to the nuclear region but occurs throughout the whole disk.
The relevant electron energy loss processes are bremsstrahlung, Compton, and synchrotron, whereas for protons the main losses are γ-ray emission from π 0 decay following p-p collisions.Bremsstrahlung losses dominate at lower energies, whereas π 0 -decay losses dominate at higher energies; in the GeV-TeV region, emission is mainly from pp-induced π 0 decay (e.g., Rephaeli et al. 2010).
For a source with gas number density n g , proton energy density U p , and volume V, the integrated hadronic emission from pp-induced π 0 decay is L [q]  ≥E = V g [q] ≥E n g U p dV s −1 , with the integral emissivity g [η] ≥ in units of photon s −1 (H-atom) −1 (eV/cm 3 ) −1 (Drury et al. 1994).Therefore, U p can be determined, when L ≥ and n gas (r) are observationally known, and the particle steadystate energy distributions are numerically determined by solving the convection-diffusion model for CRe and CRp propagation.
For the two local SB galaxies M 82 and NGC 253, GeV-TeV fluxes and spectra agree with numerical modeling with U p ∼ 200 eV cm −3 in the central SB region (Acciari et  Published U p values based on GeV-TeV data modeling are reported in Table 2.

Estimating U p from radio emission
Determining U p from the measured radio flux clearly requires knowledge of the mean magnetic field in the emitting region, B, and a way to couple CRp and CRe.To overcome this (implied) indeterminacy, the assumption of field and particle energy equipartition is commonly made.In addition, the contribution of secondary electrons (from π − decay) 1 to the (steady state) electron density has to be included.(Throughout this section we follow Persic & Rephaeli 2014.) While the exact form of the particle steady-state spectral density does not generally have a single power-law form (e.g., Rephaeli 1979, Rephaeli & Persic 2015), the radiative yields are largely by protons and electrons with energies higher than a few Gev, for which Coulomb losses (which flatten the spectral density) are subdominant.In this limit, the total electron spectral density can be approximated by where the electron Lorentz factor γ is in the range γ 1 ≤ γ ≤ γ 2 , N e,0 is a normalization factor of the primary electrons, and χ is the secondary-to-primary electron ratio.The electron spectral index is q e ≥ 2, with the minimal value of 2 corresponding to the strong-shock limit of the Fermi-I acceleration mechanism.Outside the acceleration region q e > 2. Ignoring the contribution of lowenergy electrons with γ < γ 1 , the electron energy density is U e N e,0 (1 For a population of electrons described by Eq. ( 1), traversing a homogeneous magnetic field of strength B that permeates a region of (spherically equivalent) radius r s located at a distance d from the observer, and emitting a 5 GHz synchrotron radiation flux of f 5 Jy, the standard synchrotron formula (e.g., Blumenthal & Gould 1970) yields where quantities are expressed in c.g.s.units, the factor a qe is defined in the latter paper, and ψ 5 ≡ ( rs 0.1 kpc ) −3 ( d Mpc ) 2 ( f5 Jy ).From the above expressions we derive Since U e includes both primary and secondary electrons, the rough assumption that both populations can be characterized with nearly the same power-law index (see Persic & Rephaeli 2014) means that the primary electron energy density is U e /(1 + χ).Denoting the primary p/e energy density ratio by κ(q p , q e ), the proton energy density is U p κ(q p , q e ) U e /(1 + χ) .Since tight coupling is expected in the very dense environment of SB nuclei, particle and magnetic field energy densities can be assumed to be close to equipartition.If so, we can express the field in terms of the total particle energy density; since γ 2 >> γ 1 ,  2 we obtain In general, q e , q p , χ, γ 1 , and κ need to be specified in order to compute U p .
• The value of q e is deduced from the (nonthermal) radio spectral index, α, through q e = 2 α + 1.
• The proton spectral index has been measured to be q p = 2.2 in the nearby SBGs NGC 253, NGC 3034, and NGC 4945 (Ackermann et al. 2012).Suprathermal particles injected into a supernova shock have a power-law spectrum with index q = (R + 2)/(R − 1), where R is the shock compression ratio.Taking R 3.6 in the nuclear SB regions yields q p 2.2.
• The secondary-to-primary electron ratio, χ, depends on the effectiveness of creating electron-positron pairs from pp interactions of primary protons with ambient protons.Therefore, χ will depend on the injection p/e number ratio, ζ, that sets the number of primary protons per primary electron, on the pp cross section at ∼ > 1 TeV, and on the gas density.The inverse of the product of the latter two quantities is the characteristic mean free path of energetic protons, whose pp interactions yield π ± which decay into e ± with a branching ratio of 2/3.The probability for a proton to undergo pp interactions during its 3D random walk through a region of radius r s and density n, is χ = 2 3 ζ(q) √ 3 r s n σ pp (n).In a typical SB nucleus with r s = 0.2 kpc, n = 200 cm −3 , and q = 2.2, we get of χ ∼ 1, in agreement with results from numerical models (see Fig.  • For consistency with the assumed power-law form of the electron spectrum, we take the low-energy limit γ 1 to be the value of the Lorentz factor at which the sum of the Coulomb (or electronic excitation, in ionized gas) and bremsstrahlung loss rates equals the synchrotron-Compton loss rate.This is also based on the fact that even for the relatively high values of the magnetic field in SB nuclei, the measured radio emission (upon which our normalization of the electron density is based) samples electrons with γ > 10 3 .Equating the sum of the first two loss rates with the latter yields an estimate of γ 1 .In Fig. 1 we display the energy-loss rates expressed in Eqs. ( 18)-( 20) of Persic & Rephaeli (2014) for typical SB nuclei parameters.
• An approximate expression for κ at injection is κ(q) ( mp me ) (3−q)/2 in the relativistic limit.Thus, for given values of q e , q p , χ, γ 1 , and κ, the proton energy density can be estimated if the radio synchrotron flux, source size, and distance are known (see col. 3 in Table 1).

Estimating U p from SN rates
With SN shocks the likely sites of particle acceleration, the combination of core-collapse SN rate and CR proton residence time in the galactic disk, τ res , provide a basis for another independent estimate of the proton energy density.The value of τ res is determined by two timescales: (i) energy-loss timescale for p-p interactions, τ pp = (σ pp cn p ) −1 that, for protons with kinetic energy E ∼ > 10 TeV for which σ pp 50 mb, can be written as τ pp ∼ 2 × 10 5 ( np 100 cm −3 ) −1 yr; and (ii) CRp advection timescale, τ out , that characterizes proton escape out of the disk mid-plane region in a fast (v out ∼ 2500 km s  For a homogeneous distribution of SNe within a central SB region of radius r s the latter timescale is τ out = 3 × 10 4 ( rs 0.3 kpc ) ( vout 2500 km s −1 ) −1 yr.The overall residence time is then τ −1 res = τ −1 pp (n HI ) + τ −1 out (r s , v out ).During τ res , a number ν SN τ res of SN explode and deposit the kinetic energy of their ejecta, E ej = 10 51 erg per SN (Woosley & Weaver 1995), into the interstellar medium.Arguments based on the CR energy budget in the Galaxy and SN statistics yield an estimate of an efficiency factor of η ∼ 0.05 in converting kinetic energy to particle acceleration (e.g., Peng et al. 2016).Accordingly, In Fig. 3 we show values of U p vs. SN rates for the above galaxy sample; the scaling relation in Eq.( 5) is plotted as a thin dotted line which represents the explicit scaling U p = 100 ( νSN 0.2 yr −1 ) δ eV cm −3 with δ ∼ 1.5.that incorporates the reciprocal dependences τ res ∝ ν −α SN with α ∼ 2.5, and r s ∝ ν −β SN with β ∼ 1).

Uncertainties
The values of U p listed in Table 1 do not include substantial observational and modeling uncertainties.In this  1.Results deduced from γ-ray (radio) measurements are denoted by triangles (circles).The thin dotted line represents the scaling relation U p = 100 ( ν SN 0.2 yr −1 ) 1.5 eV cm −3 , that corresponds to Eq.( 5) (as specified in the text).
section we attempt to estimate the level of precision with which U p was determined based on only limited information on the errors in the various observational parameters.
The use of γ-ray measurements to estimate U p , the only direct method to measure U p , entails a substantial uncertainty due to the poorly known spectral slope and, especially for fainter and/or more distant galaxies, the density distribution of ambient gas.
Uncertainties in estimating U p from radio measurements are more tractable due to the fact that the quantities in Eq.( 4) are usually well determined for our sample galaxies, except for the p/e energy density ratio κ, for which a CRp spectral index, q p , need to be assumed.The uncertainty on the latter parameter, δq p 0.1, translates to a ∼50% uncertainty in κ, i.e. typically an uncertainty of ∼10% in U p as deduced from Eq.( 4).The effective source radii, r s , are deduced from high-resolution optical and radio data, so their values are relatively well determined, with typical uncertainties of a few 10%.Thus, all the galaxy quantities relevant to computing the CRp density are presumed to be measured at a precision level of ∼ 30%; accordingly, U p can be estimated to within an overall uncertainty of 50%.
A precise measurement of the actual rate of corecollapse SNe is obviously crucial essential for reliable estimates of U p .This is quite difficult especially in central SB regions where optical extinction is very strong.Moreover, radio counts of SN remnants require information on their ages in order to determine actual SN rates.For galaxies in our sample galaxies, available observational results indicate that ν SN are known to within a factor ∼2. Also, uncertainties in U p can also be due to uncertainties in τ res ; the latter mostly arise from the uncertainty in the fast wind velocity, which in turn is probably known to within a factor of ∼2.

Discussion
The three methods discussed here are not independent.The γ-ray and radio methods are related through the common dependence on the p/e ratio at injection, on the secondary-to-primary electron ratio, and on the assumption of particle-field equipartition.The SN method is not independent of the γ-ray method either, since both depend on the CRp residence time in the emission region -although, unlike the γ-ray and radio methods, it does not depend on the particle radiative yield but on the statistics of core-collapse SN.Also, the γ-ray, radio, and SN methods are not on equal footing.By these methods the value of U p is either measured, inferred, or estimated, respectively.This is due to (i) π 0 -decay γ-ray emission is the most robust measure of U p when the distribution of target gas is known and the particle propagation mode and energy losses are well known, whereas (ii) radio emission enables deduction of U p from the radio spectral distribution, and assumed p/e energy density ratio and particlefield energy equipartition; and (iii) assuming a SN origin for CR, SN statistics for a given (region of a) galaxy leads to an estimate of U p there.
Nonetheless, it seems that the three methods yield roughly consistent results.Specifically, as expected, there is a very significant difference between the very high CRp energy densities (U p ∼ O(10 2 ) eV cm −3 ), deduced for a central SB region as compared with the low values (U p ∼ O(10 −1 ) eV cm −3 ) for the rest of the galactic disk.
Finally, note that the distribution of data points in Fig. 3, with scaling U p ∝ ν 1.5  SN , is consistent (within uncertainties) with the observational relation L γ ∝ SFR 1.4  (Abdo et al. 2010b, for essentially this same sample of galaxies).This clearly supports a mostly pionic, SNpowered origin of the γ-ray emission in star forming galaxies (e.g., Pavlidou & Fields 2001).

Conclusion
For a sample of galaxies for which pointed GeV-TeV data are available (from Fermi/LAT) and from Cherenkov telescopes), we compared values of U p derived from the radio, γ-ray, and SN methods.Reasonable agreement among estimates based on the three methods is reached, showing clearly that there is a very substantial gap between the very high values (U p ∼ O(10 2 ) eV cm −3 ) in central SB regions such as Arp 220, M 82, NGC 253 and NGC 4945, and the low values (U p ∼ O(10 −1 ) eV cm −3 ) in the less active SF regions in the disks of, e.g., Local Group galaxies.
The results of this study extend the validity of our previous findings: (i) The Fermi-I acceleration mechanism, assumed to be at play in the environment of SN remnants, leads to acceleration timescales for CRp in galaxies such that particles and fields attain equilibrium over typical SF timescales, in agreement with observational evidence.(ii) CR energy densities and magnetic fields, inferred from radio data under the assumption of energy equipartition, can be used as proxies of the actual quantities that are measured directly only from γ-rays.This could be particularly useful in the case of distant (z ∼ > 1) galaxies, whose (unbeamed) γ-ray fluxes are too faint to be measured but whose radio fluxes are within reach of sub-mJy surveys (e.g., Tozzi et al. 2009).

Figure 1 .
Figure 1.Steady-state spectra of primary protons (dashed line), primary electrons (solid line), and secondary electrons (dashdotted line) in the central source region of NGC 253 (from Rephaeli et al. 2010).

Figure 3 .
Figure 3. Correlation between CRp energy density and SN rate for the sources in Table1.Results deduced from γ-ray (radio) measurements are denoted by triangles (circles).The thin dotted line represents the scaling relation U p = 100 ( ν SN 0.2 yr −1 ) 1.5 eV cm −3 , that corresponds to Eq.(5) (as specified in the text).

Table 1 .
Star-forming galaxies: CRp energy densities and SN rates.