1 Introduction

# Gamma-Ray Dominated Regions: Extending the Reach of Cosmic Ray Ionization in Starburst Environments

## Abstract

Cosmic rays are appealing as a source of ionization in starburst galaxies because of the great columns they can penetrate, but in the densest regions of starbursts, they may be stopped by pion production and ionization energy losses. I argue that gamma rays are the source of ionization in the deepest molecular clouds of dense starbursts, creating Gamma-Ray Dominated Regions (GRDRs). Gamma rays are not deflected by magnetic fields, have a luminosity up to that of the injected cosmic rays, and can easily penetrate column depths of before being attenuated by pair production. The ionization rates of GRDRs, , are much smaller than in cosmic ray dominated regions, but in the most extreme starbursts, they may still reach values comparable to those in Milky Way molecular clouds. The gas temperatures in GRDRs could be likewise low, if there is no additional heating from dust or turbulence, while at high densities, the kinetic temperature will approach the dust temperature. The ratio of ambipolar diffusion time to free-fall time inside GRDRs in dense starbursts is expected to be similar to those in Milky Way cores, suggesting star-formation can proceed normally in them. The high columns of GRDRs may be opaque even to millimeter wavelengths, complicating direct studies of them, but I argue that they could appear as molecular line shadows in nearby starbursts with ALMA. Since GRDRs are cold, their Jeans masses are not large, so that star-formation in GRDRs may have a normal or even bottom-heavy initial mass function.

galaxies:starburst – galaxies:ISM – gamma rays: galaxies – cosmic rays – galaxies: star formation

## 1. Introduction

With large star-formation rates often concentrated into small volumes, starburst galaxies are filled with the many kinds of high energy radiation that accompany star-formation, from ultraviolet photons to cosmic rays (CRs). This radiation can interact with the interstellar medium (ISM) of the abundant gas to heat or ionize it. The ionization of gas in star-forming environments in particular can drive chemistry in molecular clouds and sets the dynamics of magnetic fields that may in turn regulate star formation.

However, the ionization rate from this radiation is limited by two key factors: the luminosity of the radiation, and the regions within the starburst it can penetrate into. Energy conservation ultimately limits the rate atoms can be ionized by radiation; more radiation generated means more gas can be ionized, subject to the efficiency with which that radiation is converted into ionization. A powerful source of ionizing luminosity are extreme ultraviolet (EUV; ) or Lyman continuum photons from young O stars. The O stars create H II regions around themselves, where the gas is highly ionized. Indeed, the ionizing photon rate as measured by free-free radio luminosity is correlated with the star-formation rate (e.g., Condon, 1992; Murphy, 2011). The amount of ionized material in H II regions is large because of the high efficiency that EUV light ionizes surrounding gas: starburst galaxies are highly optically thick to EUV light (e.g., Leitherer et al., 1995; Hurwitz et al., 1997; Heckman et al., 2001). Almost all of the power in EUV radiation goes into ionizing H II regions or possibly into dust absorption.

However, the high optical depth of neutral gas to EUV light also poses a problem for the ionization of most of the ISM. Most of the EUV light is downgraded before it goes very far, so that H II regions are a small fraction of the volume of starburst galaxies, and something else must be responsible for ionizing the dense molecular gas in starburst galaxies. Far ultraviolet (FUV; ) light is not absorbed by neutral hydrogen, but it can dissociate molecular hydrogen and ionize carbon: thus, regions where FUV light penetrates into are photodissociation regions (PDRs; see the extensive review by Hollenbach & Tielens 1999). PDRs include much of the neutral gas in the Milky Way, but FUV light is quickly absorbed by dust (limiting both the power going into ionization and the penetration depth), and PDRs extend only to columns of (), smaller than typical columns through starbursts (; Kennicutt 1998; Hopkins et al. 2010).

The interiors of molecular clouds are ionized by more penetrative but less luminous sources of radiation. X-rays can traverse higher columns and can ionize neutral gas, creating X-ray Dominated Regions (XDRs; Maloney et al. 1996; Meijerink & Spaans 2005). In active galactic nuclei with high X-ray luminosities, X-rays can provide most of the ionizing luminosity, although this is concentrated towards the nucleus (Papadopoulos, 2010). Even in pure starburst galaxies, hard X-rays are emitted by high mass X-ray binaries, with a typical luminosity of , where is the bolometric (largely in the infrared) luminosity from star-formation (e.g., Ranalli et al., 2003; Persic et al., 2004). However, most galaxies are optically thin to hard X-rays, meaning most of that luminosity simply leaves the starburst before ionizing neutral gas. Furthermore, even hard X-rays are stopped once the columns become Compton thick (; ), which can occur in extreme starbursts like Arp 220 (Lehmer et al., 2010) and in molecular clouds. Even higher column densities can be attained in the parsec-scale gas around active galactic nuclei (Hopkins et al., 2012, and references therein).

