The Spatially Uniform Spectrum of the Fermi Bubbles: the Leptonic AGN Jet Scenario

TL;DR

This paper demonstrates that the spatially uniform gamma-ray spectrum of the Fermi bubbles can be explained by a leptonic jet scenario, where cosmic-ray electrons cool rapidly near the GC and then advect outward, maintaining a consistent spectrum.

Contribution

The study introduces a new simulation approach with a CR spectrum tracking module, showing how leptonic jets can produce the observed uniform spectrum of the Fermi bubbles.

Findings

High-energy cutoff caused by rapid synchrotron and IC cooling near the GC.

Spectrum remains uniform due to advection after initial cooling.

Leptonic jet model matches observed gamma-ray features.

Abstract

The Fermi bubbles, two giant gamma-ray bubbles above and below the Galactic center (GC), are among the most important findings of the Fermi Gamma-ray Space Telescope. Because of the proximity, spatially resolved, multi-wavelength observations offer excellent opportunities to learn about the physical origin of the bubbles. One of the unique features of the observed bubbles is that their gamma-ray spectrum, including a high-energy cutoff at ~110 GeV and the overall shape of the spectrum, is nearly spatially uniform. The high-energy spectral cutoff is suggestive of a leptonic origin as it could be a signature of synchrotron and inverse-Compton (IC) cooling of cosmic-ray (CR) electrons; however, even for a leptonic model, it is not obvious why the spectrum should be spatially uniform. In this work, we investigate the formation of the Fermi bubble in the leptonic jet scenario using a newly implemented CRSPEC module in FLASH that allows us to track the evolution of CR spectrum on-the-fly during the simulations. We show that the high-energy cutoff is caused by fast synchrotron and IC cooling near the GC when the jets were first launched. After the jets propagate away from the GC, the dynamical timescale of the jets become the shortest among all relevant timescales, and therefore the spectrum is essentially advected with only mild cooling losses. This could explain why the bubble spectrum is nearly spatially uniform: the CRs from different parts of the bubbles as seen today all share the same origin near the GC at early stage of the bubble expansion. We find that the CR spatial and spectral distribution predicted by the leptonic jet model can simultaneously match the normalization, spectral shape, and high-energy cutoff of the observed gamma-ray spectrum and their spatial uniformity, suggesting that past AGN jet activity is a likely mechanism for the formation of the Fermi bubbles.

Equations

$$\frac{dN}{d{E}_{\gamma }d{E}_{\mathrm{ph}}dt}$$

$$\Gamma$$

$$q$$

$$E_{\mathrm{ph}}$$

$$⟨E_{\gamma }⟩=(4/3){\gamma }^2⟨E_{\mathrm{ph}}⟩.$$

$$f(p)=f_i{\left(\frac{p}{p_{i-1}}\right)}^{-q_i},$$

$$v_{\mathrm{jet}}<0.065c\left(\frac{u_{\mathrm{tot}}}{10^{-11}\mathrm{erg}\mathrm{c}{\mathrm{m}}^{-3}}\right)\left(\frac{E_{\max ,0}}{\mathrm{TeV}}\right),$$

