{'fraction_non_alphanumeric': 0.03945045667357434, 'fraction_numerical': 0.023217438022872055, 'mean_word_length': 4.369668246445498, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simulating Gamma-ray Emission in Star-forming Galaxies

Star forming galaxies emit GeV- and TeV-gamma rays that are thought to originate from hadronic interactions of cosmic-ray (CR) nuclei with the interstellar medium. To understand the emission, we have used the moving mesh code Arepo to perform magneto-hydrodynamical galaxy formation simulations with self-consistent CR physics. Our galaxy models exhibit a first burst of star formation that injects CRs at supernovae. Once CRs have sufficiently accumulated in our Milky-Way like galaxy, their buoyancy force overcomes the magnetic tension of the toroidal disk field. As field lines open up, they enable anisotropically diffusing CRs to escape into the halo and to accelerate a bubble-like, CR-dominated outflow. However, these bubbles are invisible in our simulated gamma-ray maps of hadronic pion-decay and secondary inverse-Compton emission because of low gas density in the outflows. By adopting a phenomenological relation between star formation rate (SFR) and far-infrared emission and assuming that gamma rays mainly originate from decaying pions, our simulated galaxies can reproduce the observed tight relation between far-infrared and gamma-ray emission, independent of whether we account for anisotropic CR diffusion. This demonstrates that uncertainties in modeling active CR transport processes only play a minor role in predicting gamma-ray emission from galaxies. We find that in starbursts, most of the CR energy is"calorimetrically"lost to hadronic interactions. In contrast, the gamma-ray emission deviates from this calorimetric property at low SFRs due to adiabatic losses, which cannot be identified in traditional one-zone models.


INTRODUCTION
CRs, magnetic fields and turbulent motions contribute equally to the total midplane pressure in the Milky Way (Boulares & Cox 1990). This could originate from a selfregulated feedback loop and may even suggest that CRs play an active role in shaping galaxies by altering the dynamics and driving galactic winds through their gradient pressure force, as suggested by theoretical works (Breitschwerdt et al. 1991;Everett et al. 2008) threedimensional simulations of galaxies (Uhlig et al. 2012;Hanasz et al. 2013;Booth et al. 2013;Salem & Bryan 2014;Pakmor et al. 2016b) and the interstellar medium (ISM) (Girichidis et al. 2016;Simpson et al. 2016). This idea can be scrutinized by studying CR-induced radiative processes with the goal to extract the CR pressure from the non-thermal emission at radio and gamma-ray energies.
Most of the galactic diffuse emission is generated by massive stars; preferentially in the beginning and at the end of their lives. Young massive stars emit mostly ultraviolet (UV) photons, which are absorbed by dust grains. This radiation is then re-emitted in the far infrared (FIR) provided the dust is optically thick to UV photons. This is the case in actively star-forming galaxies, supporting the phenomenological correlation of the FIR emission with the SFR (Kennicutt 1998a). After their fuel is exhausted, massive stars explode as core-collapse supernovae, whose remnants are believed to accelerate CR protons and electrons in galaxies through diffusive shock cpfrommer@aip.de acceleration.
CR electrons gyrate in the interstellar magnetic field and emit radio synchrotron radiation that closely correlates with the total FIR luminosity of galaxies over five orders of magnitude in luminosity (Bell 2003). The same radio-emitting CR electrons can generate high-energy gamma-ray emission either through free-free transitions in the vicinity of gas nuclei (bremsstrahlung) or by inverse Compton (IC) scattering off of radiation fields. Hadronic collisions between CR nuclei and the ISM also produce gamma rays alongside other decay products via pion decay: π 0 → 2γ , π ± → µ ± + ν µ /ν µ → e ± + ν e /ν e + ν µ +ν µ .
These secondary relativistic electrons and positrons (hereafter simply electrons) can emit radio-synchrotron emission in the presence of ubiquitous galactic magnetic fields as well as up-scatter ambient photon fields into the gamma-ray regime through IC interactions.
Recently, nearby starburst galaxies have been detected in high-energy gamma rays by space-based (Abdo et al. 2010) and imaging air-Cerenkov telescopes (VERITAS Collaboration et al. 2009;Acero et al. 2009). Upper limits and detections of the galactic gamma-ray emission by the Fermi telescope revealed a tight FIR-gamma-ray relation (Ackermann et al. 2012;Rojas-Bravo & Araya 2016). The gamma-ray emission enables testing whether starbursts act as proton "calorimeters" (where most of the CR energy is lost to hadronic interactions, Lacki et al. Here, we present high-resolution MHD simulations with CR physics of forming disk galaxies and compute the gamma-ray emission. This enables us to critically assess simplifying assumptions made in previous one-zone analyses. We describe our simulations in §2 and lay out our methodology to compute pion-decay and IC emission in §3. We analyze gamma-ray emission maps and the FIR-gamma-ray relation in §4 and conclude in §5.

SIMULATIONS
We simulate the formation of isolated disks in darkmatter (DM) halos that range in mass from 10 10 to 10 12 M . We use the second-order accurate, adaptive moving-mesh code Arepo (Springel 2010;Pakmor et al. 2016c). We model the ISM with an effective equation of state with radiative cooling and star formation using a probabilistic approach (Springel & Hernquist 2003). Magnetic fields are treated with ideal MHD (Pakmor & Springel 2013).
CRs are modelled as a relativistic fluid with a constant adiabatic index of 4/3 in a two-fluid approximation . We describe CR generation at remnants of core-collapse supernovae by instantaneously injecting all CR energy produced by a stellar particles population through this channel with an energy efficiency of ζ SN = 0.1 into its surroundings immediately after birth. We account for adiabatic changes of the CR energy as well as CR cooling via Coulomb and hadronic CR interactions, assuming an equilibrium momentum distribution . To bracket the uncertainty regarding CR transport, we simulate two models. In model CR adv, we only advect CRs with the gas while model CR diff additionally accounts for anisotropic diffusion relative to the rest frame of the gas with a diffusion coefficient of 10 28 cm 2 s −1 along the magnetic field and no diffusion perpendicular to it (Pakmor et al. 2016a). Our simulations do not include CR streaming (unlike Wiener et al. 2017).
DM halos are modelled as static NFW (Navarro et al. 1997) profiles with a fixed concentration parameter of c 200 = 12 across our halo mass range and one 10 12 M halo with c 200 = 7. The gas is initially in hydrostatic equilibrium with the DM potential. The halo carries a small amount of angular momentum, parametrized by the spin parameter λ = 0.05. In all cases, we adopt a baryon mass fraction of Ω b /Ω m = 0.155. We assume a uniform, homogeneous seed field along the x-axis with strength 10 −12 G and no CRs in the initial conditions. We start our simulations with 10 7 gas cells inside the virial radius, each of which has a target mass of 155 M × M 10 , where M 10 = M 200 /(10 10 M ). This target gas mass corresponds to the typical mass of a stellar population particle. We require that the mass of all cells is within a factor of two of the target mass by explicitly refining and de-refining cells. We additionally require adjacent cells to differ in volume by less than a factor of 10 and refine the larger cell if this condition is violated.

RADIATIVE PROCESSES
3.1. Pion-decay gamma rays Pion-decay gamma rays arise from inelastic CR interactions with thermal nuclei. We assume that the one-dimensional CR momentum spectrum per volume element, d 3 x, follows a power-law: where p = P p /(m p c) is the dimensionless proton momentum, θ(a) denotes the Heaviside step function, q is the low-momentum cutoff and α is the CR spectral index. The normalization is determined by solving the CR energy density (ε cr ) for C p : where E(p) = ( 1 + p 2 − 1) m p c 2 is the kinetic proton energy and B y (a, b) denotes the incomplete beta function (Abramowitz & Stegun 1965), assuming α > 2. The omnidirectional gamma-ray source function from decaying pions (s π 0 −γ ) is (Pfrommer & Enßlin 2004): where m π 0 is the pion rest mass and n n = n H + 4n He = ρ/m p is the target density of nucleons, neglecting metals. This treatment accounts for all physical processes at the pion production threshold, which are parametrized by the shape parameter δ and the effective cross section σ pp in terms of the photon index α γ = α according to δ 0.14 α −1.6 + 0.44 and (4) σ pp 32 × 0.96 + e 4.4 − 2.4 α mbarn.
The energy-integrated gamma-ray emissivity for the energy band [E 1 , E 2 ] (in units of erg s −1 cm −3 ) that results from pion decay is given by 3.2. Secondary inverse Compton emission The mean energies of isotropically scattered IC photons and scattering CR electrons are related by where the particle kinetic energy E/(m e c 2 ) = γ −1 is defined in terms of the Lorentz factor γ. We adopt cosmic microwave background (CMB) photons with a characteristic energy hν CMB 0.66 meV as the source for IC emission using Wien's displacement law. The IC cooling time of these relativistic electrons is t IC ∼ 2 Myr and thus shorter than both the CR residency time in our Galaxy and the dynamical time scales in the warm and coronal phases of the ISM. Hence, the steady-state approximation for modeling secondary IC emission is justified. Note that we neglect IC interactions with starlight and FIR photons: including them would not change any conclusions since the steady state IC emissivity does not depend on the photon energy density (in the IC-dominated scattering regime), which solely depends on the CR electron energy density. Throughout the Letter, we neglect primary IC emission, which has a subdominant contribution to the Fermi-band luminosity (but may be important outside the disk, Selig et al. 2015).
At high momenta (P e > GeV/c), the injection of secondaries is balanced by IC and synchrotron cooling, which results in an equilibrium distribution of secondary CR electrons (Pfrommer et al. 2008), f e (p) dp = C e p −αe dp (9) where we redefined p = P e /(m e c) as the dimensionless electron momentum, α e = α + 1 is the electron spectral index, σ T is the Thomson cross section, and ε B and ε ph are the energy densities of the magnetic and photon fields, respectively. The energy-integrated gamma-ray emissivity for an isotropic power law distribution of CR electrons of equation (9) that Compton-upscatters CMB photons is (derived from eq. (7.31) in Rybicki & Lightman (1979), in the case of Thomson scattering), f IC (α e ) = 2 αe+3 (α 2 e + 4 α e + 11) (α e + 3) 2 (α e + 5) (α e + 1) where α ν = (α e − 1)/2 denotes the IC spectral index, T CMB = 2.725 K is the CMB temperature at the present time, r e = e 2 /(m e c 2 ) the classical electron radius, e is the elementary charge, h is Planck's constant and Γ(a) and ζ(a) are the gamma and Riemann ζ functions, respectively (Abramowitz & Stegun 1965). The gammaray luminosities follow as a result of volume integrations of the emissivities,

CR-driven outflows in Milky Way-like galaxy
We first analyze our simulation of the Milky Waylike galaxy using the model CR diff (M 200 = 10 12 M , c 200 = 7). At the beginning of our simulation, radiative cooling diminishes pressure support of the central dense gas, starting collapse while conserving the gas' specific angular momentum. After settling into a rotationally supported disk, gas is compressed by self-gravity to sufficiently high densities for star formation. During the first collapse phase, a turbulent dynamo quickly grows a small-scale magnetic field, which is further amplified by a large-scale dynamo that preferentially grows a toroidal field in the disk (Pakmor et al. 2016b. CRs are injected into the ambient ISM surrounding stellar macro-particles, providing the gas with additional non-thermal pressure. As CRs are advected and diffused, they collect in the disk and start to dominate the pressure, especially in the outer midplane of the disk and everywhere at the disk-halo interface (Fig. 1). This result points to the importance of modeling CRs in simulations of galaxy formation that aim to capture their dynamical evolution. While our SFRs significantly vary over the simulation, our galaxy properties are not changing too quickly on the timescale over which the CRs lose energy and hence can be considered a reasonable proxy for observed systems.
However, the CR pressure falls short of the thermal pressure in the inner part of the disk. This is caused by over-cooling of our steady-state CR fluid in the densest regions since the simulations assume a softer CR spectral index (α = 2.2) than the analysis (α = 2.05). This will be improved upon with future simulations that dynamically follow the CR spectrum. Moreover, we adopt a single-phase ISM rather than a realistic multi-phase ISM. The latter would minimize CR losses since CRs spend more time in the lower-density warm and coronal phases, which dominate the volume.
After 1 Gyr, the CR pressure has increased to the point where the buoyancy force overcomes the magnetic tension of the toroidal magnetic field which bends and opens. CRs diffuse ahead of the gas into the halo and accelerate the gas, thereby driving a bubble-like outflow (Fig. 1). The total pressure, composed of thermal, CR, and magnetic pressures, declines smoothly outwards, which demonstrates that the bubble edges are contact discontinuities and not shocks.
We project the pion-decay and secondary IC emissivities along the line-of-sight into face-on and edge-on views ( Fig. 1) according to S i−γ (n) = Λ i−γ (r,n)dr where i ∈ {π 0 , IC} andn is a unit vector perpendicular to the line-of-sight. Despite the high CR load of the outflow, neither pion-decay nor IC gamma-ray emission show signatures of hadronic gamma-ray bubbles because the fast outflow (υ z 150 km s −1 ) evacuates the gas in the bubble region. This is consistent with findings by Bayesian non-parametric reconstructions of the Fermi sky which suggests that the Fermi bubbles are of leptonic origin (Selig et al. 2015). Figure 2 compares the pion decay-induced gamma-ray spectra to the IC spectra from an equilibrium secondary CR electron population in the weak-field regime (B = 0). This choice maximizes the IC yield of the cooling electron population and keeps the pion decay-to-IC photon ratio at a constant value that solely depends on the spectral index. Increasing the field strength above the CMB equivalent field strength of B CMB = 3.24 µG, causes the electrons to cool preferentially by emitting radio syn- The MHD simulations account for CR injection at supernova remnants and follow their advection with the gas and anisotropic diffusion along the magnetic fields relative to the gas. In the top panels, we show cross-sections of gas properties in the mid-plane of the disk (face-on views) and vertical cut-planes through the centre (edge-on views) of the gas density (left), CR energy density (middle), and CR-to-thermal pressure ratio (right). In the bottom panels, we show face-on and edge-on projections of observables: gas surface mass density (left), pion decay-induced gamma-ray surface brightness from 0.1 to 100 GeV (middle) and secondary IC surface brightness in the same energy band (right). Both gamma-ray maps do not show the CR-loaded outflow in the form of hadronic gamma-ray bubbles.  -Spectral distribution of the differential gamma-ray source function (top) and the energy-integrated gamma-ray emissivity, Λγ , vs. energy (middle) that result from hadronic CR interactions for different CR spectral indices, α (colored differently). We also present Λγ vs. spectral index in the Fermi-energy band (bottom) for various low-momentum cutoffs q of the distribution function (colored differently). Shown are pion decay induced gammaray spectra (solid) and IC spectra from a steady-state secondary CR electron population (dash dotted) in the weak-field regime (B = 0). The model calculations assume a fixed CR energy density of εcr = 10 −11 erg cm −3 and target nucleon density nn = 1 cm −3 , which both impact the absolute normalization of the spectra in equal measures. chrotron radiation. Less energy is then emitted via IC interactions. This happens in our Milky Way-like galaxy for radii r 3 kpc at t = 1 Gyr and the field strength increases with time to tens of µG in the center due to the ongoing magnetic dynamo. At larger radii in the low-B regime, Λ π 0 −γ /Λ IC−γ = 4.3 for α = 2.05. Due to the strong IC suppression in the center, the pion-decay luminosity (3.2 × 10 40 erg s −1 ) dominates over the secondary IC luminosity (5.6 × 10 39 erg s −1 ) by a factor of 5.7. Hence, we neglect the IC contribution to the gammaray luminosity in what follows.
Throughout this Letter, we adopt α = 2.05 and q = 0.5 which yields an average value of Λ π 0 −γ that is within a factor of 2.5 from the extreme outliers for a broad variety of physically motivated values for α and q (see bottom panel of Fig. 2). The observed photon index is expected to be somewhat steeper (matching the observed range of α γ = 2.1 to 2.4, Ackermann et al. 2012) due to a combination of energy-dependent CR streaming and diffusion in Kolmogorov turbulence.

FIR-gamma-ray relation
We now analyze our entire galaxy sample for both models CR adv and CR diff while fixing c 200 = 12 to ensure self-similarly evolving galaxies. We believe that each of our galaxy simulations are good analogues of observed galaxies as they go through different evolutionary phases: initial gas collapse is immediately followed by a starburst that transitions to an intense star-forming thick disk, which eventually settles to a quiescently starforming thin-disk galaxy.
After 100 Myr, the 10 12 M halo enters a starburst phase with a central density of 3 M pc −3 (averaged within a radius of 300 pc). Line-of-sight integration over this central region yields a central surface mass density of 1.8 × 10 3 M pc −2 or 0.4 g cm −2 , amounting to 1.6 (2.6) times the surface mass density inferred from M 82 (NGC 253), see Lacki et al. (2011). The corresponding simulated SFR density is ≈ 40 M yr −1 kpc −3 , which integrates to a global SFR of ≈ 120 M yr −1 , substantially larger than the SFRs in M 82 and NGC 253. The simulations reach a peak resolution (minimum cell radius) of ≈ 5 pc, sufficient to resolve the central starburst region.
This assumes that thermal dust emission is a calorimetric measure of the radiation of young stars, and the factor = 0.79 derives from the Chabrier (2003) initial mass function (IMF, Crain et al. 2010). While this conversion is reliable at L 8−1000 µm > 10 9 L , it becomes progressively worse at smaller FIR luminosities due to the lower metallicity and dust content, which implies a low optical depth to IR photons and invalidates the calorimetric assumption (Bell 2003). Blindly applying the conversion yields the grey data points in Fig. 3 Fig. 3 we correlate the gamma-ray luminosity in the Fermi band (0.1 − 100 GeV) to the SFR for all simulated galaxies at various times (see Table 1 and relate them to the FIR luminosity via equation 15). We find very LAT non-detected with AGN (upper limit) LAT non-detected (upper limit) LAT detected with AGN LAT detected best fit fit uncertainty calorimetric relation CR diff, 10 10 M CR diff, 10 11 M CR diff, 10 12 M CR adv, 10 11 M CR adv, 10 12 M CR adv, 10 10 M 10 7 10 8 10 9 10 10 10 11 10 12 L 8−1000 µm [L ] Fig. 3.-Correlation between the gamma-ray luminosity (L 0.1−100 GeV ) and the SFR, respectively the FIR luminosity (L 8−1000 µm ) of star forming galaxies. Upper limits on the observable gamma-ray emission by Fermi-LAT (open red symbols, Rojas-Bravo & Araya 2016) are contrasted to gamma-ray detections from star-forming galaxies only (solid black) and with AGN emission (filled red); data are taken from Ackermann et al. (2012), except for NGC 2146 (Tang et al. 2014) and Arp 220 (Griffin et al. 2016;Peng et al. 2016). See text for details on the FIR-to-SFR conversion. We overplot the emission of our simulated galaxies that only account for advective CR transport (light blue) and simulations in which we additionally follow anisotropic CR diffusion (blue). Different symbols indicate simulations of differently sized galaxy halos. In each simulation, we sample the SFR in logarithmic steps of e. Note that our simulations fall on the best-fit observational FIR-gamma-ray correlation (orange) and start to deviate from the calorimetric relation (dashed green) for small SFRs. similar gamma-ray luminosities in our models CR adv and CR diff, which both match the observed relation L γ /(erg s −1 ) = 8.9 × 10 27 (L 8−100 µm /L ) 1.12 (Rojas-Bravo & Araya 2016) for L 8−1000 µm > 10 9 L . The simulations appear to overpredict the gamma-ray luminosity at the lowest SFRs. The match may be improved by lowering the gas density in these halos or by adopting a realistic multi-phase ISM, which becomes more porous towards lower SFRs.
The relative contribution of hadronic losses to the total CR loss rate, Γ hadr /Γ tot , sets the normalization of the calorimetric relation, L γ ∝ L cr ∝ SFR (L cr is the CR luminosity). Star-forming galaxies with SFR 10 M yr −1 are close to the calorimetric relation (see Fig. 3). Instead, galaxies with lower SFRs fall below this relation, indicating that non-hadronic CR energy losses start to become relevant.
We explore different CR energy gain and loss processes in Fig. 4. In the Milky Way-like halo, the nonadiabatic CR loss processes (hadronic and Coulomb interactions) cool most of the CR energy that is injected; with cooling rates of Γ hadr = 7.44 × 10 −16 s −1 and Γ Coul = 2.78 × 10 −16 s −1 , respectively ). The (subdominant) adiabatic gains in the CR diff model are twice those of the CR adv model as CRs that are diffusing from the dense ISM into the halo are caught by the accretion flow and get advected back onto the disk. Hence, the calorimetric assumption for these high-mass galaxies is justified. This is in stark contrast to the dwarf galaxy simulation (10 10 M halo), where adiabatic and non-adiabatic losses are equally strong in the CR diff model. This is the main reason why these slowly simmering star-forming galaxies deviate from the calorimetric relation. In our dwarf-galaxy simulation in the CR adv model, the SFR levels off after 0.3 Gyr. This causes an almost linear increase of the CR injection rate that is counteracted by the increasing non-adiabatic CR cooling rates. Those dominate over the adiabatic CR losses at late times and thus move the lowest two (light-blue) simulation points in Fig. 3 closer to the calorimetric relation.
So far, we choose the canonical value for the CR injection efficiency at supernovae of ζ SN = 0.1. Varying this value by a factor of three results in a similar change in L γ (Table 1). In all these cases, the CR injection and adiabatic gains are balanced by the total CR cooling rate, indicating self-regulation. Hence, L γ and-by extension-the FIR-gamma-ray relation remain invariant if ζ SN decreases by a factor of two and the CR spectral index is increased to α = 2.15 (Fig. 2).

CONCLUDING REMARKS
For the first time, we have calculated the gamma-ray emission in galaxy simulations that span four orders of magnitude in SFRs and by modeling ideal MHD and CR physics self-consistently. We identify the influence of different CR gain and loss processes on the gammaray emission.
In agreement with previous literature, we find that the gamma-ray luminosity from decaying pions dominates over the IC emission by at least a factor of 5.7 in the Fermi-energy band 0.1 − 100 GeV. This dominance increases significantly for spectral indices steeper than α = 2.05. The continuous injection of CR energy at supernova remnants increases CR pressure to the point where the buoyancy force overcomes the magnetic tension of the dominant toroidal field after 1 Gyr. This enables CRs to diffuse into the halo and to accelerate the gas in a bubble-like, CR-dominated outflow. However, these bubbles are invisible in our simulated gamma-ray maps of hadronic pion-decay and inverse-Compton emission because of low gas density in the outflows. This suggests that morphological features such as the Fermi bubbles in the Milky Way may be generated by leptonic IC emission.
We find that most CR energy is lost to hadronic interactions at high SFRs. However, the gamma-ray emission deviates from this calorimetric property at low SFRs due to adiabatic losses, which cannot be identified in traditional one-zone models. Assuming that gamma rays mainly originate from decaying pions, we show that our simulated galaxies exactly reproduce the observed FIRgamma-ray relation at FIR luminosities L 8−100 µm > 10 9 L . This non-trivial finding comes about because we adopt a phenomenological relation between SFR and FIR emission and because L π 0 −γ ∝ f cal L cr ∝ f cal SFR where the calorimetric factor f cal smoothly decreases towards smaller SFRs due to the increasing importance of   . We contrast our model of advective CR transport (dashed) to the model in which we additionally follow anisotropic CR diffusion (solid). While CRs experience a modest adiabatic energy gain in the Milky Way-like galaxy, they suffer a substantial adiabatic loss in the dwarf galaxy which even matches the non-adiabatic CR loss in model CR diff. This causes the FIR-gamma-ray relation to deviate from the calorimetric relation. adiabatic losses. At small FIR luminosities (L 8−100 µm < 10 9 L ), the Magellanic Clouds emit less gamma rays in comparison to our simulated analogues. We speculate that this is either due to a too dense ISM or due to a missing multi-phase, porous ISM in our simulations, which would lower the hadronic gamma-ray yield at fixed gas mass, provided the CR density is (anti-)correlated with the thermal gas.