An attractive candidate for the ionization of starbursts is cosmic rays 3, which ionize the interiors of Milky Way molecular clouds. The CRs regulate the temperature and ionization fraction in molecular clouds and cores, and because they pervade Milky Way molecular clouds, they set uniform conditions for star-formation throughout the Galaxy. CRs have many advantages as ionization agents in starbursts. (1) They are produced by star-formation at similar luminosities as X-rays, with (Thompson, Quataert, & Waxman, 2007). Since they trace star-formation, they naturally should be present in star-forming regions, unlike X-rays from AGNs, which are intense only near the nucleus itself. (2) They can traverse columns of up to at GeV energies, stopped mainly by inelastic collisions with ISM atoms that produce pions. At lower kinetic energies (less than a few hundred MeV), ionization losses also stop them, becoming rapidly more effective at smaller kinetic energies. They can therefore penetrate into regions of high dust extinctions, a big advantage over UV photons in the dusty environments of starbursts. (3) Galaxies are extreme scattering atmospheres for CRs, deflecting them with magnetic inhomogeneities. This allows the CR energy density to build up to very high levels (e.g., Socrates et al., 2008). Thus, Papadopoulos (2010) proposed that molecular gas in starburst galaxies are Cosmic Ray Dominated Regions (CRDRs; see also Suchkov et al. 1993).

But the high scattering optical depth to CRs is a double-edged sword. On the one hand, it means that each CR has many chances to ionize gas atoms as it scatters across a starburst, but on the other it means that a CR has many chances to be destroyed through pionic and ionization losses. The scattering atmosphere means that CRs see a much higher column that can stop them than photons do. Indeed, in the Milky Way, there is evidence that low energy () CR components contribute to high ionization rates in the diffuse ISM but are stopped by their own ionization losses (Indriolo et al., 2009). In dense starburst galaxies, protons may even lose all of their energy to pionic (and ionization) losses, in the “proton calorimeter” limit – with little of the CR luminosity escaping (Pohl, 1994; Loeb & Waxman, 2006; Thompson, Quataert, & Waxman, 2007; Lacki et al., 2010). So while starbursts with their dense gas may efficiently convert CR energy into ionization, CRs may be unable to reach heavily shielded molecular regions in starbursts. The GeV and TeV gamma-ray observations of M82 and NGC 253 (Acero et al., 2009; Acciari et al., 2009; Abdo et al., 2010) indicate that roughly of their CR power is lost to pionic emission (Lacki et al., 2011), much higher than in the Milky Way where the ratio is only a few percent (e.g., Strong et al., 2010), implying that pionic losses may be important for CRs in starbursts (Lacki et al., 2010).

An even more worrying problem with CR ionization in starbursts are the presence of winds (e.g., Chevalier & Clegg, 1985; Heckman, 2003): these carry CRs out of their starbursts (c.f., Suchkov et al., 1993), effectively wasting the power they could provide in ionizing the ISM; but unlike diffusion, winds do not let them penetrate deeper into molecular clouds (Crocker et al., 2011). The hard gamma-ray spectra of M82 and NGC 253 (Acero et al., 2009; Acciari et al., 2009; Abdo et al., 2010), combined with the evidence that , suggest that winds are important for CR transport in them, though denser starbursts like Arp 220 may be true proton calorimeters. However, Suchkov et al. 1993 and Papadopoulos 2010 did assume that winds limited the lifetimes of CRs in starbursts when calculating the CR energy densities.

To summarize the problem for ionizing radiation, radiation is optimally effective in uniformly ionizing regions when : in media with smaller , the radiation escapes without efficiently ionizing the gas, but in media with larger , the radiation ionizes only the skin of the region and leaves a shielded core. We therefore expect a hierarchy of ionization regions in starburst galaxies, with each step having lower luminosity but higher penetration depth. The H II regions and PDRs from UV light receive the most ionizing power but only constitute a small volume of starbursts; XDRs and CRDRs receive much less power but pervade denser phases of starbursts (Figure 1). Is there anything beyond CRDRs, which could guarantee the ionization of the molecular cores of starbursts?

I argue here that this ionization could come from gamma rays. High energy gamma rays are produced in starbursts by the pionic losses of protons. Additional leptonic gamma rays come from electrons and positrons (), some of which may be secondary made by pion production. Unlike CRs, gamma rays are not deflected by magnetic fields, and they are stopped only by Bethe-Heitler pair production at columns of (). In pair production, the gamma ray produces an electron and a positron in the electric field of an atomic nucleus () (e.g., Berestetskii et al., 1979). The result in dense columns is a cascade: the either ionize the ISM directly if they are low enough energy, or effectively downgrade the energy of the initial gamma ray when they cool by radiating. Thus it is possible to use the proton calorimetry of a starburst against it: the gamma rays can act as a “second wind” for CR ionization (c.f., Umebayashi & Nakano, 1981). The result is that even if CRs are stopped in the outer reaches of a molecular cloud, the inner regions can still be ionized as a Gamma-Ray Dominated Region (GRDR; Figure 1). Of course, in most of the starburst, gamma rays escape most of the time, with most of the energy “wasted” that way, but the key point is that gamma rays provide a low level of ionization that extends extremely deep into the molecular clouds of a starburst.

This paper will consider basic properties of GRDRs in extreme starburst environments. Using the energetics of gamma-ray production in starbursts, I estimates the characteristics of GRDRs in section 2, including the ionization rate (section 2.1), ionization fraction and ambipolar diffusion time (section 2.3), and temperature (section 2.4). I then estimate when the transition from CR to gamma-ray ionization occurs and demonstrate how the extreme scattering atmospheres of starburst can turn even clouds with relatively small columns into GRDRs (section 3). To more accurately calculate how the pair ionize the ISM, I construct a simple 1D model of a cascade in a starburst in section 4. I consider how GRDRs may appear observationally, including their indirect effects on stellar initial mass functions in section 5. I finally discuss some additional speculations in the Conclusions, section 6.