$$E\sim \frac{1}{\beta t}\sim 2.5\mathrm{TeV}{\left(\frac{u_{\mathrm{tot}}}{10^{-11}\mathrm{erg}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-1}{\left(\frac{{\tau }_{\mathrm{dyn}}}{0.018\mathrm{Myr}}\right)}^{-1}.$$

$$v_{\mathrm{jet}}>0.02c\left(\frac{u_{\mathrm{tot}}}{10^{-11}\mathrm{erg}\mathrm{c}{\mathrm{m}}^{-3}}\right)\left(\frac{E_{\max }}{300\mathrm{GeV}}\right).$$

$$\frac{\partial f}{\partial t}=-\mathbf{v}\cdot \nabla f+\nabla \cdot (\kappa \nabla f)+\frac{1}{3}(\nabla \cdot \mathbf{v})p\frac{\partial f}{\partial p}+\frac{1}{p^2}\frac{\partial }{\partial p}\left[p^2\left(b_l(p)f+D_{pp}\frac{\partial f}{\partial p}\right)\right]+j(\mathbf{x},p),$$

$$f(p)=f_i{\left(\frac{p}{p_L}\right)}^{-q_i},p_L=\{\begin{array}{ll}{p}_{i-1},& \mathrm{if}{p}_{i-1}\ge p_{cutL}\\ p_{cutL},& \mathrm{otherwise}\end{array}$$

$$\frac{\partial n_i}{\partial t}$$

$$n_i$$

$$⟨\kappa {⟩}_i$$

$$Q_i$$

$$\frac{\partial n_i}{\partial t}={\left[b(p)4\pi p^{2}f(p)\right]}_{p_{i-1}}^{p_i},$$

$$b(p)\equiv \frac{dp}{dt}=\frac{1}{3}(\nabla \cdot \mathbf{v})p+b_l(p),$$

$$n_i^{t+\Delta t}$$

$${\Phi }_i^p$$

$${\Phi }_i^p=-\frac{4\pi }{\Delta t}{\int }_{p_i}^{p_u}{p}^2{f}_j(p)dp,$$

$$j=\{\begin{array}{ll}i+1,& \mathrm{if}b(p_i)>0\\ i,& \mathrm{if}b(p_i)\le 0,\end{array}$$

$$\Delta t={\int }_{p_i}^{p_u}\frac{dp}{b(p)}.$$

$$\frac{\partial e_i}{\partial t}$$

$$e_i$$

$$⟨{\kappa }_T{⟩}_i$$

$$S_i$$

$$\frac{\partial e_i}{\partial t}={\left[b(p)4\pi p^{2}f(p)T(p)\right]}_{p_{i-1}}^{p_i}-{\int }_{p_L}^{p_R}b(p)\frac{4\pi p^{3}cf(p)}{\sqrt{m_e^2{c}^2+p^2}}dp.$$

$$R_{i}=\int_{p_{L}}^{p_{R}}b(p)\frac{p^{3-q_{i}}}{\sqrt{m_{e}^{2}c^{2}+p^{2}}}dp\bigg{/}\int_{p_{L}}^{p_{R}}p^{2-q_{i}}T(p)dp.$$

$$e_i^{t+\Delta t}\left(1+\frac{\Delta t}{2}{R}_i\right)=e_i^t\left(1-\frac{\Delta t}{2}{R}_i\right)-\Delta t({\Phi }_i^e-{\Phi }_{i-1}^e),$$

$${\Phi }_i^e=-\frac{1}{\Delta t}{\int }_t^{t+\Delta t}b(p)4\pi p^2{f}_i(t^{\prime },p)T(p){|}_{p_i}d{t}^{\prime },$$

$${\Phi }_i^e=-\frac{4\pi }{\Delta t}{\int }_{p_i}^{p_u}{p}^2{f}_j(p)T(p)dp,$$

$$\Delta t={\int }_{p_{cut}^{t+\Delta t}}^{p_{cut}^t}\frac{dp}{b(p)}.$$

$$\frac{e_i}{n_i{p}_{L}c}=\frac{3-q_i}{(p_R/p_L{)}^{3-q_i}-1}\frac{(p_R/p_L{)}^{4-q_i}-1}{4-q_i},$$

$${\left(\frac{dE}{dt}\right)}_{syn+IC}=\frac{4}{3}c{\sigma }_T{\beta }^2{\gamma }^2(u_B+u_{rad}{F}_{KN}),$$

$$F_{KN}=\frac{1}{u_{rad}}{\int }_{E_{ph,min}}^{E_{ph,max}}{f}_{KN}(x)E_{ph}n(E_{ph})d{E}_{ph},$$

$$f_{KN}(x)\simeq \frac{1}{(1+\Gamma {)}^{1.5}}\mathrm{for}\Gamma =\frac{4\gamma E_{ph}}{m_e{c}^2}\lesssim 10^4$$

$${\tau }_{syn+IC}=0.97\mathrm{Myr}{\left(\frac{u_B+u_{rad}{F}_{KN}}{10^{-12}\mathrm{erg}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-1}{\left(\frac{\gamma }{10^6}\right)}^{-1}.$$

Full text

1. Introduction

The Fermi bubbles, two giant bubbles extending 50 degrees above and below the Galactic center (GC), are among the most important findings of the Fermi Gamma-ray Space Telescope (Su et al., 2010; Ackermann et al., 2014; Narayanan & Slatyer, 2017). The observed gamma-ray bubbles have many unique charateristics, including the spatially uniform hard spectrum, nearly flat surface brightness distribution, sharp edges, and smooth surface. The bubbles are also spatially coincident with features in other wavelengths, such as the microwave haze (Finkbeiner, 2004; Planck Collaboration, 2013), X-ray properties of the Galactic halo (e.g., Snowden et al., 1997; Bland-Hawthorn & Cohen, 2003; Kataoka et al., 2013; Tahara et al., 2015; Kataoka et al., 2015; Miller & Bregman, 2016), UV absorption lines (Fox et al., 2015; Bordoloi et al., 2017), and polarized lobes (Carretti et al., 2013). Because of the proximity, the spatially resolved, multi-wavelength observational data provides unprecedented opportunities for studying the physical origin of the bubbles as well as cosmic ray (CR) propagation, Galactic magnetic field, and past activity at the GC.

Many theoretical models have been proposed to explain the formation of the bubbles. The hard spectrum of the observed bubbles implies that the CRs, if transported from the GC, must reach large distances before they have time to cool. This consideration gives a constraint on the age of the bubbles to be a few Myr if the gamma rays are produced by CR electrons (CRe, i.e., the leptonic model). In order to satisfy the age constraint, the theories can be divided into three categories: hadronic transport (e.g., Crocker & Aharonian, 2011; Mou et al., 2014; Crocker et al., 2015), leptonic transport (e.g., Guo & Mathews, 2012; Guo et al., 2012; Yang et al., 2012, 2013, the latter two are abbreviated as Y12 and Y13, respectively), and in-situ acceleration models (e.g., Mertsch & Sarkar, 2011; Cheng et al., 2011, 2015; Sarkar et al., 2015; Sasaki et al., 2015). In the hadronic transport models, the gamma rays are generated by inelastic collisions between CR protons (CRp) and thermal nuclei. The CRp are produced at the GC by nuclear starburst or activity of the central active galactic nucleus (AGN), and they are subsequently transported via starburst or AGN driven winds. The hadronic models can successfully reproduce the properties of the observed gamma-ray bubbles; however, to model the microwave haze is nontrivial (Ackermann et al., 2014) and requires an additional population of primary CRe (Crocker et al., 2015). In the leptonic transport models, CRe are injected at the GC via past jet activity of the central supermassive black hole (SMBH) and transported by fast AGN jets. Previous simulations have shown that the bubbles can be inflated within a few Myr (Guo & Mathews, 2012, Y12). Also, the key features of the gamma-ray bubbles as well as the microwave and polarization signatures are in good agreements with the observational data (Y12, Y13). Some observational studies of the thermal and kinematic properties of the Galactic halo (e.g., Kataoka et al., 2013; Sarkar et al., 2017) suggest that the Fermi bubbles are triggered by milder outflows, which could potentially be in tension with the jet model. However, more data is needed to draw a conclusion because the Galactic halo in the vicinity of the bubbles is extremely complex (Tahara et al., 2015; Kataoka et al., 2015). Also, there are discrepancies among observationally derived kinematics (e.g., Sarkar et al., 2017; Bordoloi et al., 2017), possibly due to modeling uncertainties such as the assumptions of geometry and injection pattern of the outflows. For the in-situ acceleration models, CRs are assumed to be produced by shocks or turbulence near the edges of the bubbles. Although these models could bypass the age constraints, it has been challenging for the simplest models to produce the flat gamma-ray intensity profile (Mertsch & Sarkar, 2011; Cheng et al., 2011) as well as the microwave haze emission (Fujita et al., 2014; Cheng et al., 2015).

One unique and important feature of the observed bubbles that has not been investigated in detail is the spatially uniform gamma-ray spectrum. Ackermann et al. (2014) showed that the bubble spectrum can be well fit by a power law with an exponential cutoff at $\sim 110$ GeV. Remarkably, both the shape of the spectrum and the cutoff energy are almost independent of Galactic latitude (see also Narayanan & Slatyer, 2017). The high-energy cutoff is suggestive of a leptonic origin because CRe can cool more easily due to synchrotron and inverse-Compton (IC) energy losses. However, even for a leptonic model, it is unclear why the spectrum should be spatially uniform.

In Y12 and Y13, we investigated the leptonic AGN jet scenario using three-dimensional (3D) magnetohydrodynamic (MHD) simulations including relevant cosmic-ray (CR) physics. As mentioned above, the leptonic jet model is a promising mechanism for explaining the origin of the bubbles as it is the simplest model that could simultaneously explain the gamma-ray bubbles and the microwave haze. However, in the previous works, the CRs are treated as a single fluid without distinguishing their energies, and therefore comparisons with observations have to rely on assumptions of the CR spectrum. In this study, we implement a new CRSPEC module in the FLASH code (Fryxell et al., 2000; Dubey et al., 2008) that could handle CRs of different energy channels and follow their spectral evolution on-the-fly during the simulations. We apply it to simulate the spectral evolution of the CRs within Fermi bubbles and generate the gamma-ray spectrum self-consistently. Our objectives are to answer the questions: (1) what physical mechanisms are responsible for the $\sim 110$ GeV cutoff in the observed gamma-ray spectrum? (2) why is the bubble spectrum spatially uniform, including both the overall spectral shape and the cutoff energy?

The structure of the paper is as follows. We first discuss expectations of the CR spectra as hinted by the gamma-ray data in § 2. In § 3 we outline the simulation setup and describe key aspects of the CRSPEC code. In § 4 we present results from the simulations, including the simulated distribution of CR energies (§ 4.1), general spectral evolution of the bubbles (§ 4.2), the spatial dependence of the spectrum (§ 4.3), and constraints on the AGN jet speed, magnetic field strength, and energy density of the interstellar radiation field (ISRF) derived from our model (§ 4.4). Finally we summarize our findings in § 5.

2. Hints from the observed spectrum

The observed gamma-ray spectrum of the Fermi bubbles is strikingly latitude independent, including its shape and the high-energy cutoff (e.g., Fig. 33 in Ackermann et al. (2014)). The observed bubble spectrum can be best fit by a power law with an exponential cutoff term $\exp(-E/E_{\rm cut}$), where $E_{\rm cut}\sim 110$ GeV. This energy scale represents where the gamma-ray spectrum has a turnover; however, there could still be gamma rays generated beyond this energy (see Figure 5). In order to connect with the energies of the underlying CRs more directly, hereafter we define a “maximum energy of the observed gamma rays” as $E_{\rm max,obs}$, which relates to $E_{\rm cut}$ by $\exp(-E_{\rm max,obs}/E_{\rm cut})=0.1$ (i.e., the energy scale where the gamma-ray intensity is dimmer by a factor of 10). Given $E_{\rm cut}=110$ GeV, we have $E_{\rm max,obs}\sim 250$ GeV. In this section we show that the data alone can readily provide some clues about the underlying CR spectra and their latitude dependence. Specifically, for a given latitude bin, the shape of the gamma-ray spectra can inform the characteristic energy of the CRe, and the observed cutoff energy is related to the maximum energy of the CRe.

In the leptonic scenario, the gamma rays originate from IC scattering of the ISRF by CRe. The spectrum of the up-scattered photons per one electron of Lorentz factor $\gamma$ is given by Blumenthal & Gould (1970),

[TABLE]

where $E_{\rm ph}$ is the initial photon energy, $\gamma=E_{\rm e}/(m_{e}c^{2})$ is the Lorentz factor of the CR electron, $E_{\gamma}$ is the energy of the up-scattered gamma-ray photon, and $n(E_{\rm ph})$ is the energy distribution of the photon number density.

In the Thomson limit ($\Gamma\ll 1$), the average energy of the up-scattered photons is given by

[TABLE]

In the Klein-Nishina (KN) limit ($\Gamma\gg 1$), almost all the energy of the CR electron is carried away by the up-scattered photons, i.e., $\langle E_{\gamma}\rangle\sim\langle E_{\rm e}\rangle$. Assuming the CR spectrum is a power law with spectral index $\alpha$, it can be shown that the spectral index of the up-scattered gamma-ray photons is $(\alpha+1)/2$ in the Thomson limit, and approximately $\alpha+1$ in the KN limit (Blumenthal & Gould, 1970). The observed bubble spectrum is best fit by a power-law CR distribution with spectral index of $\sim 2$, and therefore one may expect to see changes in the spectral indices from 1.5 to 3 in the gamma-ray spectrum as the IC scattering goes from the Thomson limit to the KN regime.

The observed spectrum of the bubbles is nearly latitude independent, characterized by a broad bump that roughly peaks around $E_{\rm bump}\sim 10$ GeV (see Figure 5). This is not straightforward to obtain because the ISRF is dominated by the cosmic microwave background (CMB) at high latitudes and optical starlight at low latitudes. If the underlying CRe had identical spectra across all latitudes, then the resulting gamma-ray spectra would peak at lower energies at higher latitudes. This implies that the average energy of the CR population must be latitude dependent. In fact one could estimate the average energy of CRe for different latitude bins in the Thomson limit using Eq. 5 because $\Gamma<1$ for $E_{\gamma}\sim 10$ GeV and $E_{\rm ph}<10$ eV.

For high latitudes (e.g., $b=40^{\circ}-60^{\circ}$), the ISRF peaks at $\langle E_{\rm ph}\rangle\sim 7\times 10^{-4}$ eV for CMB photons. It would therefore require an average CR energy of $\langle E_{\rm e}\rangle\sim 2$ TeV in order to produce a $\sim 10$ GeV bump. For intermediate latitudes ($b=20^{\circ}-40^{\circ}$), the intensity of the ISRF is more uniform across all wavelengths. Assuming the gamma-ray bump primarily comes from infrared (IR) photons ($E_{\rm ph}\sim 10^{-2}$ eV), one would obtain $\langle E_{\rm e}\rangle\sim 200$ GeV. Similarly, the average CR energy can be estimated to be $\langle E_{\rm e}\rangle\sim 20$ GeV at low latitudes where optical light (assuming $\langle E_{\rm ph}\rangle\sim 5$ eV) dominates the ISRF. Therefore, generally speaking, the spatially uniform spectra of the Fermi bubbles require that the average energy of CRe to be higher at higher latitudes. The exact magnitude of the energy gradient, though, may be different from the above estimate because the observed gamma-ray bump is broad and $E_{\rm bump}$ does not have to be close to 10 GeV. In fact, the gradient of CR energies should be smaller in order to be consistent with the maximum energy of CRe as estimated below.

On the other hand, the maximum energy of the observed gamma-ray spectrum, $E_{\rm max,obs}$, comes from up-scattered optical light in the ISRF and provides information about the maximum energy of the underlying CR electron population, $E_{\rm max}$. For optical photons ($\langle E_{\rm ph}\rangle\sim 5$ eV), the IC scattering would be in the KN limit for CRe with energies greater than $\sim 100$ GeV. Therefore, for high and intermediate latitudes, $E_{\rm max}\gtrsim E_{\rm max,obs}\sim 250$ GeV as it is in the KN regime. For low latitudes, the average CR energy is smaller and hence the formula in the Thomson limit applies. Assuming $E_{\rm max,obs}=250$ GeV, Eq. 5 gives $E_{\rm max}\sim 100$ GeV. The observed cutoff energy is almost independent of latitudes, with a slight tendency of higher $E_{\rm max,obs}$ for higher latitudes (Figure 33 of Ackermann et al. (2014)). Therefore, the estimates above imply that $E_{\rm max}$ is also nearly spatially uniform, on the order of a few hundred GeV, and may be somewhat greater at higher latitudes. Our simple estimates of the CR electron cutoff energy for different latitudes are consistent with best-fit values to the observed gamma-ray spectrum obtained by Narayanan & Slatyer (2017).

3. Methodology

We simulate the spectral evolution of the Fermi bubbles in the leptonic AGN jet scenario using 3D hydrodynamic simulations including CRs. The simulation setup is essentially identical to that of Y12 and Y13, to which we refer the readers for details including the initial conditions for the Galactic halo as well as parameters for the AGN jets. Here we briefly summarize our approach and emphasize the differences from the previous works.

The simulations are performed using the adaptive-mesh-refinement (AMR) code FLASH (Fryxell et al., 2000; Dubey et al., 2008). Same as in Y12, the CRs are injected at the GC during a short (0.3 Myr), active phase of the central SMBH about 1.2 Myr ago and then the CRs are advected with the AGN jets. CR diffusion is omitted in the present work since it is slow enough that it only affects the sharpness of the bubble edges but not the overall dynamics and CR distribution (Y12). However, we return to this point in § 4.3. As in Y12, we assume that CRs are scattered by extrinsic turbulence rather than self-excited Alfven waves, and therefore the effects of CR streaming is not included. Since only a small fraction of injected CRs is needed to reproduce the gamma-ray signal (Y13), we assume $3\times 10^{-3}$ of the injected CRs to be CRe 111We note that the normalization factor is different from that adopted in Y13 because the injected the CR spectrum has a different spectral range. and follow their spectral evolution due to adiabatic compression, adiabatic expansion, synchrotron losses, and IC losses. The rest of the injected CR energy density 222Though called CR energy density, it is degenerate with the thermal energy density since the dynamics is determined by the total energy density of the jets (Y12). does not cool and is the dominating component in terms of dynamics. Passively evolving tracer particles are injected along with the jets in order to track the evolution of the bubble spectrum. We note that Su & Finkbeiner (2012) and Ackermann et al. (2014) have found substructures in the intensity distribution of the south bubble (i.e., the “cocoon”), which might be related to a second event of energy injection from the GC. However, in this work we do not consider CR injections from a second AGN outburst in order to avoid introduction of an additional set of jet parameters that are not uniquely constrained. To this end, we refrain ourselves from interpreting the substructures and only focus on the primary bubble emission that has nearly flat intensity distribution.

In Y13, we demonstrated that the magnetic field within the bubbles has to be amplified to values comparable to the ambient field in order to simultaneously produce the microwave haze. Therefore, we do not include magnetic fields in the current simulations but simply assume the default magnetic field distribution as in GALPROP (Strong et al., 2009), $|B|=B_{0}\exp(-z/z_{0})\exp(-R/R_{0})$, for the computation of synchrotron losses, where $R$ is the projected radius to the Galaxy’s rotational axis. We adopt $z_{0}=2$ kpc and $R_{0}=10$ kpc, which are best-fit values in the GALPROP model to reproduce the 408 MHz synchrotron radiation in the Galaxy. The normalization of the magnetic field strength $B_{0}$ is treated as a free parameter and the simulations presented in this paper has a fiducial value of 10 $\mu{\rm G}$. As we will discuss in § 4, the cutoff energy of the gamma-ray spectrum is very sensitive to $B_{0}$ and therefore could be used to put constraints on the initial conditions. Note that $B_{0}$ represents the magnetic field strength at the GC right after the initial injection, and therefore does not need to be the same as the field strength as observed today. In fact, $B_{0}$ is likely smaller than the present observed field strength at the GC (e.g., Crocker et al., 2010) because it takes time for the magnetic field within the bubbles to amplify after the initial adiabatic expansion caused by the jets (Y13).

For IC losses, we adopt the ISRF model from GALPROP v.50 (Strong et al., 2007) and compute the CR energy losses and gamma-ray emissivity including the KN effects (Jones, 1968). While there exist other ISRF models that are more general to all spiral galaxies (e.g., Popescu & Tuffs, 2013), we chose the GALPROP model because it is calibrated using stellar and dust distributions specific to the Milky Way. Also, it is adopted by all previous studies of the Fermi bubbles, allowing us to make comparisons to previous results directly. The ISRF model in GALPROP v.50 provides a 3D distribution of photon energy densities for discrete values of $(x,y,z)$ with spacings of 0.1 kpc. We therefore used 3D linear interpolations to obtain the photon energy densities on our simulation grid. The adopted ISRF decreases away from the GC. Specifically, the values range from $\sim 19\ {\rm eV\ cm}^{-3}$ near the GC to $\sim 1\ {\rm eV\ cm}^{-3}$ at 5 kpc away from the Galactic plane near the rotational axis.

3.1. Modeling CR spectral evolution

The core of this work is the newly implemented CRSPEC module in FLASH. Advection of CRs, anisotropic diffusion (though not used in this work), and dynamical coupling between the CRs and the gas are done in the same way as in the previous version (see equations in Y12). But instead of having a single equation for the evolution of the total CR energy density, the CRs are divided into $N_{p}$ logarithmically spaced momentum bins. Equations are solved for the CR number densities $n_{i}$ and CR energy densities $e_{i}$ in each bin with index $i$. The algorithm for solving this set of equations is based on the method for fast cooling electrons in the COSMOCR code (Miniati, 2001), and we made modifications in order to handle finite spectral ranges. In the adopted approach, the CR distribution function as a function of momentum, $p$, is approximated with a piece-wise power law,

[TABLE]

where $f_{i}$ and $q_{i}$ are the normalization and logarithmic slope for the $i$th momentum bin. Fluxes of $n_{i}$ and $e_{i}$ across different momentum bins are computed according to adiabatic processes and synchrotron and IC energy losses. For completeness we summarize the details, relevant equations, and test cases in the Appendix.

We inject CRe with a power-law spectrum from the GC in the beginning of the simulations for a duration of 0.3 Myr. The initial spectrum ranges from 10 GeV to 10 TeV, with a spectral index of $q_{i}=4.1$. The normalization factor in each bin $f_{i}$ is chosen so that the energy density of CRe is $3\times 10^{-3}$ of the total CR energy density of the jets, or equivalently, $7.5\times 10^{-12}\ {\rm erg\ cm}^{-3}$. A representative value for the normalization factor at 10 GeV is $2.24\times 10^{-4}\ {\rm cm}^{-2}\ {\rm s}^{-1}\ {\rm GeV}^{-1}$. To simulate the spectrum we have five logarithmically spaced momentum bins between 0.1 GeV and 10 TeV. The range is chosen so as to cover the energy shifts of the injected CRs due to the adiabatic and cooling processes. We use a relatively small number of momentum bins in order to minimize computational costs. This is adequate for our application because the CR spectrum, except for the high-energy end, only experiences advection and adiabatic processes, and therefore the spectrum can be well approximated by a power law. In order to accurately simulate the cutoff energy of the CR spectrum at the high-energy end due to fast electron cooling, which is one of the main purposes of the paper, we store extra variables for the minimum and maximum of the CR spectrum ($p_{cutL}$ and $p_{cutR}$, respectively) and track their evolution separately. In order to account for fast synchrotron and IC cooling of CRe accurately, the simulation timestep is set to 0.1 times the cooling timescale. When this timestep is the shortest among all relevant timescales in the simulation, we subcycle over the CR spectral evolution in order to accelerate computations. As a test of our algorithm, we performed a simulation of the Fermi bubble spectral evolution including only adiabatic processes (synchrotron and IC losses are turned off). We verified that the total CR number density is conserved after the jets are shut off at $t=0.3$ Myr, the spectrum is shifted but the shape is unaltered, and that the total CR energy density distribution at the end of the simulations, $t=1.2$ Myr, is identical to what was obtained using the energy-integrated version of the code as in Y12.

4. Results

4.1. Distribution of CR energies

Figure 1 shows distributions of the CRe at the end of the simulation, $t=1.2$ Myr, including the total CR electron energy density, total CR electron number density, maximum CR electron energy 333The maximum CR electron energy, $E_{\rm max}$, is solved separately for all cells with nonzero CR energy densities, and therefore the values are computed for an extended region beyond the bubbles where a tiny but nonzero amount of CRs exist due to numerical diffusion. We determined that this is a numerical artifact and therefore only plotted $E_{\rm max}$ for regions where the total energy density of CRe is greater than $10^{-16}\ {\rm erg}\ {\rm cm}^{-3}$, which corresponds to a minuscule level of $\sim 6\times 10^{-5}\ {\rm eV\ cm^{-3}}$., and CR electron energy densities in different energy channels. The first and last energy bins (i.e., 0.1-1 GeV and 1-10 TeV) are not shown because they contain little amount of CRe at the end of the simulations.

Because the CRe only contribute to $3\times 10^{-3}$ of the total CR energy density of the jets, they are not dynamically dominant and therefore the overall distribution of the CR electron energy density is similar to that in the adiabatic simulation (Y12, Figure 1). Though some of the details (e.g., structures close to the GC) are slightly different due to cooling of CRe, the main characteristics such as the bubble morphology and the edge-brightened distribution are recovered. We recall that the edge-enhanced CR distribution, which is a result of compression of jet materials during the active phase of the AGN injections, plays a crucial role in reproducing the flat surface brightness distribution of the observed bubbles after line-of-sight projection. In particular, the CR number density near the top of the bubbles must be greater by the right amount so that, after convolving with the ISRF whose intensity decreases with Galactic latitudes, the projected gamma-ray intensity is almost spatially uniform (Y13).

One thing to note from Figure 1 is that the total energy density (upper left panel) and total number density (upper middle panel) of CRe have similar distributions. Since $e_{\rm cr}\sim\langle E_{\rm e}\rangle n_{\rm cr}$, where $\langle E_{\rm e}\rangle$ is the characteristic energy of the CRe, this implies that $\langle E_{\rm e}\rangle$ is largely spatially uniform. This is indeed suggested by the map of maximum CR electron energy $E_{\rm max}$ within the bubbles (upper right panel): while both the CR electron energy and number densities differ by about two orders of magnitudes from minimum to maximum, the variation in $E_{\rm max}$ is relatively small, with $E_{\rm max}$ on the order of a few hundreds of GeV for most regions within the bubbles.

Although $\langle E_{\rm e}\rangle$ and $E_{\rm max}$ are generally quite spatially uniform, they do exhibit some gradients. The top right panel in Figure 1 shows that $E_{\rm max}$ varies from $\sim 100$ GeV at low Galactic latitudes to $\sim 1$ TeV at high latitudes within the bubbles. There is also a thin shell of very high energy CRs at the edges of the bubbles. However, since they occupy only a very small volume, and that the CR number density is low in this shell, their contribution to the projected CR spectra is negligible (see Figure 4). The gradient in CR energies is also evident by comparing the CR energy densities divided into different energy channels (bottom row in Figure 1). The overall uniformness and mild gradient toward higher energies at higher latitudes are consistent with the expectations derived from the observed bubble spectrum (see § 2). In the following section, we describe the evolution of the CR spectrum in order to understand the final distribution of CR energies as seen in Figure 1.

4.2. Spectral evolution of the Fermi bubbles

Figure 2 shows the evolution of one representative tracer particle that was injected at early stage of the bubble expansion. Panels from top to bottom show the evolution of the maximum energy of the CR spectrum ($E_{\rm max}$), minimum spectral energy ($E_{\rm min}$), relevant timescales, and divergence of the velocity field. The dotted line overplotted in the top panel represents the result from the adiabatic simulation without synchrotron and IC cooling. For this adiabatic case, the change in $E_{\rm max}$ and $E_{\rm min}$ are directly proportional to each other, meaning that the CR spectrum is only shifted without changing the spectral shape. Right after the particle was injected ($t\sim 0.1$ Myr), there was a brief increase in $E_{\rm min}$ and $E_{\rm max}$ due to adiabatic compression, i.e., $\nabla\cdot{\bm{v}}<0$. Afterwards, the CRe propagate outward and the only cooling mechanism is adiabatic expansion ($\nabla\cdot{\bm{v}}>0$), and therefore both $E_{\rm min}$ and $E_{\rm max}$ decrease with time monotonically.

In contrast, the evolution of $E_{\rm max}$ is quite different for the simulation including synchrotron and IC cooling. As shown in the top panel in Figure 2 (solid line), the maximum energy of the CRe drops rapidly from the injected energy of 10 TeV to $\sim 1$ TeV before $t\sim 0.3$ Myr. This fast change in $E_{\rm max}$ is owing to synchrotron and IC energy losses, as during this early phase of evolution the cooling timescale for synchrotron and IC radiation near the GC is much shorter than all other relevant timescales (see the third panel). After $t\sim 0.4$ Myr, the dynamical timescale ($\tau_{\rm dyn}\equiv(1\ {\rm kpc})/v$) becomes shorter/comparable to the synchrotron and IC cooling timescale, while the timescale for adiabatic processes ($\tau_{\rm AD}\equiv 1/(\nabla\cdot{\bm{v}}$)) is subdominant. That is, at later stage of the bubble expansion, the CRe experience advection (which does not cause energy losses) and synchrotron plus IC cooling at the same time (which now occurs on longer timescales compared to the beginning). Therefore, $E_{\rm max}$ only decreases slightly after $t\sim 0.4$ Myr and reaches a value about 700 GeV at $t=1.2$ Myr. In other words, the value of $E_{\rm max}$ at the present day is closely related to fast cooling of CRe near the GC at early stage of the bubble expansion when the jets were first launched.

4.3. Why is the spectrum spatially uniform?

In this section, we provide explanations as to why the gamma-ray spectrum of the Fermi bubbles is almost spatially uniform, including the maximum energy $E_{\rm max,obs}\sim 250$ GeV and the overall shape of the spectra for different latitude bins.

Figure 3 shows the evolution of a selection of tracers that have different final locations within the bubbles. Their distance to the Galaxy rotational axis, vertical height to the Galactic disk, maximum energy of the CR spectrum, and relevant timescales are plotted in the panels from top to bottom. We find that, although the particles all have distinct trajectories, their evolution of $E_{\rm max}$ is very similar, marked with a fast decay within the early $\sim 0.3$ Myr after injection and a subsequent mild decrease to several hundreds of GeV at the end of the simulation. Similar to the tracer shown in Figure 2, the CRe encounter significant energy losses due to synchrotron and IC cooling near the GC soon after they are injected with the AGN jets ($\tau_{\rm syn+IC}\ll\tau_{\rm dyn}$ as shown in the bottom panel). Afterwards, the dynamical time of the jets becomes shorter than adiabatic, synchrotron and IC cooling timescales, and hence the CR spectrum is essentially advected with only mild cooling losses. Though somewhat dependent on the time of injection and degree of initial compression, the final value of $E_{\rm max}$ for each tracer particle is similar because the energy scale is mainly set by fast cooling near the GC where the particles all had the same initial conditions. This is the reason why $E_{\rm max}$ is nearly spatially uniform at the present day (top right panel of Figure 1).

Figure 4 shows the spectra of the CRe and their evolution after line-of-sight projections. The red, green, and blue curves represent the CR spectra projected into a longitude range of $l=[-10^{\circ},10^{\circ}]$ and latitude ranges of $[40^{\circ},60^{\circ}],[20^{\circ},40^{\circ}]$, and $[10^{\circ},20^{\circ}]$, respectively. The top left panel shows the spectra for different latitudes at the present day. This plot confirms our expectation that the maximum energy of the CRe, $E_{\rm max}$, is only mildly varying with latitudes, ranging from $\gtrsim 100$ GeV at low latitudes to $\sim 1$ TeV at higher latitudes (see also top right panel of Figure 1). Again, the spatial uniformity of $E_{\rm max}$ is resulted from initial fast cooling and subsequent mild adiabatic losses. This process can be seen from the other three panels of Figure 4, in which we plot the CR spectral evolution for the three latitude bins. At early times when the jets just shut off ($t=0.3$ Myr), only the low latitude bin is populated with CRs (see the long-dashed curve in the lower right panel). Due to the initial cooling, the spectra at this time already showed an exponential cutoff at $\sim 1$ TeV. Afterwards, the CRs are propagated to higher latitudes, but $E_{\rm max}$ only mildly shifts to lower energies due to adiabatic cooling (the shift is strongest for the lowest latitude bin because the CRe also suffer from stronger synchrotron and IC losses). In terms of the amplitudes of the spectra, in general they decrease with time owing to cooling; this trend is only inverted when the CRs first entered the lowest and highest latitude bins.

Because of the spatially uniform distribution of $E_{\rm max}$, we expect the gamma-ray spectrum to have a spatially uniform high-energy cutoff at similar energies ($E_{\rm max,obs}\lesssim E_{\rm max}$ for intermediate and high latitudes; see § 2). In the upper left panel of Figure 5 we plot the simulated gamma-ray spectra of the Fermi bubbles calculated for a longitude range of $l=[-10^{\circ},10^{\circ}]$ for different latitude bins. For a given longitude and latitude range, the simulated spectrum is computed by projecting the gamma-ray emissivities as a function of energy along finely sampled lines of sight (with resolutions of 0.5 degrees), and then we average the spectra over all the sightlines within the region. Indeed, we find that the simulated spectra for all latitudes exhibit a spectral cutoff at similar energies around several hundreds GeV, consistent with the observed cutoff energy.

The top left panel of Figure 5 also shows that, not only is the high-energy cutoff similar across all latitudes, but the general shape of the simulated spectra are also latitude independent as observed. As discussed in § 2, the energy of CRe must be slightly higher at higher latitudes because the CMB (optical) photons dominates the ISRF at high (low) latitudes. Here we note that although we did not account for diffusion, CR diffusion in the solar neighborhood is known to be energy dependent, scaling as $E^{0.3-0.6}$ (Strong et al., 2007). Advection expands the bubbles to $\sim$ 6 kpc in $\sim$ 1 Myr; for CRs to diffuse a comparable distance their diffusivity $D$ would have to be of order $1.1\times 10^{31}$ cm2s*-1*. This is 220 times larger than $D$ for GeV particles ($D\sim 5\times 10^{28}$ cm2s*-1*), but if $D$ scales with energy as $E^{0.5}$, diffusion would dominate advection for $E>50$ TeV. While diffusion is likely to be only a small effect, it goes in the right direction to explain an excess of high energy particles at large heights. In order to see this effect more clearly, we decompose the spectra for each latitude bin into three components, which are calculated from the CMB, IR, and optical photons in the ISRF. At high latitudes ($b=40^{\circ}-60^{\circ}$), the simulated CRe have an average energy of $\langle E_{\rm e}\rangle\sim 530$ GeV, and therefore the gamma-ray spectrum has a bump with $E_{\rm bump}\sim 1$ GeV after up-scattering the CMB photons (see Eq. 5). Because for this latitude bin the IC scattering goes from the Thomson regime to the KN regime, one can also see the change in spectral indices from $\sim 1.5$ to $\sim 3$ from low to high gamma-ray energies (see § 2). For intermediate latitudes ($b=20^{\circ}-40^{\circ}$), the three components make comparable contributions to the spectrum. For the low latitude bin ($b=10^{\circ}-20^{\circ}$), the gamma-ray emission is dominated by IC scattering of the optical starlight, with $E_{\rm bump}\sim 40$ GeV and an average CR energy of $\langle E_{\rm e}\rangle\sim 40$ GeV. Because at low latitudes the scattering is in the Thomson limit, the spectral index is $\sim 1.5$ up to the cutoff energy.

In short summary, we demonstrated that the spectra of the Fermi bubbles are nearly latitude independent because the CRe from different parts of the bubbles at the present day all originate from the GC where they suffer from fast synchrotron and IC cooling soon after injections. We reproduced the latitude-independent cutoff energy and spectral shape of the gamma-ray spectra despite the complex convolution of CR energies and the latitude-dependent ISRF. We also note that the normalizations of the simulated spectra for different latitude bins are comparable to one another (top left panel in Figure 5), indicating the flat surface brightness of the observed bubbles is also recovered. This is quite a remarkable result since one must get the CR distribution right both spatially and spectrally in order to successfully reproduce the flat intensity and latitude-independent spectra simultaneously.

4.4. Constraints on the initial conditions

Because the maximum energy of the CRe at the present day, $E_{\rm max}$, is largely determined by fast cooling of CRe near the GC, it could be used to constrain the initial conditions at injection, including the initial speed of the AGN jets and the energy densities of the ISRF and the magnetic field. In this section we discuss the parameter space allowed to build a successful model, and how it would be influenced by improved measurements of the cutoff energy from future observational data. In deriving these constraints, we assume that no significant re-acceleration of CRs took place near the GC.

Two criteria need to be satisfied at early stage of the bubble evolution in order to generate a spatially uniform bubble spectrum in the scenario described in § 4.3. First, the initial cooling must be fast enough to act on the jets before they propagate away from the GC. Therefore, the cooling timescale of CRe must be shorter than the dynamical time of the jets, i.e., $\tau_{\rm syn+IC}<\tau_{\rm dyn}$. Using the expression for the synchrotron and IC cooling time (Eq. A28) and the definition of $t_{\rm dyn}\equiv(1\ {\rm kpc})/v_{\rm jet}$, we obtain an upper limit on the initial jet velocity,

[TABLE]

where $c$ is the speed of light, $E_{\rm max,0}$ is the characteristic maximum energy of CRe near the GC, $u_{\rm tot}=u_{\rm B}+u_{\rm rad}F_{\rm KN}$ is the summation of the energy density of the magnetic field and the ISRF with the correction factor for the KN effect (Moderski et al., 2005). Note that the strengths for both the magnetic field and the ISRF rapidly decay away from the GC, and hence $u_{\rm tot}$ in the above equation represents an average value near the GC (roughly within the central kpc). For the following discussion, we assume $f_{\rm cool}\equiv E_{\rm max}/E_{\rm max,0}=0.3$ to account for the difference between the characteristic CR energy near the GC ($E_{\rm max,0}$) and that observed today ($E_{\rm max}$).

Another criterion comes from the fact that the initial cooling cannot be so strong that the energy of the CRe cools below the energy required to produce the observed high-energy cutoff today. In other words, the energy of CRe after the initial cooling losses has to be greater than the maximum energy of the CRe today, i.e., $E>E_{\rm max}$. The CR energy after going through synchrotron and IC losses is given by $E=E_{0}/(1+\beta tE_{0})$ (Kardashev, 1962), where $E_{0}$ is the initial CR energy and $\beta=(4/3)(\sigma_{T}/m_{e}^{2}c^{3})u_{\rm tot}$. For very large $E_{0}$, the CR energy after cooling is approximately

[TABLE]

The requirement of $E>E_{\rm max}$ gives a lower limit on the jet speed,

[TABLE]

In Figure 6 we plot the permitted values of $v_{\rm jet}$ as a function of $u_{\rm tot}$ bracketed by the above two criteria (Eq. 7 and 9) in the shaded region, assuming $E_{\rm max}=300$ GeV. The color shows the value of $E_{\rm max}$ for given values of $v_{\rm jet}$ and $u_{\rm tot}$ (Eq. 8). The parameter set adopted in the current simulation (plotted using the star symbol) lies within the permitted parameter space and is therefore able to successfully reproduce the spatially uniform spectrum of the bubbles. However, this figure illustrates that the solution is not unique. 444Though not unique, the jet parameters adopted in the current simulations have been shown to satisfy many other observational constraints (see Y12 for detailed discussion), in additional to those presented here. For example, for the current observational constraint of $E_{\rm max}\gtrsim E_{\rm max,obs}\sim 300$ GeV (near the lower solid line), if we were to use an average energy density for the magnetic field and the ISRF of $u_{\rm tot}=2\times 10^{-11}\ {\rm erg\ cm^{-3}}$, the initial velocity of the AGN jets would have to be in the range of $0.04c-0.13c$ in order to have a successful model. Generally speaking, in order to produce CRe with energy $E_{\rm max}\gtrsim 300$ GeV, the initial jet velocity must be faster (slower) for larger (smaller) initial strength for the magnetic field and ISRF.

Figure 6 also shows that, assuming $u_{\rm tot}$ at the time of injection is not significantly smaller than the value adopted in the current simulation, the required initial velocity of the outflow that transports the CRe must be at least $\sim 0.01c$ or $3000\ {\rm km\ s^{-1}}$. Such a fast speed is easily achievable by AGN jets but not by winds driven by nuclear starburst, for example. Therefore, the mechanism for generating spatially uniform spectrum proposed in this work would not be applicable for models that are based on outflows with lower velocities.

Finally, we discuss the influence on the allowable parameter space by improved constraints on $E_{\rm max}$ in the future from GeV and TeV observations such as Fermi, High Altitude Water Cherenkov (HAWC) (Abeysekara et al., 2017), Cherenkov Telescope Array (CTA), Large High Altitude Air Shower Observatory (LHAASO), and the Hundred Square km Cosmic ORigin Explorer (HiScore). In Figure 6 we plot the permitted parameter space assuming $E_{\rm max}=3$ TeV and 30 TeV (bracketed by the dashed and dotted lines, respectively), assuming a constant $f_{\rm cool}=0.3$. If the constraints from future data finds $E_{\rm max}$ to be greater than 300 GeV, the average energy density of the magnetic field and the ISRF near the GC, $u_{\rm tot}$, would have to be smaller, and/or the initial jet velocity must be higher, in order to be consistent with the observed $E_{\rm max}$. As the $E_{\rm max}$ gets bigger and bigger, the limits on $v_{\rm jet}$ and $u_{\rm tot}$ become more and more stringent. Therefore, if $E_{\rm max}$ approaches tens of TeV, the leptonic jet model would become less favorable because it is difficult for the jets to avoid cooling and keep the CRe at such high energies, unless significant re-acceleration of CRs occured near the GC to compensate for the cooling losses. If there were CR re-acceleration, it would have similar effects as slower cooling and shift the permitted parameter space shown in Figure 6 downward.

5. Conclusions

One of the unique features of the Fermi bubbles is the spatially uniform spectrum, including the spectral shape and the high-energy cutoff above 110 GeV. Because reproducing the latitude-independent spectrum requires the correct CR distribution both spatially and spectrally, it provides stringent constraints on the theoretical models proposed to explain the origin of the bubbles. In this work, we investigate the spectral evolution of the Fermi bubbles in the leptonic AGN jet scenario using 3D hydrodynamic simulations that include modeling of the CR spectrum. The simulations are done using the newly implemented CRSPEC module in FLASH, which allows us to track the spectral evolution of CRe due to adiabatic processes and synchrotron plus IC cooling after they are injected with the AGN jets from the GC. Our main findings are summarized as follows.

(1) The high-energy cutoff in the observed gamma-ray spectrum of the bubbles is a signature of fast synchrotron and IC cooling of CRe near the GC when the jets were first launched.

(2) After the initial phase of fast cooling near the GC, the dynamical time of the jets becomes the shortest among all other cooling timescales and therefore the CR spectrum is essentially advected with only mild cooling losses. This could explain why the bubble spectrum is nearly spatially uniform, because the CRe from different parts of the bubbles today all share the same origin near the GC at early stage of the bubble expansion.

(3) The simulated distribution of CR energies, despite being quite uniform, has a slight gradient toward higher energies at higher Galactic latitudes. We show that this is essential for reproducing the latitude-independent shape of the gamma-ray spectrum because the ISRF is dominated by lower-energy CMB photons at high latitudes and optical starlight at low latitudes.

(4) Because the observed cutoff energy of the gamma-ray spectrum today is closely related to the early phase of fast cooling, it can be used to constrain the initial conditions near the GC, such as the initial speed of AGN jets and the energy density of the magnetic field and the ISRF. The permitted parameter space for building a successful model and its dependence on the future measurements of the cutoff energy are summarized in Figure 6.

Finally, we note that in addition to the above spectral features, the simulated 3D CR distribution is edge-brightened (Figure 1), which is key for recovering the flat surface brightness of the observed bubbles, or the latitude-independent normalization of the observed spectrum (top left panel of Figure 5). It is remarkable that the leptonic jet model predicts the right spatial and spectral distribution of CRe that simultaneously matches the normalization, overall spectral shape, and high-energy cutoff of the observed gamma-ray spectrum and their spatial uniformity. Together with the fact that the microwave haze is more easily explained by the leptonic jet model, we conclude that past AGN jet activity is a likely mechanism for the formation of the Fermi bubbles. Future data from multi-messenger observations, particularly improved measurements of the cutoff energy of the gamma-ray spectrum by GeV and TeV observatories including Fermi, HAWC, CTA, LHAASO, and HiScore, will provide crucial verification of the scenario proposed in this work.

The authors would like to thank Jay Gallagher for helpful discussion. We also thank Ellen Zweibel and an anonymous referee for useful comments that helped to improve the manuscript. HYKY acknowledges support by NASA through Einstein Postdoctoral Fellowship grant number PF4-150129 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. MR acknowledges support from NASA grant NASA ATP 12-ATP12-0017 and NSF grant AST 1715140. The simulations presented in this paper were performed on the Deepthought2 cluster, maintained and supported by the Division of Information Technology at the University of Maryland College Park. FLASH was developed largely by the DOE-supported ASC/Alliances Center for Astrophysical Thermonuclear Flashes at the University of Chicago. Data analysis presented in this paper was partly performed with the publicly available yt visualization software (Turk et al., 2011).

Appendix A The CRSPEC module

The evolution of CR particles is described by the diffusion-advection equation (Skilling, 1975),

[TABLE]

where $f({\bm{x}},p)$ is the isotropic part of the particle distribution function, $\kappa(p)$ and $D_{pp}(p)$ are the diffusion coefficients in spatial coordinates and in momentum space, respectively, $b_{l}(p)$ describes mechanical and radiative losses (see Eq. A29), and $j({\bm{x}},p)$ is the source term accounting for CR injections at shocks or production of secondary particles. In order to save computational costs, the momentum space is divided into relatively sparse $N_{p}$ logarithmically spaced momentum bins, bounded by $p_{0}$,…,$p_{N_{p}}$. The width of the bins on a log scale, $\Delta\log p=\log(p_{i}/p_{i-1})$, is assumed to be a constant for convenience. The modeled CRs have a finite spectrum with minimum and maximum momenta of $p_{cutL}$ and $p_{cutR}$, respectively, and one has to ensure that $p_{0}$ and $p_{N_{p}}$ bracket the evolution of $p_{cutL}$ and $p_{cutR}$. For momentum bins that at least partially intersect with the CR spectrum, the distribution function $f(p)$ is approximated by a piece-wise power-law distribution,

[TABLE]

where $f_{i}$ and $q_{i}$ are the normalization and logarithmic slope for the $i$th bin. For momentum bins that do not contain any CRs, i.e., if $p_{i-1}>p_{cutR}$ or $p_{i}<p_{cutL}$, $f_{i}$ is assigned to be zero. After integrating over the $i$th momentum bin of Eq. A1 multiplied by $4\pi p^{2}$, we obtain equations for the CR number densities,

[TABLE]

The advection and spatial diffusion terms (i.e., first and second terms in Eq. A3) are implemented in the same way as in the energy-integrated version of the CR module (Y12). The terms that represent second order Fermi acceleration ($\propto D_{pp}(p)$) and the source term ($Q_{i}$) are neglected hereafter in order to focus on processes relevant for this paper. The evolution of $n_{i}$ due to advection in momentum space becomes

[TABLE]

where $b(p)$ includes adiabatic compression or expansion as well as other energy loss terms $b_{l}(p)$ (see Eq. A29). Note that $b(p)$ is greater (smaller) than zero when CRs are cooling (heating). Integration of Eq. A7 over time gives

[TABLE]

where $\Phi_{i}^{p}$ is the time-averaged flux evaluated at the cell boundary $p_{i}$. Using Eq. A8 to rewrite the above equation, we have

[TABLE]

where $p_{u}$ is the upstream momentum, which can be solved using the equation

[TABLE]

Note that because of the finite CR spectrum, when CRs are cooling ($b(p)>0$), Eq. A11 is integrated up to $\min(p_{u},p_{cutR})$; when CRs are heating ($b(p)<0$), the upper integration limit is $\max(p_{u},p_{cutL})$.

For CR ions, at each timestep one can solve for $f_{i}$ and $q_{i}$ given the updated $n_{i}$ based on Eq. A7 assuming that the curvature of the spectrum is constant (see Miniati et al. (2001) for detailed discussion). However, for fast cooling CRe, this assumption is not valid and therefore it is necessary to employ other constraints. This can be achieved by evolving the CR energy densities, $e_{i}$, in addition to the number densities, $n_{i}$. The equations for $e_{i}$ can be derived by taking one moment of Eq. A1. That is, we multiply Eq. A1 by $4\pi p^{2}T(p)$, where $T(p)=(\gamma-1)m_{e}c^{2}$ is the particle kinetic energy and $\gamma$ is the Lorentz factor, and integrate over the $i$th momentum bin. This yields

[TABLE]

where Eq. A2 is used and sub-relativistic contribution is ignored in order to derive the last expression in Eq. A15. Same as before, we neglect the advection, diffusion, and source terms and only focus on the energy transfer in momentum space, i.e.,

[TABLE]

Using Eq. A2, the second term in the above equation can be rewritten as $e_{i}R_{i}$, where

[TABLE]

Integrating Eq. A18 over time gives

[TABLE]

where $\Phi_{i}^{e}$ is the time-averaged flux evaluated at the cell boundary $p_{i}$. Again, we could rewrite the equation using Eq. A8 and obtain

[TABLE]

where $j$ and $p_{u}$ are defined in Eq. A12 and A13, respectively. Similar to Eq. A11, the upper integration limit is $\min(p_{u},p_{cutR})$ and $\max(p_{u},p_{cutL})$ when CRs are cooling and heating, respectively. The minimum and maximum momenta of the CR spectrum ($p_{cutL}$ and $p_{cutR}$, respectively) are solved explicitly using

[TABLE]

In the simulations, in addition to hydrodynamic variables, extra $2N_{p}+2$ variables are stored for $n_{i}$, $e_{i}$, $p_{cutL}$, and $p_{cutR}$. At each simulation timestep, after accounting for CR advection and diffusion, we first convert $n_{i}$ and $e_{i}$ into $f_{i}$ and $q_{i}$, where $q_{i}$ could be solved using the following equation for each momentum bin (assuming $p_{i}\gg m_{e}c$):

[TABLE]

and $f_{i}$ could be computed directly using Eq. A4 or A15. We then update $n_{i}$ and $e_{i}$ from $t$ to $t+\Delta t$ using Eq. A9 and A20. Finally the minimum and maximum momenta of the CR spectrum are updated using Eq. A23.

In general, in order to achieve numerical accuracy, the simulation timestep $\Delta t$ has to satisfy $|\log\frac{p_{u}}{p_{i}}|\leq\epsilon\log\frac{p_{i}}{p_{i-1}}$, where $\epsilon\leq 1$ is similar to the Courant number. However, when fast cooling CRe are included, $\Delta t\lesssim 0.1\ \tau_{cool}$ should be adopted in order to accurately follow the high-energy end of the spectrum. In the current implementation in FLASH, when timestep constraints due to CR cooling become the most limiting among all relevant timesteps, we subcycle over the CR spectral evolution in order to speed up the computation. Because of the block-structured nature of the AMR architecture in FLASH, neighboring blocks are likely assigned to the same processor. This causes significant load imbalance among processors when a particular region in the simulation domain (e.g., near the GC for the current simulation) suffers from fast CRe cooling. To this end, we pay special attention to load balancing in order to achieve good parallel performance.

For simulations presented in this paper, only synchrotron and IC losses for CRe are relevant. We refer the readers to Strong & Moskalenko (1998) for the expressions for other energy losses of CRe and CRp. The energy loss rate of CRe due to synchrotron and IC losses is given by

[TABLE]

where $\beta=\sqrt{1-v^{2}/c^{2}}$, $\gamma=E/(m_{e}c^{2})$ and $u_{B}$ and $u_{rad}$ are energy densities in magnetic field and the radiation field in units of ${\rm erg\ cm^{-3}}$, respectively. The factor $F_{KN}$ accounts for the reduced IC cross section in the KN regime:

[TABLE]

where

[TABLE]

is an analytical approximation for the general KN formula (Moderski et al., 2005). The synchrotron and IC cooling time ($\equiv E/(dE/dt)$) is

[TABLE]

By defining $\hat{p}\equiv p/m_{e}c$, one can write the energy of CR particles as $E=\sqrt{\hat{p}^{2}+1}m_{e}c^{2}$. Therefore, the momentum loss rate as used in Eq. A8 is related to the energy loss rate by

[TABLE]

Figure 7 shows a test of the CRSPEC module including synchrotron losses of CRe. The initial CR spectrum ranges from $10^{2}$ to $10^{6}$ GeV with constant spectral indices of $q_{i}=5$. The initial spectrum is normalized such that the number density of the first momentum bin $n_{1}=10^{3}\ {\rm cm^{-3}}$. The results are in excellent agreement with the analytical solution (Kardashev, 1962).