Galaxy Simulations


Background and Motivation

Cosmic rays (CRs) are extremely energetic particles that propagate throughout the universe, consisting of protons, electrons, heavy nuclei, etc. They can be accelerated to energies beyond \(10^{20}~{\rm eV}\) with possible sources including relativistic jets and powerful astrophysical shocks. Much effort has been made starting from the early twentieth century to present day in understanding their origin, properties, and role in various environments in the universe, from the interstellar medium (ISM), to galaxies, and to galaxy clusters. The study of CRs has aided in advances in a variety of fields, including particle physics, cosmology, high-energy astrophysics, and plasma astrophysics. Open questions still remain to this day, such as the origin of the ultra-high-energy CRs and the nature of particle acceleration. Continual advances, however, in numerical methods, high-performance computing architectures, observational technology, and laboratory-astrophysics techniques show promise for further unraveling the mysteries behind CRs.

The two main deviations from a pure power-law are the knee, where the energy spectrum steepens slightly near \(10^{15}~{\rm eV}\), and at the ankle, where the spectrum slightly flattens near \(10^{18}~{\rm eV}\). These deviations in the power-law may indicate a transition of different populations of CRs, where it is theorized that CRs below the knee originate in the Galaxy and CRs above the ankle are extragalatic in origin (Aharonian 2004, Beatty and Westerhoff 2009).

In this work, we will focus on the role of CRs in galactic environments. Galaxies are gravitationally bound systems of stars, interstellar gas, dust, and dark matter. They host a variety of physical processes and structures, such as star formation, supernovae, the ISM, and super-bubbles, which are expected to play a major role in CR acceleration (Bykov 2001). Super-bubbles are particularly interesting to test CR physics since they produce different CR spectra than supernovae (Ferrand and Marcowith 2010). The ISM contains on the order of \(10\%\) of the total mass of the gaseous disc and plays a significant role in CR physics. The global magnetic field within the ionized part of the ISM is responsible for CR diffusion. This process will be detailed in later sections. The interactions of CRs with the ISM produces secondary CRs of lower mass, which is an important phenomenon for testing CR propagation models. The gamma-ray emission from CR interactions traces the spatial distribution of the ISM. Of particular interest to CR physics is starburst galaxies, where the star-formation rates (SFRs) are significantly higher than the long-term average rate of the systems. The elevated SFR would result in high CR fluxes due to the subsequent supernovae in such galaxies. Also, starburst galaxies contain large amounts of dense, molecular gas, strong magnetic fields, and intense radiation fields. This combination of high CR populations and extreme environment leads to high radio and gamma-ray fluxes that are observed (Ackermann et al. 2012).


The dynamics of galactic evolution involves a variety of coupled physical processes, including star formation, gas dynamics, magnetic fields, and CRs. The propagation of CRs is an important mechanism in galactic environments for global evolution and for observations of galaxies at multiple wavelengths. The standard picture that is commonly accepted by the scientific community is summarized in this paragraph. CRs from \(1~{\rm GeV}\) to the knee originate in our Galaxy and are accelerated in supernova remnants (SNRs) via diffusive shock acceleration (DSA), where the energy spectrum that arises is a power-law. The accelerated particles propagate in the interstellar medium (ISM) and are deflected by the galactic magnetic field. Provided that the magnetic field is strong enough to keep CRs on small gyro-radii relative to galactic scales, CRs do not stream freely through the ISM. Also, since the galactic field has a turbulent component whose amplitude is comparable with the local field, CRs follow an erratic path and their motion is well-described by a diffusion equation. The correlation length is generally smaller than the Larmor radius since CRs diffuse more quickly at higher energies. This effect steepens the spectrum. While CRs propagate, they may be stochastically re-accelerated by the interaction with Alfvén waves and advected by winds. These effects are greater for CRs at lower energies. Observations will be important for constraining the diffusion coefficients of CRs in different galactic environments.

Once CRs are accelerated, they can propagate along magnetic-field lines. CR protons interact with the interstellar medium (ISM) via pion production \[{\rm p}_{\rm CR}^+ + {\rm p}_{\rm ISM}^+ \to {\rm p}^+ + {\rm p}^+ + \pi^0 \ ,\] where the pions rapidly decay into gamma rays (\(\pi^0 \to \gamma + \gamma\)) with energies \(> 100~{\rm MeV}\). CR electrons interact with both the ISM and interstellar radiation field (ISRF) via Bremsstrahlung and inverse Compton scattering, respectively