## 2. Basic Characteristics of GRDRs

### 2.1. Estimates of GRDR Ionization Rates

Suppose a starburst is a disk with radius and a midplane-to-edge scale height . It has a gamma ray luminosity of , giving a gamma-ray energy flux of . The volumetric ionization power that gamma rays deposit in a gas cloud that is optically thin to gamma rays () is

 Γγion=FγτγZfγion/ℓ=FγNHσγZfγion/ℓ=FγnHσγZfγion (1)

where is the gas column length, is the number density of hydrogen atoms4, is the column density, is the cross section for pair production, and is the fraction of energy in pair that goes into ionization losses. The quantity depends on the physical conditions within the cloud and the spectrum of the gamma rays; here I will estimate it as .

When a typical high energy proton ionizes an atom, the ejected electron will have a typical kinetic energy of (e.g., Cravens & Dalgarno, 1978); the CR proton therefore typically loses per ionization event. On the other hand, the initial ionized electron in turn ionizes, on average, additional electrons (e.g., Cravens & Dalgarno, 1978). Therefore, if ionizations by high energy are similar to those by high energy protons, each requires per ionization event, whether primary or secondary. Here I scale to 30 eV for simplicity.

The ionization rate from gamma rays is then simply :

 ζγ≈FγσγZfγionEion. (2)

Few starbursts are gamma-ray detected, but we can estimate the gamma-ray flux using the star-formation rate, the Schmidt law, and the time scale for pionic energy losses. The gamma-ray luminosity is

 Lγ≈fπLCR/3, (3)

where is the fraction of the CR power that goes into pionic losses: at high energies. The pionic loss time is roughly

 tπ≈50 Myr(⟨nH⟩cm−3)−1, (4)

where is the mean density of gas in the entire starburst (Mannheim & Schlickeiser, 1994); this is equivalent to saying pionic losses stop CR protons over a column . In terms of the mean column density of gas in the starburst , we have . The wind escape time is . Combining these timescales, I have:

 fπ≈⎡⎣1+0.16(vwind300 km %s−1)(⟨Σg⟩g cm−2)−1⎤⎦−1. (5)

The luminosity of injected CRs scales with the supernova rate (or more generally, the massive star formation rate). A CR luminosity of per supernova is expected. For a Salpeter IMF, we expect a supernova rate of (Thompson, Quataert, & Waxman, 2007). In turn, I can relate the star-formation rate to the gas density using the Schmidt law, (Kennicutt, 1998). Putting these relations all together, I find the gamma-ray flux:

 Fγ ≈ 0.0036 ergsec−1cm−2 fπ⎛⎝SFRM\sun yr−1⎞⎠(R250 pc)−2 (6) Fγ ≈ 0.0252 ergsec−1cm−2 fπ(⟨Σg⟩g cm−2)1.4 (7)

and finally,

 Fγ≈⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩0.157 ergsec−1%cm−2(⟨Σg⟩g cm−2)2.4(⟨Σg⟩≲0.16 g cm−2)0.0252 ergsec−1cm−2(⟨Σg⟩g cm−2)1.4(⟨Σg⟩≳0.16 g cm−2) (8)

Evaluating the cross section for pair production (Berestetskii et al., 1979),

 σγZ=289αFSr2e(ln2Eγmec2−10942) (9)

at to be , I find

 Unknown environment '% (10)

The results are plotted in Figure 2. These ionization rates are small, though they do rise to for the densest ULIRGs with . In principle, if is 1, these ionization rates can be ten times higher. In any case, though, they remain much higher than the ionization rate from long-lived radioactive elements such as K in the Milky Way, (Umebayashi & Nakano, 1981, 2009) even for mean column densities of .

In starbursts, there may be additional radioactive ionization from relatively short-lived isotopes like Al because of the high supernova rate (c.f., Diehl et al., 2006); for example, in the young Solar System, Al alone sustained an ionization rate of (Umebayashi & Nakano 2009; see also Stepinski 1992). In a companion paper (Lacki, 2012), I find that the mean Al ionization rate in starbursts quite high because of their rapid specific star-formation rates, of the order . If true, gamma-ray ionization only dominates in the most extreme starbursts with . However, whether these high ionization rates are actually attained in the molecular gas of starbursts depends on whether the Al is actually well-mixed with the gas; as with CRs, the propagation of Al from its sources matters. This is not an issue with the gamma rays.

Some starbursts are in fact gamma-ray detected, allowing us to directly calculate . The brightest starburst in the gamma-ray sky in terms of observed flux is M82, with a GeV-TeV gamma-ray luminosity of (Abdo et al., 2010) and a starburst radius of 250 pc (e.g., Goetz et al., 1990; Williams & Bower, 2010). Therefore, using eqn. 2, I find

 ζγ=1.1×10−19 sec−1(fion0.1)(Eion30 eV)−1, (11)

much smaller than even in Milky Way molecular clouds, but still much larger than the ionization rate from long-lived radioisotopes like K.

