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

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.


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 Gammaray Space Telescope (Su et al. 2010;Ackermann et al. 2014;Narayanan & Slatyer 2017). The observed gammaray 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), Xray 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., Yang et al. 2012Yang et al. , 2013, the latter two are abbreviated as Y12 and Y13, respectively), and in-situ acceleration models (e.g., Mertsch & Sarkar 2011;Cheng et al. 2011Cheng et al. , 2015Sarkar 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 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 ac-celeration 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 ∼ 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 cosmicray (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 ∼ 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.

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 cut ), where E cut ∼ 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 max,obs , which relates to E cut by exp(−E max,obs /E cut ) = 0.1 (i.e., the energy scale where the gamma-ray intensity is dimmer by a factor of 10). Given E cut = 110 GeV, we have E max,obs ∼ 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 γ is given by Blumenthal & Gould (1970), where E ph is the initial photon energy, γ = E e /(m e c 2 ) is the Lorentz factor of the CR electron, E γ is the energy of the up-scattered gamma-ray photon, and n(E ph ) is the energy distribution of the photon number density. In the Thomson limit (Γ 1), the average energy of the up-scattered photons is given by In the Klein-Nishina (KN) limit (Γ 1), almost all the energy of the CR electron is carried away by the up-scattered photons, i.e., E γ ∼ E e . Assuming the CR spectrum is a power law with spectral index α, it can be shown that the spectral index of the up-scattered gamma-ray photons is (α + 1)/2 in the Thomson limit, and approximately α + 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 ∼ 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 bump ∼ 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 Γ < 1 for E γ ∼ 10 GeV and E ph < 10 eV.
For high latitudes (e.g., b = 40 • −60 • ), the ISRF peaks at E ph ∼ 7×10 −4 eV for CMB photons. It would therefore require an average CR energy of E e ∼ 2 TeV in order to produce a ∼ 10 GeV bump. For intermediate latitudes (b = 20 • − 40 • ), the intensity of the ISRF is more uniform across all wavelengths. Assuming the gammaray bump primarily comes from infrared (IR) photons (E ph ∼ 10 −2 eV), one would obtain E e ∼ 200 GeV. Similarly, the average CR energy can be estimated to be E e ∼ 20 GeV at low latitudes where optical light (assuming E ph ∼ 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 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 max,obs , comes from upscattered optical light in the ISRF and provides information about the maximum energy of the underlying CR electron population, E max . For optical photons ( E ph ∼ 5 eV), the IC scattering would be in the KN limit for CRe with energies greater than ∼ 100 GeV. Therefore, for high and intermediate latitudes, E max E max,obs ∼ 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 max,obs = 250 GeV, Eq. 5 gives E max ∼ 100 GeV. The observed cutoff energy is almost independent of latitudes, with a slight tendency of higher E max,obs for higher latitudes ( Figure  33 of Ackermann et al. (2014)). Therefore, the estimates above imply that E 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 gammaray spectrum obtained by Narayanan & Slatyer (2017).

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 adaptivemesh-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 × 10 −3 of the injected CRs to be CRe 4 and follow their spectral evolution due to adiabatic compression, adiabatic expansion, synchrotron losses, and IC losses. The rest of the injected CR energy density 5 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 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 µ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 GAL-PROP 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 ∼ 19 eV cm −3 near the GC to ∼ 1 eV cm −3 at 5 kpc away from the Galactic plane near the rotational axis.

Modeling CR spectral evolution
The core of this work is the newly implemented CR-SPEC 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, where f i and q i are the normalization and logarithmic slope for the ith 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 × 10 −3 of the total CR energy density of the jets, or equivalently, 7.5 × 10 −12 erg cm −3 . A representative value for the normalization factor at 10 GeV is 2.24 × 10 −4 cm −2 s −1 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. 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 6 , 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.

Distribution of CR energies
Because the CRe only contribute to 3×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-ofsight 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 cr ∼ E e n cr , where E e is the characteristic energy of the CRe, this implies that E e is largely spatially uniform. This is indeed suggested by the map of maximum CR electron energy E 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 max is relatively small, with E max on the order of a few hundreds of GeV for most regions within the bubbles.
Although E e and E max are generally quite spatially uniform, they do exhibit some gradients. The top right panel in Figure 1 shows that E max varies from ∼ 100 GeV at low Galactic latitudes to ∼ 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. How-6 The maximum CR electron energy, Emax, 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 Emax for regions where the total energy density of CRe is greater than 10 −16 erg cm −3 , which corresponds to a minuscule level of ∼ 6 × 10 −5 eV cm −3 . ever, 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. 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 max ), minimum spectral energy (E 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 max and E 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 ∼ 0.1 Myr), there was a brief increase in E min and E max due to adiabatic compression, i.e., ∇ · v < 0. Afterwards, the CRe propagate outward and the only cooling mechanism is adiabatic expansion (∇ · v > 0), and therefore both E min and E max decrease with time monotonically.

Spectral evolution of the Fermi bubbles
In contrast, the evolution of E 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 ∼ 1 TeV before t ∼ 0.3 Myr. This fast change in E 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 rel-evant timescales (see the third panel). After t ∼ 0.4 Myr, the dynamical timescale (τ dyn ≡ (1 kpc)/v) becomes shorter/comparable to the synchrotron and IC cooling timescale, while the timescale for adiabatic processes (τ AD ≡ 1/(∇ · 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 max only decreases slightly after t ∼ 0.4 Myr and reaches a value about 700 GeV at t = 1.2 Myr. In other words, the value of E 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.

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 max,obs ∼ 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 max is very similar, marked with a fast decay within the early ∼ 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 (τ syn+IC τ 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 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 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 • , 10 • ] and latitude ranges of [40 • , 60 • ], [20 • , 40 • ], and [10 • , 20 • ], 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 max , is only mildly varying with latitudes, ranging from 100 GeV at low latitudes to ∼ 1 TeV at higher latitudes (see also top right panel of Figure 1). Again, the spatial uniformity of E 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 ∼ 1 TeV. Afterwards, the CRs are propagated to higher latitudes, but E 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 max , we expect the gamma-ray spectrum to have a spatially uniform high-energy cutoff at similar energies (E max,obs E 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 • , 10 • ] 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 ∼ 6 kpc in ∼ 1 Myr; for CRs to diffuse a comparable distance their diffusivity D would have to be of order 1.1 × 10 31 cm 2 s −1 . This is 220 times larger than D for GeV particles (D ∼ 5 × 10 28 cm 2 s −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 • − 60 • ), the simulated CRe have an average energy of E e ∼ 530 GeV, and therefore the gamma-ray spectrum has a bump with E bump ∼ 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 ∼ 1.5 to ∼ 3 from low to high gamma-ray energies (see § 2). For intermediate latitudes (b = 20 • − 40 • ), the three components make comparable contributions to the spectrum. For the low latitude bin (b = 10 • − 20 • ), the gammaray emission is dominated by IC scattering of the optical starlight, with E bump ∼ 40 GeV and an average CR energy of E e ∼ 40 GeV. Because at low latitudes the scattering is in the Thomson limit, the spectral index is ∼ 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 latitudedependent 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.

Constraints on the initial conditions
Because the maximum energy of the CRe at the present day, E 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 . The other three panels show decomposition of the simulated spectra into different components of the ISRF, namely the CMB (dashed-triple-dotted), IR (dashed), and optical (dotted) radiation field. The grey band represents the observational data of Ackermann et al. (2014). The leptonic jet model successfully reproduced the latitude independence of the observed spectra, including the normalization, overall spectral shape, and the spectral cutoff above ∼ 110 GeV, despite the complex convolution of CR energies and the latitude-dependent ISRF.
it would be influenced by improved measurements of the cutoff energy from future observational data. In deriving these constraints, we assume that no significant reacceleration 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., τ syn+IC < τ dyn . Using the expression for the synchrotron and IC cooling time (Eq. A28) and the definition of t dyn ≡ (1 kpc)/v jet , we obtain an upper limit on the initial jet velocity, v jet < 0.065c u tot 10 −11 erg cm −3 where c is the speed of light, E max,0 is the characteristic maximum energy of CRe near the GC, u tot = u B + u rad F 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 tot in the above equation represents an average value near the GC (roughly within the central kpc). For the follow-ing discussion, we assume f cool ≡ E max /E max,0 = 0.3 to account for the difference between the characteristic CR energy near the GC (E max,0 ) and that observed today (E 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 max . The CR energy after going through synchrotron and IC losses is given by E = E 0 /(1 + βtE 0 ) (Kardashev 1962), where E 0 is the initial CR energy and β = (4/3)(σ T /m 2 e c 3 )u tot . For very large E 0 , the CR energy after cooling is approximately The requirement of E > E max gives a lower limit on the jet speed, v jet > 0.02c u tot 10 −11 erg cm −3 E max 300 GeV .
In Figure 6 we plot the permitted values of v jet as a function of u tot bracketed by the above two criteria (Eq. 7 and 9) in the shaded region, assuming E max = 300 GeV. The color shows the value of E max for given values of v jet and u 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. 7 For example, for the current observational constraint of E max E max,obs ∼ 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 tot = 2 × 10 −11 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 max 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 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 ∼ 0.01c or 3000 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 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 max = 3 TeV and 30 TeV (bracketed by the dashed and dotted lines, respectively), assuming a constant f cool = 0.3. If the constraints from future data finds E max to be greater than 300 GeV, the average energy density of the magnetic field and the ISRF near the GC, u tot , would have to be smaller, and/or the initial jet velocity must be higher, in order to be consistent with the observed E max . As the E max gets bigger and bigger, the limits on v jet and u tot become more and more stringent. Therefore, if E 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.

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, 7 Though 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. -Allowable parameter space for successful models. The x axis represents an average value of the summation of energy densities from the ISRF and magnetic field near the GC, and the y axis is the initial velocity of the jets. The color shows log(Emax/GeV), where Emax is the value of maximum CR energy, for given values of v jet and utot. Parameters within the shaded region satisfy the upper and lower limits of v jet given by Eq. 7 and Eq. 9, respectively, assuming Emax = 300 GeV and f cool = 0.3 (see the text for definition). The star symbol shows the parameters used in the current simulation. The region bracketed by the dashed and dotted lines is the permitted parameter space assuming Emax = 3 TeV and 10 TeV, respectively, indicating that future observational limits of Emax, if bigger than 300 GeV, would shift the allowable parameter space to the upper-left corner. These constraints are derived assuming there is no significant re-acceleration of CRs near the GC.
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 lati-tudes and optical starlight at low latitudes.
(4) Because the observed cutoff energy of the gammaray 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 edgebrightened (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). where p u is the upstream momentum, which can be solved using the equation ∆t = pu pi dp b(p) . (A13) 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πp 2 T (p), where T (p) = (γ − 1)m e c 2 is the particle kinetic energy and γ is the Lorentz factor, and integrate over the ith momentum bin. This yields 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., m 2 e c 2 + p 2 dp.
Using Eq. A2, the second term in the above equation can be rewritten as e i R i , where p 3−qi m 2 e c 2 + p 2 dp p R p L p 2−qi T (p)dp.
Integrating Eq. A18 over time gives where Φ e i is the time-averaged flux evaluated at the cell boundary p i . Again, we could rewrite the equation using Eq. A8 and obtain Φ e i = − 4π ∆t pu pi p 2 f j (p)T (p)dp, 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 ∆t = p t cut p t+∆t cut dp b(p) .
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 m e c): and f i could be computed directly using Eq. A4 or A15. We then update n i and e i from t to t + ∆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 ∆t has to satisfy | log pu pi | ≤ log pi pi−1 , where ≤ 1 is similar to the Courant number. However, when fast cooling CRe are included, ∆t 0.1 τ 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 where β = 1 − v 2 /c 2 , γ = E/(m e c 2 ) and u B and u rad are energy densities in magnetic field and the radiation field in units of erg cm −3 , respectively. The factor F KN accounts for the reduced IC cross section in the KN regime: where f KN (x) 1 (1 + Γ) 1.5 for Γ = 4γE ph m e c 2 10 4 (A27) is an analytical approximation for the general KN formula (Moderski et al. 2005 By definingp ≡ p/m e c, one can write the energy of CR particles as E = p 2 + 1m e c 2 . Therefore, the momentum loss rate as used in Eq. A8 is related to the energy loss rate by b l (p) ≡ dp dt = √p + 1 cp dE dt .
(A29) 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 cm −3 . The results are in excellent agreement with the analytical solution (Kardashev 1962).