\[\begin{aligned} {\rm e}_{\rm CR}^- + {\rm p}_{\rm ISM}^+ &\to {\rm e}^- + {\rm p}^+ + \gamma \ , \\ {\rm e}_{\rm CR}^- + \gamma_{\rm ISRF} &\to {\rm e}^- + \gamma \ .\end{aligned}\]

CR electrons also produce synchrotron emission as they propagate along magnetic-field lines. The energy-loss rates for these processes are as described by Schlickeiser 2002:

\[\begin{aligned} -\left. \frac{{\rm d} E}{{\rm d} t} \right|_{\rm brem} &\propto n_{\rm ISM} E_{\rm e} \ , \\ -\left. \frac{{\rm d} E}{{\rm d} t} \right|_{\rm IC} &\propto U_{\rm rad} E_{\rm e}^2 \ , \\ -\left. \frac{{\rm d} E}{{\rm d} t} \right|_{\rm synch} &\propto U_{\rm mag} E_{\rm e}^2 \ . \end{aligned}\]

Thus, CR electrons at \({\rm GeV}\) energies lose the majority of their energy to bremsstrahlung. At higher energies, inverse Compton and synchrotron processes are comparable in significance, with the dominant one determined by the energy density in the magnetic and radiation fields.

Determining CR populations in our own Galaxy and in others is a challenging endeavor due to several difficulties, including the partial or complete confinement of CRs to galaxies and our inability to trace CRs back to their original sources. CRs at energies below the knee propagate along galactic magnetic field lines, which obscures any information about particle trajectories (Strong et al. 2007). For ultra-high-energy CRs, the angle at which they are deflected by galactic and extragalactic magnetic fields is extremely small (Beatty and Westerhoff 2009). While this may allow us to trace these CRs back to their sources, further steepening of the CR spectrum at around \(10^{19}~{\rm eV}\) due to the Greisen-Zatsepin-Kuzmin cutoff results in a detection of such particles only once per square kilometer per century. Due to these difficulties, we must rely on a variety of messengers to study galactic and extragalactic CR populations.


Generally, simulations of CRs either consist of kinetic particle-in-cell (PIC) simulations at small scales \(({\sim} 1~{\rm cm})\) or hydrodynamic simulations on galactic scales. We review the advances made in both regimes in this section.

Astrophysical particle acceleration is a long-standing problem in plasma astrophysics. DSA has been the main acceleration model (Reynolds 2008 and references therein); however, a major criticism of this model is that the maximum energy that CRs achieve via this mechanism is about 1 magnitude less than where the knee in the spectrum is (Lagage and Cesarsky 1983). Nonlinear effects, magnetic field amplification, etc., may also play a significant role in accelerating particles to relativistic energies (Aharonian 2004, Zweibel 2013).

Direct simulation of DSA requires expensive kinetic techniques and is still an active area of research. Hybrid particle-in-cell (PIC) simulations of CRs self-consistently coupled to magnetized fluids have emphasized their role in exciting flow instabilities. While it is computationally prohibitive to perform a PIC simulation on galactic scales, appropriate CR-MHD coupling can be achieved via a moment-equation approach. Results suggest that while the nonlinear back-reaction due to CRs modifies the shock structure, in realistic shocks with Mach numbers \({\sim} 10\), the shocks and CR spectra achieve steady state. Previous work, such as Skillman et al. 2008, has treated shock acceleration by using an algorithm that detects shocks, estimates the upstream Mach numbers and adds a local source term to the CR energy density. We do not consider DSA in this work, but future investigation will involve using astronomical observations to constrain subgrid acceleration models.

Previous studies of galactic evolution have focused on the dynamical effects of these physical processes, while simulated observations have not received as much attention. Modeling CR-driven galactic outflows has received significant attention in the context of investigating. Whether CRs can actively regulate processes in the ISM or drive winds and outflows strongly depends on their spatial distribution as well as their coupling with the gas. CR diffusion results from particle scattering on magnetohydrodynamic (MHD) waves and discontinuities. The resulting spatial diffusion is anisotropic locally and goes mostly along the magnetic-field lines; however, strong fluctuations of the magnetic field on large scales \(L\sim 100~{\rm pc}\), where the strength of the random field is several times higher than the average field strength, which makes CR diffusion on global Galactic scales primarily isotropic (Strong et al. 2007). The existence of galactic winds observed in several galaxies suggests that convective transport could play a significant role in CR transport (Strong et al. 2007). In addition to spatial diffusion, the scattering of CR particles on randomly moving MHD waves leads to stochastic acceleration which can be described by diffusion in momentum space.