Arp 220, which if star-formation powered has a star-formation rate of (Torres 2004 and references therein; Anantharamaiah et al. 2000), is a more likely home for GRDRs. Most of its radio emission comes from two 50 - 100 pc radius nuclei (Norris, 1988; Downes & Solomon, 1998), although there are varying estimates of the literature of how much of the infrared luminosity is actually from the nuclei (Downes & Solomon, 1998; Downes & Eckart, 2007; Sakamoto et al., 2008). If I assume that each nucleus has a star-formation rate of (for a bolometric luminosity of each) and that Arp 220 is a proton calorimeter, then I find:

 ζγ≈4.8×10−17 sec−1(R50 % pc)−2(fion0.1)(Eion30 eV)−1. (12)

Thus, even deep molecular cores in Arp 220’s nuclei should have an ionization rate comparable to the molecular clouds of the Milky Way.

A further effect which may enhance the GRDR ionization rate is clumpiness in the gamma-ray emission of the starburst. I assumed that the gamma-ray emission was evenly distributed across the starburst. However, some molecular clouds in starbursts may be located near gamma-ray sources. The outer layers of the molecular cloud likely themselves generate gamma rays as CRs are stopped by them, enhancing the interior gamma-ray flux and gamma-ray ionization rate (e.g., Umebayashi & Nakano, 1981).

### 2.2. Comparison with CRDR Ionization Rates

How does the ionization rate in GRDRs compare with that in the surrounding CRDRs? One way of estimating the ionization rate from CRs is to scale the Milky Way molecular cloud ionization rate by the energy density of CRs. A potential problem with this approach, however, is that the CR spectrum has a different shape in starbursts than in the Milky Way. Whereas the Milky Way has a gamma-ray spectrum because CR protons escape through diffusion (e.g., Ginzburg & Ptuskin, 1976), starbursts are observed to have hard gamma-ray spectra from strong advection and/or pionic losses for protons (Acero et al., 2009; Acciari et al., 2009; Abdo et al., 2010). Therefore, some of the larger CR energy density in starbursts is in very high energy particles which are not responsible for ionization.5 Furthermore, this approach requires a CR ionization rate for the Milky Way and a Milky Way CR energy density.

Another way to estimate the CR ionization rate is to use the injected power in CR protons that is lost to ionization:

 Lion=∫dQdEKfpion(K)dE, (13)

where is the fraction of energy of a CR with kinetic energy that is lost to ionization and is the injection spectrum of CR protons. The fraction of energy going into ionization is roughly:

 fpion≈tlifetion≈[1+tpiontπ+tpiontwind]−1. (14)

Supposing the pionic losses to have a threshold of , the rest energy of a pion,

 fpion≈⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩[1+tpiontwind]−1(K<140 MeV)[1+tpiontπ1fπ]−1(K>140 MeV). (15)

The lifetime of protons to ionization losses is approximately

 tpion≈159 Myr(nHcm−3)−1(KGeV)β (16)

where is the speed of the proton divided by (Torres, 2004, and references therein). This gives us, roughly, an ionization to pionic lifetime ratio of .

Now suppose the injection spectrum is a power-law in total energy, extending from () to infinity: . The normalization of this power law is set by , giving . When and , I find that . Another commonly assumed injection spectrum is a power law in momentum. For and , with minimum kinetic energies ranging from 1 MeV to 100 MeV, I find . I therefore adopt as a fiducial value in the calculations below.

The CR ionization rate is

 ζCR=ηionLCREionV⟨nH⟩=ηionLCRmHEionMH (17)

after defining . In this equation, is the gas mass in the starburst. Scaling to values typical in Arp 220’s nuclei, I find

 ζCR=4.4×10−15 sec−1⎛⎝SFR50 M% \sun yr−1⎞⎠(MH4×108 M% \sun)−1(Eion30 eV)−1. (18)

This is times that of Milky Way molecular clouds or an Arp 220 GRDR.

A convenient expression for the ratio of the ionization rates in a GRDR to a CRDR is the ratio of equations 2 and 17:

 ζγζCR≈fπσγZfγionh⟨nH⟩3ηion≈fπσγZfγion⟨Σg⟩6ηionmH. (19)

where I have assumed a disk geometry for the volume (). In Arp 220, this ratio reaches

 ζγζCR≈0.010fπ(⟨Σg⟩10 g cm−2)(fγion0.1)(ηion0.1)−1. (20)

Note this ratio increases even once proton calorimetry is attained. This is because at high densities, the pionic lifetime continues to decrease as , so that the energy density of CRs per supernova decreases; the energy density of gamma rays per supernova instead remains constant. Only when for the entire starburst does the gamma-ray flux also start to diminish.

### 2.3. The Ionization Fraction of GRDRs: Consequences for Ambipolar Diffusion

McKee (1989) gives the ionization fraction of a cloud with density as

 xe=1×10−7r−1gd(nchnH)1/2[(1+nch4nH)1/2+(nch4nH)1/2]. (21)

Here, is a characteristic density, and is the gas-to-dust ratio divided by 100 (see also Papadopoulos 2010). For most GRDRs, the ionization rate is low enough that and this can be approximated as

 xGRDRe=3.2×10−9(ζγ10−17 sec−1)1/2(nH106 cm−3)−1/2 (22)

The ionization fraction when is plotted in Figure 3. In Arp 220, this fraction rises to for . However, in using equation 21, I have ignored the possible appearence of different chemical effects in the very low ionization environments of GRDRs, such as charge being carried by dust grains (McKee, 1989).

The ionization fraction is important for determining the ambipolar diffusion time, the time for the magnetic field to slip away from the neutral gas. If the magnetic field remains strong in a cloud, free-falling gravitational collapse cannot proceed (e.g., Mouschovias & Spitzer, 1976). When the ambipolar diffusion time is longer than the free-fall time and the magnetic fields are strong, the collapse instead proceeds quasi-statically on the ambipolar diffusion time (e.g., Mestel & Spitzer 1956; see also the discussion in Crutcher 1999). In most Milky Way starless cores, the ambipolar diffusion timescale is roughly an order of magnitude longer than the free-fall time (e.g., Caselli et al., 1998; Williams et al., 1998; Maret & Bergin, 2007; Hezareh et al., 2008). There is some recent evidence that at least some of the very deepest Milky Way cores may have low ionization fractions and ambipolar diffusion times approximately equal to the free-fall time (e.g., Caselli et al., 2002; Bergin & Tafalla, 2007). Given that strong magnetic fields are observed by Zeeman splitting in dense OH maser regions in ULIRGs (Robishaw et al., 2008), it seems plausible that magnetic fields regulate star formation in starbursts, so understanding the ambipolar diffusion rate is important.

McKee (1989) also gives the ambipolar diffusion time as . In the low ionization environments of GRDRs with , we have

We can compare this with the local free-fall time :

 tff=5.1×104 yr(nH106 cm−3)−1/2. (24)

The interesting point for star-formation is that both the ambipolar diffusion time and free-fall time scale with . Thus all GRDRs in a given starburst galaxy will have the same ratio of ambipolar and free-fall times, regardless of their overdensity. This ratio is

The ratio is also plotted in Figure 3; from the Schmidt Law, we expect ambipolar diffusion in GRDRs to be slower than free-fall in galaxies with . Even in M82, this ratio is , so that ambipolar diffusion is no quicker than free-fall time. But in a starburst like Arp 220, this ratio reaches .

The ratio of ambipolar diffusion timescale to free-fall time in dense starburst GRDRs is therefore similar to those observed in Milky Way cores, . This means that star-formation may proceed similarly in dense starbursts and normal galaxies. By contrast in CRDRs, the predicted ambipolar diffusion time is vastly longer than the free-fall time; star-formation would effectively be halted unless some mechanism, perhaps turbulence, speeds up ambipolar diffusion (Papadopoulos, 2010). It is interesting that in Arp 220, approaches , the empirically observed timescale for star-formation (e.g., Kennicutt, 1998; Krumholz et al., 2012). This suggests that in Arp 220, ambipolar diffusion in GRDRs may alter the star-formation efficiency slightly.

### 2.4. The Temperature of GRDRs

Papadopoulos (2010) gives a convenient expression for the minimum gas kinetic temperature of a CRDR which also applies for GRDRs, assuming no dust or turbulent heating:

 Tmink=6.3 K [(0.707n1/24ζ−17+0.1862n34)1/2−0.186n3/24]2/3 (26)

where and . The temperatures for GRDRs and CRDRs are plotted in Figure 4.

Simpler expressions can be found for in the limit of high density and low ionization rate, where a Taylor series expansion yields:

 Tmink=9.7 K n−2/34ζ2/3−17. (27)

Thus, using the expressions for above, I derive approximate gas kinetic temperatures:

 Unknown environment '% (28)

The minimum temperatures of GRDRs are clearly very low: only if (and ) is very large and the local density is relatively small () do the temperatures exceed 1 Kelvin. For the values of derived for M82 in section 2.1, I calculate to be 0.5 K if and 0.02 K if . For Arp 220’s nuclear starbursts, I find is 9 K if and 1.3 K if . For comparison, the CRDR temperature in Arp 220 is for (using eqn. 18 for ). It seems that in most starbursts GRDRs would be much colder than molecular clouds in the Milky Way.

Of course, the temperatures derived in equation 26 are minimum temperatures derived by assuming dust has a temperature of 0, so that the gas is cooled by collisions with dust grains. In reality, clouds will be embedded in an intense far-infrared radiation field (and the CMB), heating its dust grains. If the dust temperature is not cooled significantly by dust-grain interactions, if there is no source of dust heating – whether light from young stars or absorption of molecular line emission – within the GRDR, and if the dust is in thermal equilibrium with the external radiation field, there will be no net flux in dust thermal emission through the cloud because of energy conservation, so that the dust temperature should be equal to the external effective temperature.

For a given dust temperature, dust-grain interactions will set a gas kinetic temperature floor where the dust-grain heating rate equals the molecular line cooling rate . Papadopoulos (2010) quotes

 Γg−d=2.5×10−26(nH104 cm−3)2(TkK)1/2(Tdust−TkK) erg cm−2 sec−1 (29)

and

 Λline=4.2×10−24(nH104 cm−3)1/2(Tk10 K)3 erg cm−2 sec−1. (30)

Equating these two rates, I find a temperature floor of

 Tfloor=8.9 K(nH104 cm−3)3/5(Tdust40 K)2/5 (31)