Previous theoretical and computational modeling consists of steady-state and dynamic models, which have progressively increased in complexity. Breitschwerdt et al. (1991) considered 1-D steady wind models that incorporated CR streaming along large-scale coherent magnetic fields and pressure forces due to CR scattering from self-excited Alfvén waves. In their model, CRs stream along the large-scale flux tubes of coherent magnetic fields, oriented vertically close to the disc and radially far from the galaxy. The CRs exert a pressure on the thermal gas via scattering with the field and an Alfvén-wave pressure is also explicitly treated. This model suggested that CRs play an important role in accelerating thermal gas beyond the disc plane, with slowly outflowing to increasingly high speeds far above the star-forming disc.

Subsequent work has enhanced this simple model to include effects in the disc-halo interface, slightly modified magnetic-field geometries, disc rotation/magnetic tension, and more sophisticated damping mechanisms (Breitschwerdt et al. 1993, Zirakashvili et al. 1996, Everett et al. 2008). More recently, Dorfi and Breitschwerdt (2013) further extended these models to include time-dependent outflows and CR diffusion in 1-D. This extended approach found bursting, transient winds that featured forward- and backward-propagating shock fronts, with implications for DSA of CRs. These 1-D models have made a persuasive case for the dynamical importance of CRs, but can only treat aspects of the flow inherent to a uniform, coherent wind along ordered magnetic-field lines. More realistically, the dynamics are likely to involve multiphase, turbulent gas flows and magnetic-field lines with more complex topologies. In addition, all these models treat the flow of gas and CRs beyond the disc as an inner boundary condition and dod not explore how gas and CRs are produced within and rise out of the patchward star-forming regions of a real disc. Girichidis et al. (2016) performed 3-D MHD simulations of zoomed-in regions of galactic discs that treated anisotropic CR diffusion and a sophisticated treatment of chemical network to compute gas cooling. Their work demonstrated that CRs can drive a smooth cold wind and thicken galactic discs provided that the CRs are coupled to the gas.

CRs have an observed energy density of \({\sim}10^{-12}~{\rm ergs}~{\rm cm}^{-3}\) in the solar neighborhood (Wefel 1987), which makes their dynamical effect comparable to turbulent pressure and magnetic fields. The CR pressure is contributed primarily by protons with energies of the order \({\rm GeV}\) that scatter off of magnetic-field inhomogeneities. This explains the observed isotropic distribution despite the rarity of supernova, which are presumed to be sites of CR acceleration. CRs stream along their own pressure gradient, but streaming that occurs faster than the Alfven speed will excite waves in the magnetic field and heat the gas (Cesarsky 1980, Wentzel 1969,). It has often been argued that the CR-pressure contribution is unimportant to ISM dynamics since CR diffusion ought to rapidly diminish pressure gradients that could drive gas flows on global scales. Models, however, suggest that high star-formation rates can produce and maintain a local enhancement of CR pressure which can drive large-scale galactic winds that regulate star formation (McKenzie et al. 1987, Breitschwerdt et al. 1991, Socrates et al. 2008, Samui et al. 2010, Dorfi and Breitschwerdt 2013). Also, the role of diffusion diminishes on larger length scales, suggesting that CRs may be important for suppressing large-scale disc fragmentation, while still allowing smaller-scale molecular clouds to form.

CRs have recently been incorporated into 3-D global hydrodynamics simulations of galaxy evolution. Salem and Bryan (2014). Hydrodynamic simulations of CR diffusion in model galaxies have demonstrated the role of CRs in driving galactic outflows and winds (Booth et al. 2013, Salem and Bryan 2014). MHD simulations of CR diffusion in model galaxies have suggested the impact of CRs in driving a global magnetic dynamo in the host galaxy (Hanasz et al. 2004, 2006, 2009). These works have shown that weak magnetic fields can be amplified up to a few \(\mu {\rm G}\) within 1-2 Gyr by the dynamo process, which is driven by CR buoyancy as proposed by Parker (1992). Their model describes CRs as with the diffusion-advection equation supplemented to the resistive MHD equations. The CRs are supplied by SNRs and it is assumed that \(10\%\) of the supernova kinetic-energy output contributes to the CR energy density. The latest work on CR propagation in model galaxies was performed by Ruzkowski et al. (2016), where they have added streaming into their model and have obtained galactic outflows in agreement with Salem and Bryan (2014). Their model includes MHD, star formation, self-gravity, radiative cooling, anisotropic CR diffusion, and CR streaming along magnetic-field lines.