when the density is low and , and

 Tfloor=Tdust−1.7 K(nH106 % cm−3)−3/2(Tdust40 K)5/2 (32)

when the density is high and . GRDRs with will have a kinetic temperature near the dust temperature when . At lower densities, if the heating and cooling rates are accurate, the gas will be much colder, though given the extreme conditions, the low temperatures should be approached with caution.

Another possible source of heating is the dissipation of turbulence. Most of the molecular gas in starbursts exists in a highly turbulent phase (e.g., Downes & Solomon, 1998). Supersonic turbulence rapidly dissipates, partly because the supersonic motions generates shocks (e.g., Stone et al., 1998). The cores in molecular clouds and protostars achieve densities much higher than most of the mass in the turbulent ISM – they are said to be “decoupled” from the surrounding medium – although such collapsed structures can be initially generated by turbulence (e.g., Klessen et al., 2000). It is these clouds which have the high columns and densities where GRDRs are most likely to be present. Bergin & Tafalla (2007) argues that in these clouds, turbulence is unimportant, based on the observed relationship between turbulent line width and cloud size in Milky Way molecular clouds (Myers, 1983; Larson, 2005).

The cold temperatures of GRDRs have implications for their Jeans mass, which may set the minimum mass of stars in starbursts and affect their initial mass function (IMF). The Jeans mass is

 MJ=(kBTkGμmH)3/2(nHmH)−1/2. (33)

where is the mean particle mass in the molecular ISM, in terms of the atomic mass of hydrogen (Papadopoulos et al., 2011). In dust-heated gas with temperatures near , I find the Jeans mass is

 MJ ≈ 1.2 M\sun(Tk9 K)3/2(nH104 cm−3)−1/2 (34) ≈ 1.1 M\sun(Tk40 K)3/2(nH106 cm−3)−1/2. (35)

For a dust temperature of , the Jeans mass reaches a maximum of when . Thus, the Jeans mass of GRDRs is similar or perhaps slightly higher than that in normal galaxies, unlike in CRDRs, where it may be much higher () in starbursts (c.f. Papadopoulos et al., 2011).

I conclude that the temperature of GRDRs is more likely set by the dust grains being illuminated by infrared radiation or possibly turbulence than gamma-ray heating, except perhaps in the lower density regions of the most extreme starbursts.

## 3. The Transition from Cosmic Rays to Gamma Rays: The Perils of Cosmic Ray Scattering

At first glance, the ionization rates and equilibrium temperatures of GRDRs seem surprisingly low compared to predictions for CRDRs, given the high gamma-ray luminosities of starbursts. In a proton calorimetric starburst, is about one-third of , so why doesn’t equation 1 imply that the CRDR heating rate and ionization rate are also low? The reason is that galaxies are extreme scattering atmospheres for CRs, at least for clouds of sufficiently low CR absorption optical depth. The flux in equation 1 is not the net flux, which is relatively small in galaxies, but the total flux, which is much larger. Because of the large scattering optical depth of a galaxy, a CR can cross a cloud multiple times: this increases the odds that it will ionize an atom in it, and therefore increases the CR ionization and heating rate. However, the disadvantage of this is that a CR has a much greater chance of being destroyed by pion production or ionization cooling, precisely because it can cross the clouds multiple times.

The slower CRs traverse a molecular cloud, the quicker that the CR flux is attenuated (assuming no CR accelerators within the cloud), and the more likely that gamma rays provide the ionization. At best, however, CRs can traverse the cloud at , with no deflection by magnetic fields within the cloud. I consider the 1D version of this case in the Appendix, where a molecular cloud with an absorption optical depth to CRs is embedded in a medium with CR absorptivity and scattering coefficient . Far from the cloud, the CR ionization rate is . In general, : galaxies are scattering atmospheres for CRs. I find that even with free travel within a molecular cloud, the CR flux is attenuated once the pionic absorption optical depth of the cloud equals (much less than 1). To take a specific example, the radio nuclei of Arp 220 have an average density of , so the pionic lifetime of CRs within them is according to eqn. 4. If the (unknown) scattering mean free path within them is , then a cloud effectively becomes optically thick to pion production in this model when

 Σabs≈√λdiff/(ctπ)Σπ≈0.6 g cm−2. (36)

The rapid attenuation of CRs is shown in Figure 5. The short effective attenuation length arises simply because CRs cross the cloud many times after being scattered by the surrounding galaxy.

On the other hand, it takes several optical depths before the CR ionization rate drops below the gamma-ray ionization rate. Equation 14 is an expression for the optical depth of a cloud that attenuates the CR flux enough so that gamma-ray ionization dominates in the center (). These values are plotted in Figure 6: I find CR absorption optical depths of are sufficient for gamma-ray ionization to become more important than CR ionization. Plugging in values for Arp 220’s nuclei from eqn. 20 (), I find gamma-ray ionization starts to dominate once . For a typical pionic absorbing column of , this comes out to . A cloud with this column and density will have a mass . Turbulence in the molecular ISM does readily generate transient density fluctuations, but simulations show that typical column densities should still be near (Ostriker et al., 2001), which is generally less than . However, turbulence can also generate Jeans-unstable, collapsing structures such as protostars, and these can easily achieve such columns (Klessen et al., 2000). Protostars with are very likely to be GRDRs in Arp 220’s nuclei, since at these densities, even if they are not shielded by a surrounding molecular cloud.

Moreover, the necessary for a GRDR grows only logarithmically with . When CR scattering in the surrounding medium is slow enough (), the necessary optical depth is approximately

 τGRDR≈−ln[ζγ/ζextCR2√αout/σout]. (37)

The larger CR lifetimes of weaker starbursts also means that CRs cross a cloud more times before being destroyed outside of the cloud: the probability of the CR being destroyed inside of the cloud rather than outside goes up. For M82 and NGC 253, with and a typical CR lifetime of , I find , or , when . Of course, gamma rays themselves start to be attenuated over such columns, but as will be seen in the next section, cascade emission from the pair can ensure the gamma-ray flux is preserved, to order of magnitude, out to these columns.

Of course, since these calculations are 1D, they are only illustrative. In reality, the 1D calculations apply only for clouds that are essentially large sheets. In three dimensions, CRs have more directions to go that do not intersect clouds. Therefore, CRs are more likely to miss small, spherical clouds, and would cross it fewer times; the CR flux reduction need not be so drastic.

However, the above calculations are conservative in assuming that CRs can freely enter into a cloud and cross it at the speed of light. If instead diffusion is slow inside a molecular cloud, the CR flux within it will be much smaller. The column through which CRs can penetrate is then . The CRs can only diffuse into clouds with a column of less than

 Σ≈0.038 g cm−2(D1028 cm2 sec−1)1/2(n104 cm−3)1/2. (38)

Thus, if CRs diffuse slowly in molecular clouds, the GRDRs may encompass large fractions of the molecular material within a starburst. Note that hard X-rays are also not attenuated under columns , so they would contribute to ionization at such low column densities. The increase of the column with density arises because the effective speed at which CRs diffuse increases at smaller scales. Below the mean free path, the diffusion approximation breaks down, and a treatment similar to the one above, where the CRs are assumed to free-stream, must be used.

The actual propagation of CRs in molecular clouds is poorly understood. It is possible that CRs essentially free-stream through molecular clouds (Cesarsky & Volk, 1978; Morfill, 1982), which could occur if the plasma waves that normally scatter CRs are quickly damped by collisions (Kulsrud & Cesarsky, 1971). However, there are reasons to expect slow diffusion in or near molecular clouds. Skilling & Strong (1976) argued that because molecular clouds can absorb CRs, there is a net flux of CRs into the cloud; the net flux causes plasma instabilities in the ionized material surrounding the cloud, generating waves that scatter ionizing CRs away from the cloud. More simply, CRs may be excluded from molecular clouds simply because the magnetic field is higher, which might lead to a smaller diffusion constant (e.g., Gabici et al., 2007). Observationally, there is radio evidence that diffusion is slow inside the Galactic Center molecular cloud Sgr B (Protheroe et al., 2008; Jones et al., 2011). The diffusion constant also appears to be small in the vicinity of the molecular clouds surrounding the supernova remnants W28 and IC 443, although regions outside the molecular clouds are included in these measurements (e.g., Fujita et al., 2009; Gabici et al., 2010; Torres et al., 2010); the small diffusion constant could be a result of plasma waves generated by the CR flux from the nearby supernova remnants in the ionized ISM (Fujita et al., 2010). Crocker et al. (2011) interpret the relative gamma-ray dimness of the Galactic Center as being caused by the exclusion of CRs from molecular material, although M82 and NGC 253 seem instead to be gamma-ray bright.

## 4. The Development of the Pair Cascade: The Long Reach of Gamma Rays

Gamma rays do not directly ionize gas; instead they pair produce that ionize the gas. The largest uncertainty in the estimates of the GRDR ionization rate is the efficiency that the energy in these pair is converted into ionization of the ISM. Pair undergo a cascade, emitting high energy radiation, some of which may themselves pair produce . At low energies, ionization losses (with an energy-independent ) always dominate over the other losses, so low energy gamma rays are converted into ionization power relatively efficiently. Above a few hundred MeV, bremsstrahlung losses dominate over ionization losses at all densities, since for bremsstrahlung. Bremsstrahlung radiation from an electron of energy produces gamma rays with a typical energy (Schlickeiser, 2002). In optically thin clouds, the gamma rays escape, so that bremsstrahlung losses waste the pair energy that could go into ionization. In clouds of high enough column, though, the bremsstrahlung gamma rays can themselves pair produce and have another chance to ionize the ISM. At high energy, synchrotron and Inverse Compton (IC) losses with are the main cooling mechanism for pair . The synchrotron radiation from of relevant energies is in the form of radio and infrared emission, which either escapes or is absorbed by dust grains. Synchrotron emission therefore does not contribute to the ionization of GRDRs. IC radiation is ultraviolet light, X-rays, and gamma rays. Ultraviolet and X-rays can directly ionize atoms. As with bremsstrahlung, gamma rays either escape or pair produce depending on the cloud column density.

Estimating requires modeling of all of these cooling processes. Furthermore, the gamma-ray spectral features of starburst galaxies must be taken into account. The gamma rays of starbursts are thought to be mostly pionic; because of the kinematic threshold of pion production, the pionic spectrum should drop away at energies below 70 MeV in (e.g., Stecker, 1970), or a few hundred MeV in (for example, for the Galactic pionic spectrum, as in Kamae et al. 2006). At high energies, the spectrum will be a hard power law. As a result, most of the pair will be at high energy, where they will not directly ionize the gas but instead cool by radiation. Estimating the ionizing effects of the resulting gamma-ray cascade is the subject of this section.