In terms of observations, CRs can either be observed indirectly through their radiation signatures or directly when they penetrate the Earth’s atmosphere and create particle showers. For this work, we focus on CR radiation signatures, which consist of radio, hard x-ray, and gamma-ray wavelengths.These non-thermal emission signatures provide insight into CR propagation and provide a test-bed for validating CR theories. Ground observations with HESS and VERITAS and space observations with the Fermi Gamma-Ray Space Telescope are making pioneering measurements of high-energy gamma-rays from galaxies, including nearby starburst galaxies (Acero et al. 2009, Abramoski et al. 2012, Ackermann et al. 2012). Observations from the Fermi have suggested that a significant portion of the galactic CRs originate from galactic supernova remnants (SNRs), which supports previous theoretical suggestions (Ginzburg and Syrovatskii 1964).

Analysis of these multi-wavelength observations provide unique insights into interactions between CRs and the ISM, which in turn lead to a better understanding of the distribution of feedback energy supplied by supernovae (Paglione et al. 1996, Torres 2004, de Cea del Pozo et al. 2009b, Lacki et al. 2012, Paglione and Abrahams 2012). Further progress in assessing the high-energy particle content of galaxies comes from radio measurements of the synchrotron emission produced by CR electrons. Radio observations are also advancing in terms of dynamic range, improved angular resolution, and sensitivity, especially at low frequencies. A combination of gamma-ray and radio spectra can thus be analyzed to yield information on the high-energy particle content and magnetic fields in galaxies (Thompson et al. 2006).

Laboratory experiments

The field of laboratory astrophysics, especially in the high-energy-density regime, is a relatively new field that has proven to be a powerful window into the transient dynamics of various astrophysical systems, including supernovae, jets, and accretion discs. Experimental advances have shown promise in advancing CR physics. Collisionless shocks have been studied in the laboratory, which has helped to improve our understanding of their role and impact in astrophysical observations (Fiuza et al. 2012). Experiments to study DSA have been thought out (Bell et al. 2012), with estimates for detection and what experimental designs would be optimal for such a detection.

More recently, particle acceleration at collisionless shocks has started to be investigated in the laboratory (Stockem et al. 2014). Such experiments propose to generate magnetized collisionless shocks and inject particles into them to try to initiate shock acceleration. Not only might these experiments provide useful validation benchmarks for PIC codes and CR models, they may allow us to understand the evolution of particle acceleration at collisionless shocks at a level which has previously not been possible.


CR physics has been a long-lived problem in the scientific community and major breakthroughs will likely require the collaboration of theory, simulation, observation, and experiment. In light of the advances in understanding CR propagation, we will focus on the non-thermal emission signatures of galaxies. To interpret multi-wavelength observations of galaxies in the context of galactic evolution requires detailed simulations that incorporate and couple the aforementioned physical processes. Ideally, such simulations would capture the fine-scale to large-scale effects of self-gravitating MHD coupled to energy-dependent multi-species CR transport.

Using the FLASH code, we will perform hydrodynamic simulations of CR propagation in model galaxies and produce simulated observations at gamma-ray wavelengths. As a starting point, we follow the model of Salem and Bryan (2014), where we assume that the CRs can be treated as a relativistic plasma that is advected by the thermal plasma and diffuses isotropically. The FLASH code currently treats CRs as a single relativistic fluid and evolves its total energy by solving an advection-reaction-diffusion equation (Rasera and Chandran 2008). A CR pressure terms is added to the MHD equations to account for the dynamical effects of CRs on the thermal gas. This model has been used for previous studies of the Fermi bubbles in the Milky Way (Yang et al. 2012) and has demonstrated efficacy.

The outline of this paper is as follows. In Section 2, we present our methodology, describe our model and discuss the physical processes that we account for in our simulations. We present and discuss our results in Section 3 and conclude in Section 4.