I make many simplifying assumptions. Once again, I use a 1D two stream model for the gamma rays. In this case, I assume that neither diffuse nor are advected from the regions they are created (an on-the-spot approximation), which would be the case if in the in the dense environments of molecular clouds.6 I can phrase the equation of radiative transfer in terms of the column instead of the physical depth . Each generation of gamma rays produces pair , which in turn generate the gamma rays of the next generation . The equation of radiative transfer for these assumptions is:

 dIi±dNH=±(−σγZIi±+Qi−1brems+Qi−1IC2nH), (39)

where for the primary gamma rays (), , and is the total cross section for pair production. The intensity is a number intensity per unit energy here rather than energy intensity. The 1/2 factor for the terms takes into account that half of the emitted radiation would be emitted towards higher columns (in ), and half towards lower columns (in ), since the direction of are likely to be scrambled by magnetic fields.

Pair production and the radiation from the pairs is calculated using -function approximations. For pair production, the pair from a photon of energy are assumed to all be of energy :

 Qie=4αγZ(2Ee)∫Ii(2Ee)dΩ=4(Ii+(2Ee)+Ii−(2Ee))αγZ(2Ee). (40)

I then estimate the pair spectrum as

 Nie=Qie×tlife (41)

where includes losses from ionization, bremsstrahlung, synchrotron, and Inverse Compton. Defining the loss times for each process as , I use the Bethe-Bloch energy losses for ionization from Strong & Moskalenko (1998), the bremsstrahlung losses in a neutral hydrogen gas from Strong & Moskalenko (1998), and the standard Thomson energy loss formulas for synchrotron and IC emission from Rybicki & Lightman (1979).

Bremsstrahlung emission is assumed to come out at energies :

 Qi+1brems=4Nie(2Eγ)tbrems. (42)

Finally IC radiation is assumed to come out at an energy (Rybicki & Lightman, 1979), where is the energy of photons in the background radiation field. Since IC spreads the energy of one dex in energy into two dex of upscattered photon energy, I have :

 Qi+1IC=E2eNie(Ee)2E2γtIC. (43)

The volumetric energy input in ionization losses is calculated as

 Γion=∫EeNietiondEe. (44)

In addition, I assume that Inverse Compton photons with energies between 13.6 eV and 1 MeV are immediately absorbed, and contribute all of their energy to ionization. EUV and soft X-ray () photons can directly ionize neutral atoms through photoabsorption, though in practice, they can also heat the ISM by exciting atoms without ionizing it. Hard X-rays lose energy through Compton scattering in large column depths (). The energy lost per scattering is hundreds of eV to many keV, so these scatterings eject electrons from atoms, which then can act as secondary electrons to ionize other atoms. I calculate the volumetric low energy IC energy input as

 ΓLE−IC=∫1 MeV13.6 eVEγQICdEγ. (45)

Then the ionization rate is just

 ξ=ΓLE−IC+ΓionnHEion. (46)

I applied this method to an input gamma-ray spectrum from the fiducial models of Lacki & Thompson (2010) of M82, which is detected in gamma rays, and Arp 220’s east nucleus. The western nucleus of Arp 220 would likely provide roughly similar results (c.f. Torres, 2004), but it is unclear how much it is powered by an active galactic nucleus (e.g., Downes & Eckart, 2007). The input gamma-ray spectrum is in terms of unabsorbed volumetric power emitted ; I convert it to an intensity by multiplying by a midplane-to-edge scale height of , set to 100 pc for M82 and 50 pc for Arp 220’s east nucleus: . I assume the modeled column is completely opaque, so that . The spectrum extends down to 1 MeV, and I cut it off at 10 TeV, since higher energy gamma rays are likely absorbed by pair production processes (Torres, 2004; Inoue, 2011; Lacki & Thompson, 2010). I adopt a density of and a magnetic field strength of 10 milliGauss, similar to those found by Zeeman splitting measurements of dense OH masers (Robishaw et al., 2008). The radiation field has a blackbody spectrum with temperature of 40 K for M82 and 50 K for Arp 220, but is scaled to a characteristic radiation energy density . M82 is optically thin in the far-infrared, so I calculate , where is the total IR luminosity from Sanders et al. (2003) and is the radius. Arp 220’s nuclei are optically thick in the far-infrared (e.g., Downes & Eckart, 2007; Sakamoto et al., 2008; Papadopoulos et al., 2010), so I just use the blackbody energy density.

The resulting ionization rates for a large column of material are displayed in Figure 7. I find that, to order of magnitude, the gamma-ray ionization rate at low column depths is similar to that derived in section 2.1 using . In M82, , about three times higher than predicted; in Arp 220 East, is about , about twice the estimated value. In reality, these values will vary for GRDRs with different densities, magnetic field strengths, and radiation fields.

In both M82 and Arp 220 East, the ionization rate is maintained to enormous column densities. In the M82 GRDR, the ionization rate remains above to (), above through (), and above the K radioactivity-induced ionization rate of