Sources and Radiations of the Fermi Bubbles

Abstract

Two enigmatic gamma-ray features in the galactic central region, known as Fermi Bubbles (FBs), were found from Fermi-LAT data. An energy release, (e.g., by tidal disruption events in the Galactic Center, GC), generates a cavity with a shock that expands into the local ambient medium of the galactic halo. A decade or so ago, a phenomenological model of the FBs was suggested as a result of routine star disruptions by the supermassive black hole in the GC which might provide enough energy for large-scale structures, like the FBs. In 2020, analytical and numerical models of the FBs as a process of routine tidal disruption of stars near the GC were developed; these disruption events can provide enough cumulative energy to form and maintain large-scale structures like the FBs. The disruption events are expected to be ${10}^{-4}\sim {10}^{-5}{\mathrm{yr}}^{-1}$, providing an average power of energy release from the GC into the halo of $\dot{\mathcal{E}}\sim 3\times {10}^{41}$ erg ${\mathrm{s}}^{-1}$, which is needed to support the FBs. Analysis of the evolution of superbubbles in exponentially stratified disks concluded that the FB envelope would be destroyed by the Rayleigh–Taylor (RT) instabilities at late stages. The shell is composed of swept-up gas of the bubble, whose thickness is much thinner in comparison to the size of the envelope. We assume that hydrodynamic turbulence is excited in the FB envelope by the RT instability. In this case, the universal energy spectrum of turbulence may be developed in the inertial range of wavenumbers of fluctuations (the Kolmogorov–Obukhov spectrum). From our model we suppose the power of the FBs is transformed partly into the energy of hydrodynamic turbulence in the envelope. If so, hydrodynamic turbulence may generate MHD fluctuations, which accelerate cosmic rays there and generate gamma-ray and radio emission from the FBs. We hope that this model may interpret the observed nonthermal emission from the bubbles.

1. Introduction: Sources of the Fermi Bubbles

In this article, we present our interpretation of the origin of the Fermi Bubbles (FBs). The discussion includes the energy release in the Galactic Center (GC) to the hydrodynamic envelope in the galactic halo, the excitation of MHD turbulence that accelerates cosmic rays (CRs) in the halo, the processes of nonthermal emissions which are observed in X-ray, gamma-ray, and radio ranges, and the high-energy CRs escaping from the FBs to the galactic disk. We present a puzzle of the FB picture, where many fragments are still missing in the mosaic. Many of them are interpreted but not completely understood yet. The goal of the article is to find a way to find a proper solution to this problem.

The origin of the energy release in the FBs in the GC is still an open question. This kiloparsec-scale structure was interpreted as a manifestation of past activity of the central supermassive black hole (SMBH) Sgr A* in the GC (see Figure 1). Observations in the GC showed structures above and below it of gamma-rays, microwaves, and X-rays. The eROSITA [1] found giant bubbles in the X-ray range 0.1∼2.4 keV, extending approximately 14 kiloparsecs above and below the GC. The estimated energy of the bubbles is around ${10}^{56}$ erg. The total luminosity in X-rays is about ${10}^{39}$ erg ${\mathrm{s}}^{-1}$, which could be the result of past activity in the GC. The temperature of the envelope is about $0.3$ keV, the velocity of the shock is about 340 km ${\mathrm{s}}^{-1}$ or of Mach number $\approx 1.5$, and the energy-release rate of the gas envelope is roughly ${10}^{41}$∼${10}^{42}$ erg ${\mathrm{s}}^{-1}$.

The inner radius of the X-ray shell (about 7 kpc) coincides spatially with the region of GeV gamma rays in the range 1∼100 GeV with the luminosity of $F_{\gamma }\approx 4\times {10}^{37}$ erg $s^{-1}$ [2]. The bubble structure in the GC was also revealed in the range of microwaves which coincides nicely with that of gamma rays [3]. The flux is in the range 23∼61 GHz, and the luminosity is ${\Phi }_{\nu }\approx 1$∼$5\times {10}^{36}$ erg ${\mathrm{s}}^{-1}$.

Similar giant structures near the GC were found earlier in the radio range of hundreds of MHz (the North Polar Spur; see [4]) and in 1.5 keV X-ray emission (see [5]). These structures were postulated as bipolar supershells which were produced by starbursts. A shock front was supposed to reach a radius of 10 kpc in the polar regions, which could be consistent with the GC explosions. This model required an energy release of about ${10}^{55}$ erg at the GC, and periodic activity on a time scale of 10∼15 Myr.

Mou et al. [6] suggested that the nature of the North Polar Spur in the GC agreed with the eROSITA bubble of an age of about 20 Myr [1].

Bubbles have also been discovered in other galaxies (see review by [7]). Figure 2 shows superbubbles from the galaxy NGC 3079.

The total energy needed to generate large galactic outflows is assumed to be in the range of up to about ${10}^{56}$ erg. This energy release in the GC may be compelling evidence for a huge energetic explosion having occurred in the GC a few (2∼8) million years ago (see, e.g., [5] and references therein). For examples, Nayakshin & Zubovas [8] assumed a capture of a giant molecular cloud of mass ∼${10}^5{M}_{\odot }$ in the GC about one Myr ago; and Yang et al. [9], Yang et al. [10] suggested a model of FBs as a result of past activity at the GC.

Alternative models of the bubbles were suggested by, e.g., Cheng et al. [11], Cheng et al. [12], Cheng et al. [13], Zubovas & Nayakshin [14], Mertsch & Sarkar [15], Ko et al. [16], etc. They suggested that the source of energy of the bubbles is sporadic energy releases in the GC by stellar tidal disruption events (TDEs) near the central SMBH (see also other suggestions such as active star formation near the GC, e.g., [17]). The motion of nearby stars orbiting around Sgr A* (the GC) has been observed for more than two decades [18,19,20]. Analysis of the motions gave an estimate of about $4.4\times {10}^6{M}_{\odot }$ for the central SMBH. For an illustration of orbits of stars around Sgr A*, see Figure 16 of Gillessen et al. [20]. For the motion of stars orbiting Sgr A*, the reader is referred to https://www.eso.org/public/videos/eso1825e (accessed on 10 September 2024) (or https://www.youtube.com/watch?v=TF8THY5spmo (accessed on 10 September 2024)); and for animation of stellar orbits around Sgr A* to https://www.eso.org/public/videos/eso1825f (accessed on 10 September 2024) (or https://www.youtube.com/watch?v=wyuj7-XE8RE (accessed on 10 September 2024)), and https://www.youtube.com/watch?v=tMax0KgyZZU (accessed on 10 September 2024).

A TDE occurs when a star is coming too close to an SMBH (closer than the tidal radius). The classic picture is that the star is disrupted by the tidal force, and after half of the stellar debris is in unbound orbits and the other half in bound orbits, it falls back towards the black hole [21]. The problem is how and how much energy is released into the host galaxy. For instance, how much energy is carried away by the unbound debris, and how much binding energy is released by the bound debris, say through accretion. Theoretically, it is possible to have $5\%$ of the rest mass energy of the star being released (for a solar type star this is about ${10}^{53}$ erg). However, the results from observations are mixed, from several ${10}^{51}$ to ${10}^{53}$ erg, It is an active area of research to study TDEs from different perspectives (to name a few, e.g., [22,23,24,25,26,27,28,29,30,31,32]). The reader is also referred to reviews like Dai et al. [33] and Gezari [34]. Another issue is the rate of TDEs at the center of a galaxy with an SMBH. It is conceivable that the rate depends on the type of galaxy and the environment near the black hole. A typical estimation is roughly ${10}^{-4}$∼${10}^{-5}$ ${\mathrm{yr}}^{-1}$ (e.g., [35,36,37]).

Ko et al. [16] adopted the TDE model of Dai et al. [23] and obtained the outflow energy of an event as about ${10}^{52}$∼${10}^{53}$ erg (the outflow velocity is about $0.1$∼$0.3$ the speed of light). Together with an event rate of about ${10}^{-4}$ ${\mathrm{yr}}^{-1}$, this provides an average power of about $\dot{\mathcal{E}}\sim 3\times {10}^{41}$ erg ${\mathrm{s}}^{-1}$, which would be sufficient to power the FBs. Miller & Bregman [38] inferred a bubble expansion rate of 490 km ${\mathrm{s}}^{-1}$, an age of 4.3 Myr, and a luminosity $2.3\times {10}^{42}$ erg ${\mathrm{s}}^{-1}$ (see also the review of FBs in [7]).

Recently, evidence of the energy release at Sgr A* was interpreted as the result of the latest stellar disruption. Two elongated chimneys of about 150 pc near the GC were found in the X-ray [39,40] and radio [41] ranges. In both cases, the total energy of the chimneys was estimated to be about ${10}^{53}$ erg or below, which might be a result of the latest TDE by the central SMBH.

Similar processes of stellar disruption at the galactic centers of other galaxies have been observed. For example, X-ray transient Swift J164449.3+573451 (also known as GRB 110328A) was detected by Swift in the direction of the constellation Draco with a peak luminosity ${10}^{48}$ erg ${\mathrm{s}}^{-1}$ (see Figure 3). Observations showed that the transient originated from the center of a galaxy at cosmological distances involving an SMBH in the galaxy nucleus. It was concluded that Swift J164449.3+573451 most likely originated from the central SMBH, and the X-ray and radio emissions were interpreted as a result of stellar capture by the black hole (see, e.g., [22,32,42]).

2. Structure of the Fermi Bubbles

A sudden (sporadic) energy release by a TDE in the GC creates a cavity with a shock which expands into the surrounding nonuniform medium of the halo. For example, the gas distribution in the halo above (and below) the galactic plane decays exponentially with a scale height $H=2$ kpc and the density at the plane is $n_0=4\times {10}^{-3}{\mathrm{cm}}^{-3}$ (see, e.g., [44]). The exponential gas distribution is

$$\mathcal{R}\left(z\right)=exp\left(-{\displaystyle \frac{z}{H}}\right),$$

(1)

where $\mathcal{R}\left(z\right)=n\left(z\right)/n_0$. Alternatively, Miller & Bregman [38] suggested a so-called $\beta$-model of the gas density profile in the halo from the intensity of absorption lines:

$$\mathcal{R}\left(z\right)={\left({\displaystyle \frac{z}{z_c}}\right)}^{-3\beta },$$

(2)

where $z_c=0.26$ kpc and $n_0=0.5{\mathrm{cm}}^{-3}$.

The formalism of envelope propagation was developed as a solution for a strong explosion (see [45]), and was elaborated by Kompaneets [46] for explosion in a nonuniform atmosphere. The reader is referred to the review of Bisnovatyi-Kogan & Silich [47] and the monograph of Zel’dovich & Raizer [48]. This model was derived for a thermonuclear explosion in the terrestrial atmosphere. Figure 4 shows an example of a thermonuclear explosion test in the terrestrial atmosphere.

Kahn [49], Baumgartner & Breitschwerdt [50], Ko et al. [16], and Schulreich & Breitschwerdt [51] developed analytical solutions of a hydrodynamic model for shock wave propagation in nonuniform atmospheres or halos for different energy input rates for single and successive explosions. The shock envelope generated has a double-bubble structure in the halo (see Figure 5).

Following the Kompaneets formalism (see details in [47]), the shock front is described as

$${\left({\displaystyle \frac{\partial r}{\partial y}}\right)}^2-{\displaystyle \frac{1}{\mathcal{R}\left(z\right)}}\left[{\left({\displaystyle \frac{\partial r}{\partial z}}\right)}^2+1\right]=0,$$

(3)

where $y\left(t\right)$ is a transformed time (in units of length),

$$y={\int }_0^t\sqrt{{\displaystyle \frac{\left({\gamma }_{\mathrm{g}}^2-1\right)}{2}}{\displaystyle \frac{2\mathcal{E}\left(t\right)}{3{\rho }_{0}V\left(t\right)}}}dt,$$

(4)

and the bubble volume $V\left(t\right)$ is

$$V\left(t\right)=\pi {\int }_0^{z_u}{r}^2(z,t)dz,$$

(5)

with $r(z,t)$ or $r(z,y)$ as the bubble radius at the altitude z, ${\rho }_0=n_0{m}_p$ the mass density corresponding to $n_0$ (the number density at the base $z=0$), $\mathcal{E}$ the energy released by the central source into the bubble, and ${\gamma }_{\mathrm{g}}$ the adiabatic index of the gas.

For the exponential gas distribution given by Equation (1), the top of the bubble $z_u$ is a function of time t:

$$z_u=-2Hln\left(1-{\displaystyle \frac{y}{2H}}\right).$$

(6)

In this model of Baumgartner & Breitschwerdt [50] and Schulreich & Breitschwerdt [51], the velocity at the top of the bubble $v_u$ for the total energy release of $\mathcal{E}$ is

$$v_u={\displaystyle \frac{d{z}_u}{dt}}=exp\left({\displaystyle \frac{z_u}{2H}}\right){\displaystyle \frac{dy}{dt}}.$$

(7)

In the early phase, when $z_u\lesssim H$, the expanding cavity can be described by the Sedov solution [45], where the gas density is almost uniform and the velocity of the shock envelope decreases with time. When $z_u>H$, the shock propagates in the exponential halo with subsequent acceleration (see [50]).

The propagation of the shock envelope is derived under the strong shock assumption. In reality, if the velocity of the envelope is below the sound speed $c_s$ of the halo gas, then the shock or the envelope will decay and be absorbed in the halo. On the other hand, if this velocity is higher than the sound speed, the shock will be able to penetrate into the halo and transfer the energy from the initial central source into the exponential halo. The velocity of the top of the bubble is the fastest; Baumgartner & Breitschwerdt [50] defined a condition of shock penetration into the exponential halo: $v_u\left(y_{\mathrm{acc}}\right)>3c_s$, where $v_u\left(y_{\mathrm{acc}}\right)$ is the minimum of $v_u$, which occurs at $y=y_{\mathrm{acc}}$, i.e., ${\dot{v}}_u\left(y_{\mathrm{acc}}\right)={\ddot{z}}_u\left(y_{\mathrm{acc}}\right)=0$. The acceleration at the top of the bubble is (see [50])

$${\ddot{z}}_u={\displaystyle \frac{H}{t_{\mathrm{SN}}^2}}{\displaystyle \frac{\left({\gamma }_{\mathrm{g}}^2-1\right)}{4(1-\tilde{y}/2)\tilde{V}}}\left[{\displaystyle \frac{1}{\left(1-\tilde{y}/2\right)}}-{\displaystyle \frac{1}{\tilde{V}}}{\displaystyle \frac{d\tilde{V}}{d\tilde{y}}}\right],$$

(8)

where

$$\tilde{y}={\displaystyle \frac{y}{H}},\tilde{V}={\displaystyle \frac{V}{H^3}},t_{\mathrm{SN}}=\sqrt{{\displaystyle \frac{3{\rho }_0{H}^5}{2{\mathcal{E}}_{\mathrm{SN}}}}}.$$

(9)

As an example, in the left panel of Figure 6 we present the development of the shock velocity at the top of the bubble from a single energy explosion $\mathcal{E}$. The energy of the explosion occupies more and more volume of the exponential atmosphere, finally approaching infinity in finite time (provided that the velocity is always larger than $3c_s$).

For the parameters in the GC, a single star disruption event provides no more than ${10}^{52}$∼${10}^{53}$ erg (see [23,31,33]). This is not enough energy for the FBs or similar structures. An unusually huge single energy release in the past, say exceeding $\mathcal{E}>{10}^{54}$ erg, may explain the origin of the Fermi Bubbles.

Alternatively, this huge energy could be supplied by a series of many weaker disruption events with an effective power input of $\dot{\mathcal{E}}\ge {10}^{40}$ erg ${\mathrm{s}}^{-1}$; see the right panel of Figure 6. This may be interpreted as routine TDEs, each of which produced an energy of ${10}^{52}$∼${10}^{53}$ erg, with the average rate of stellar capture about ${10}^{-5}$∼${10}^{-4}$ ${\mathrm{yr}}^{-1}$ (see [16]).

The envelope shell is composed of the swept-up gas of the bubble, and it is much thinner in comparison to the size of the bubble. The shell thickness is defined as

$$d\left(y\right)={\displaystyle \frac{M_s\left(y\right)}{2\pi {\rho }_{\mathrm{sh}0}{\int }_0^{z_u\left(y\right)}{e}^{-z/H}r(z,y)\sqrt{1+{\left(\frac{\partial r}{\partial z}\right)}^2}dz}},$$

(10)

where $M_s$ is the total mass of the FB,

$$M_s\left(y\right)=\pi {\rho }_0{\int }_0^{z_u\left(y\right)}{e}^{-z/H}{r}^2(z,y)dz,$$

(11)

and ${\rho }_{\mathrm{sh}}={\rho }_{\mathrm{sh}0}exp(-z/H)$ is the density within the shell (see [51]).

Figure 7 shows some examples of numerical simulations of the FB envelope from sporadic star disruptions or from a single huge explosion.

3. Envelope Disruption by Rayleigh–Taylor Instability

The interface between a denser fluid supported by a lighter fluid in a gravitational field is susceptible to Rayleigh–Taylor (RT) instability. The amplitude of an infinitesimal perturbation will grow exponentially in the early phase or the linear phase (see, e.g., [52]). The growth time of the instability in the linear phase is

$${\tau }_{\mathrm{RT}}=\sqrt{{\displaystyle \frac{\lambda }{2\pi g}}{\displaystyle \frac{({\rho }_2+{\rho }_1)}{({\rho }_2-{\rho }_1)}}},$$

(12)

where $\lambda$ is the wavelength of the perturbation, g is the gravitational acceleration, and ${\rho }_1$ and ${\rho }_2$ are the densities of the lighter and denser fluids, respectively. For an illustration of the time evolution of the RT instability, the reader is referred to Figure 4 of Schulreich & Breitschwerdt [51].

In the case of a superbubble [50,51], the RT instabilities are excited between the dense shell and the hot interior when the envelope is accelerating into the exponential halo (see Equation (8)). Identifying the gravitational acceleration with the acceleration at the top of the bubble ${\ddot{z}}_u\left(y\right)$ gives the growth time of the instability (in these coordinates):

$${\tau }_{\mathrm{RT},z_u}\left(y\right)=\sqrt{{\displaystyle \frac{d\left(y\right)}{2\pi {\ddot{z}}_u\left(y\right)}}{\displaystyle \frac{[{\rho }_{\mathrm{sh}}\left(y\right)+{\rho }_{\mathrm{in}}\left(y\right)]}{[{\rho }_{\mathrm{sh}}\left(y\right)-{\rho }_{\mathrm{in}}\left(y\right)]}}},$$

(13)

where y is the transformed time (see Section 2). When the wavelength $\lambda$ of the RT fluctuations is about the envelope shell thickness d, these instabilities may destroy the bubble; see Figure 13 of Schulreich & Breitschwerdt [51].

The temporal evolution of the RT instability during the nonlinear regime is obtained by numerically solving the ordinary differential equation for the RT fluctuations of $\lambda$,

$$\dot{\lambda }\left(y\right)=2\sqrt{\alpha {\ddot{z}}_u\left(y\right)\lambda \left(y\right)},$$

(14)

the parameter $\alpha$ is estimated from the initial condition $\lambda \left(y_0\right)\sim 0.01d\left(y_0\right)$ (see [50]).

4. Energy and Spectrum of Hydrodynamic Fluctuations

From Equation (14), we can estimate the fraction of the total energy of the FBs, $\dot{\mathcal{E}}\sim 3\times {10}^{41}$ erg ${\mathrm{s}}^{-1}$, that is transformed into the hydrodynamic turbulence in the envelope excited by the RT instabilities there. From Landau & Lifshitz [53], we obtain the rate of energy dissipation in the turbulent flux:

$$\epsilon \left(t\right)={\displaystyle \frac{v_{\lambda }^3}{\lambda }}.$$

(15)

Instead of $\lambda$, we can introduce a wavenumber $k=2\pi /\lambda$. (the Kolmogorov–Obukhov spectrum of turbulence). Then, the kinetic energy spectrum $W\left(k\right)$ of the turbulence is

$$W\left(k\right)\sim {\epsilon }^{2/3}{k}^{-5/3},$$

(16)

and

$${\int }_k^{\infty }W\left(k\right)dk\sim v_{\lambda }^2.$$

(17)

In this case, a universal energy spectrum $W\left(k\right)$ is developed in the inertial range, as in Equation (16), i.e., the Kolmogorov–Obukhov spectrum of turbulence Landau & Lifshitz [53].

The energy losses of RT, $\epsilon \left(t\right)$, in the envelope, transferred into the turbulence, is

$$\epsilon \left(t\right)={\displaystyle \frac{{\dot{\lambda }}_0^3\left(t\right)}{{\lambda }_0\left(t\right)}}=8{\alpha }^{3/2}{\ddot{z}}_u^{3/2}{\lambda }_0^{1/2},$$

(18)

where ${\lambda }_0=2\pi /k_0$ is the pumping scale (see below).

The spectrum of hydrodynamic fluctuations $W\left(k\right)$ has a wide inertial interval, where the energy is supplied on the initial scale of $k_0<k$, in which only energy transfer along the spectrum is realized. In this energy range, the spectrum at $\lambda <{\lambda }_0=2\pi /k_0$ is determined only by the power of energy pumping at the scale ${\lambda }_0$.

The growth of this large-scale structure may be understood at each stage in terms of ring vortex pairing, mutual orbiting, and merging, followed by these processes repeating with the just-merged eddies on a larger scale. During this development, ambient material is entrained and intense smaller-scale turbulence is generated in the regions between the vortices, presumably establishing the turbulent cascade to higher wavenumbers, which is eventually dominated by viscosity on the Kolmogorov microscale, $l_K$. The small-scale of the stretched turbulence-generating regions between the vortices can be associated with the Taylor microscale, $l_T$.

The total fluid turbulence energy is given by $E_t=\rho v_t^2$, where $\rho$ is the fluid density and $v_t$ is the root mean square of the velocity of the turbulence (or simply the velocity of the turbulence) at the largest scale ${\lambda }_0$ or $l_0$. The range is taken to extend from $l_0$ down to the effective damping (or Kolmogorov) scale $l_K$. The Reynolds number in this case is

$$Re={\left({\displaystyle \frac{l_0}{l_K}}\right)}^{4/3},$$

(19)

where the smallest scale, $l_K$, is defined by the dissipation of turbulence.

There is no absolute definition of the scale $l_T$ of this transition from the hydrodynamic spectrum to that of MHD (upper limit of MHD turbulence scale). Following Eilek & Hendriksen [54], the Taylor length, $l_T$, was estimated as

$$l_T\sim l_0{\left({\displaystyle \frac{15}{Re}}\right)}^{1/2}.$$

(20)

Therefore, the fluid cascade extends from the Taylor scale $l_T$, on which the transition from a large-scale ordered turbulence to smaller-scale disordered motion occurs. The cascade proceeds down to the smallest scale, $l_K$, determined by dissipation:

$$l_K\approx l_0{\left(Re\right)}^{-3/4}.$$

(21)

We adopt the hydrodynamic view that the Kolmogorov equilibrium cascade exists between $l_T$ and $l_K$ as our first approximation to the complex interactions expected in the regime. The relation between (the root mean square of) the turbulence velocity $v\left(l\right)$ and the turbulence scale l is

$$v\left(l\right)=v\left(l_T\right){\left({\displaystyle \frac{l}{l_T}}\right)}^{1/3}.$$

(22)

5. Particle Acceleration by Alfvén Fluctuations and the Lighthill Radiation

Contrary to shock acceleration, turbulent resonant acceleration does not require strong shocks. The fluid turbulence in the interstellar medium or intercluster medium is another possible origin of particle acceleration, e.g., in galaxy clusters [55], and in radio jets in turbulent mixing regions [56]. Turbulent motions act as a source of waves and are presented as a hierarchy of eddies. The Lighthill mechanism acts as a direct source of energy to the MHD waves over the range of wavenumbers corresponding to the fluid turbulent spectrum. A turbulent eddy has kinetic energy, which it releases when it mixes with its surroundings. Most of this energy is returned to the ambient medium, but a small fraction is transformed into propagating waves (see, e.g., [56,57,58,59,60] and others). A strong coupling between particle energy and turbulent energy spectra can be expected, and the hydrodynamic turbulence in the medium accelerates particles through wave–particle resonance. Alfvén waves are an alternative source of charged-particle acceleration via resonant interaction of MHD waves with relativistic particles.

In a pioneering paper, Lighthill [61] developed the model of acoustic waves, which are excited by hydrodynamic turbulence in the absence of magnetic fields. This is known as the Lighthill radiation. The radiated power of the waves is roughly the energy density in the turbulent motion, $\epsilon$, divided by the decay time scale of waves $\tau$, and the compactness of the eddy (which is measured by the ratio of the size of the eddy, l, to the wavelength, $\lambda =2\pi /k$).

Kulsrud [62] showed that this radiation of MHD turbulence is excited if there is an external constant magnetic field. If there is no external magnetic field, the magnetic turbulence generates sound waves only via the Lighthill mechanism. If there is a constant external magnetic field, hydromagnetic waves are generated instead by Alfvén waves, unless the energy density of hydrodynamic turbulence prevails over the energy of magnetic density. The central idea of a coupling between the hydrodynamic eddy cascade and the MHD waves through the process of Lighthill radiation is presented in Kulsrud [62], Parker [63], and Kato [64].

Kato [64] developed the Lighthill theory of MHD radiation for strong and weak magnetic fields B, which is characterized by the magnetic Mach number, the ratio of the (root mean square) velocity of the hydrodynamic turbulence v to the Alfvén speed $v_A=B/\sqrt{4\pi \rho }$:

$$M_A={\displaystyle \frac{v}{v_A}}.$$

(23)

For small Mach numbers ($M_A\ll 1$), a small fraction of the power is emitted in the form of Alfvén waves, while for large Mach numbers ($M_A\gg 1$) the power is emitted as sound waves and the radiation of Alfvén waves is insignificant. The turbulence decay time scale is a nonlinear cascade time, which is the eddy turnover time, i.e., the eddy size l divided by its velocity v:

$$\tau \approx {\displaystyle \frac{l}{v}}.$$

(24)

For Alfvén waves in a strong magnetic field, and the frequency $\omega =\overline{k}{v}_A$,

$$\overline{k}l\approx {\displaystyle \frac{v}{v_A}},$$

(25)

which corresponds to the resonance $\tau \left(l\right)\approx 1/\omega \left(\overline{k}\right)$.

From radio polarization measurements, Zhang et al. [17] showed that there are large-scale magnetic fields in the Fermi and eROSITA bubbles. Figure 8 shows several kpc-scale magnetized structures in the bubbles.

If the turbulent magnetic field dominates the motion ($M_A\ll 1$), then the power is

$$P_A\sim {\displaystyle \frac{\rho v^3}{l}}{M}_A,$$

(26)

and a small fraction of the hydrodynamic flux is transformed into Alfvén waves. In the following, we present some details of the radiation of MHD waves in the space medium in the limit of strong external magnetic fields.

We assume that fluid turbulence is induced by the motion of a smaller cluster in a larger cluster and its energy spectrum is described by a power law (see [54,59]):

$$W_f\left(k\right)=W_f^0{k}^{-m},$$

(27)

where $k=2\pi /l$ is the wavenumber corresponding to the scale l, $W_f\left(k\right)\delta k$ is the energy per unit volume in turbulence with wavenumbers between k and $k+\delta k$, and $W_f^0$ and m are constants. If one expresses the turbulent spectrum in terms of eddy size, the spectrum is represented by $W\left(l\right)\propto l^{m-2}$. The cascade of the fluid turbulence extends from a largest eddy size $l_0=2\pi /k_0$ down to a smallest scale determined by dissipation: $l_K\sim l_{0}R{e}^{-3/4}$ (where $Re$ is the Reynolds number). Since most of the energy of fluid turbulence resides in the largest scale, the total energy density of fluid turbulence is presented in the form $E_t\sim \rho v_t^2$, where $\rho$ is the fluid density and $v_t$ is the turbulent velocity of the largest scale $l_0$. The normalization $W_f^0$ can be derived from the relation $E_t={\int }_{k_0}^{k_T}{W}_f\left(k\right)dk$,

$$W_f^0={\displaystyle \frac{E_t}{R}}{k}_T^{(m-1)},$$

(28)

where

$$\begin{array}{ccc}\hfill & R& ={\displaystyle \frac{1}{(m-1)}}\left[{\displaystyle \frac{k_0{W}_f\left(k_0\right)}{k_T{W}_f\left(k_T\right)}}-1\right]\hfill \\ \hfill & & \approx {\displaystyle \frac{1}{(m-1)}}{\displaystyle \frac{k_0{W}_f\left(k_0\right)}{k_T{W}_f\left(k_T\right)}}.\hfill \end{array}$$

(29)

Here $k_T=2\pi /l_T$ and $l_T$ is the wavelength below which Alfvén waves are driven.

A fluid eddy of size l has a velocity

$$\begin{array}{ccc}& v\left(l\right)& \approx {\left[{\displaystyle \frac{l{W}_f\left(l\right)}{\rho }}\right]}^{1/2}={\left[{\displaystyle \frac{k{W}_f\left(k\right)}{\rho }}\right]}^{1/2}\hfill \\ & & ={\left({\displaystyle \frac{E_t}{\rho R}}\right)}^{1/2}{\left({\displaystyle \frac{k}{k_T}}\right)}^{(1-m)/2}.\hfill \end{array}$$

(30)

Turbulence on a scale k will radiate Alfvén waves at the wavenumber

$$\overline{k}=\left[{\displaystyle \frac{v\left(l\right)}{v_A}}\right]k.$$

(31)

Here, we recall that k is the wavenumber of hydrodynamic turbulence and $\overline{k}$ is the wavenumber of the Alfvén waves.

Let $v\left[k\left(\overline{k}\right)\right]$ be the fluid velocity on the fluid scale $k\left(\overline{k}\right)$ that drives Alfvén waves of wavenumber $\overline{k}$. From Equations (30) and (31), we obtain

$$v\left[k\left(\overline{k}\right)\right]=v_A{\left({\displaystyle \frac{E_t}{\rho v_A^{2}R}}\right)}^{1/(3-m)}{\left({\displaystyle \frac{\overline{k}}{k_T}}\right)}^{(1-m)/(3-m)}.$$

(32)

Assume that the energy going into Alfvén waves at wavenumber $\overline{k}$ has an energy flux

$$I_A\left(\overline{k}\right)=I_0{(\overline{k}/k_T)}^{-s_t},$$

(33)

where $I_A\left(\overline{k}\right)\delta \overline{k}$ is the power per unit volume going into Alfvén waves with wavenumbers in the range $\overline{k}\to \overline{k}+\delta \overline{k}$ and $I_0$ and $s_t$ are constants. In this case, the power per unit volume going into the Alfvén mode from fluid turbulence is

$$P_A={\int }_{\overline{k}}^{{\overline{k}}_{\max }}{I}_A\left(\overline{k}\right)d\overline{k}\approx {\displaystyle \frac{I_0{k}_T}{(s_t-1)}}{\left({\displaystyle \frac{\overline{k}}{k_T}}\right)}^{(1-s_t)},$$

(34)

where $\overline{k}\ll {\overline{k}}_{\max }$ and $s_t>1$.

On the other hand, according to the Lighthill theory, $P_A$ is given by Equation (26):

$$P_A={\eta }_A\left[{\displaystyle \frac{v\left(l\right)}{v_A}}\right]\rho v^3\left(l\right)k,$$

(35)

where ${\eta }_A$ is an efficiency factor of order unity [54,56,64]. By comparing Equations (34) and (35), and using Equations (31) and (32), we obtain (see [54,59])

$$s_t={\displaystyle \frac{3(m-1)}{(3-m)}},$$

(36)

$$I_0={\eta }_A(s_t-1)\rho v_A^3{\left({\displaystyle \frac{E_t}{\rho v_A^{2}R}}\right)}^{3/(3-m)},$$

(37)

The power radiated in the form of Alfvén waves $P_A$ is dominated by small $\overline{k}$ and the smallest is ${\overline{k}}_T$. With Equation (30) and ${\overline{k}}_T=[v\left(l_T\right)/v_A]k_T$, the total power in the form of Alfv’en waves is approximately

$$P_A={\eta }_A{\left({\displaystyle \frac{E_t}{\rho v_A^{2}R}}\right)}^2\rho v_A^3{k}_T.$$

(38)

For Kolmogorov turbulence, the spectral index is $m=5/3$ (see Equation (27)), and the spectral index of the Alfvén wave flux is $s_t=3/2$ (Equation (33)).

6. Spectrum of MHD Turbulence in the Fermi Bubble Envelope

The evolution of the spectrum of the Alfvén waves, $W_k\left(t\right)$, is described by the equation of nonlinear diffusion presented in, e.g., Brunetti et al. [55] and Brunetti & Blasi [65]:

$${\displaystyle \frac{\partial W_k(\overline{k},t)}{\partial t}}={\displaystyle \frac{\partial }{\partial \overline{k}}}\left[D_{kk}{\displaystyle \frac{\partial W_k(\overline{k},t)}{\partial \overline{k}}}\right]-\Gamma \left(\overline{k}\right)W_k(\overline{k},t)+I_A(\overline{k},t).$$

(39)

The first term on the right-hand side of Equation (39) describes the nonlinear MHD wave–wave cascade. The diffusion coefficient for the Kolmogorov and the Iroshnikov–Kraichnan spectra is (e.g., [66])

$$\begin{array}{ccc}\hfill D_{kk}& =& v_A\left\{\begin{array}{c}{\displaystyle {\overline{k}}^{7/2}{\left[\frac{W_k(\overline{k},t)}{2W_B}\right]}^{1/2}},\mathrm{Kolmogorov}\hfill \\ {\displaystyle {\overline{k}}^4\left[\frac{W_k(\overline{k},t)}{2W_B}\right]},\mathrm{Iroshnikov}\text{-}\mathrm{Kraichnan}\hfill \end{array}\right.\hfill \end{array}$$

(40)

where $W_B=B_0^2/8\pi$.

The second term on the right-hand side of Equation (39) describes the damping of MHD waves by collisions of relativistic and thermal particles in the interstellar or intercluster medium (see [67]):

$$\begin{array}{ccc}& {\Gamma }_k& \simeq {\displaystyle \frac{4{\pi }^3{e}^2{v}_A^2}{\overline{k}{c}^2}}{\int }_{p_{\min }}^{p_{\max }}{p}^2(1-{\mu }_{\alpha }){\displaystyle \frac{\partial F(p,t)}{\partial p}}dp\hfill \\ & & ={\displaystyle \frac{{\pi }^2{e}^2{v}_A^2}{\overline{k}{c}^2}}{\int }_{p_{\min }}^{p_{\max }}(1-{\mu }_{\alpha })\left[{\displaystyle \frac{\partial N(p,t)}{\partial p}}-{\displaystyle \frac{2N(p,t)}{p}}\right]dp,\hfill \end{array}$$

(41)

where $F(p,t)$ is the particle distribution function and $N(p,t)=4\pi p^{2}F(p,t)$, and

$${\mu }_{\alpha }={\displaystyle \frac{v_A}{c}}\pm {\displaystyle \frac{m\Omega }{p\overline{k}}}.$$

(42)

Here, the upper and lower signs are for negatively and positively charged particles, respectively.

The third term on the right-hand side of Equation (39), $I_A$, describes the injection of Alfvén waves by the fluid turbulence through the Lighthill mechanism.

The time scale of the damping with the thermal pool is considerably shorter than the cascade time scale for $\overline{k}/{\overline{k}}_{\max }\gg 0.1$. Thus, a break or a cutoff in the wave spectrum is expected at large wavenumbers.

In the following, we present an alternative formulation of the steady-state equation for $W_k\left(\overline{k}\right)$. For simplicity, it is described in a compact form (see [68,69] and references therein):

$${\displaystyle \frac{\partial }{\partial \overline{k}}}\left({\displaystyle \frac{\overline{k}{W}_k}{T_{\mathrm{NL}}}}\right)=2{\Gamma }_{\mathrm{CR}}{W}_k+I_A\left(\overline{k}\right).$$

(43)

Here, the source of MHD fluctuations, $I_A\left(\overline{k}\right)$, is given by Equation (33).

The rate of damping of MHD waves by cosmic rays, ${\Gamma }_{\mathrm{CR}}\left(\overline{k}\right)$, is (see, e.g., [70])

$${\Gamma }_{\mathrm{CR}}\left(\overline{k}\right)={\displaystyle \frac{\pi Z^2{e}^2{v}_A^2}{2\overline{k}{c}^2}}{\int }_{p_{\mathrm{res}}\left(\overline{k}\right)}^{\infty }{\displaystyle \frac{dp}{p}}F\left(p\right),$$

(44)

where $F\left(p\right)$ is the CR distribution, p is the particle momentum, and $p_{\mathrm{res}}\left(\overline{k}\right)=ZeB/c\overline{k}$.

In the form of the Iroshnikov–Kraichnan spectrum, the term on the left-hand side of Equation (43) is

$${\displaystyle \frac{\partial }{\partial \overline{k}}}\left({\displaystyle \frac{\overline{k}{W}_k}{T_{\mathrm{NL}}}}\right)={\displaystyle \frac{d}{d\overline{k}}}\left[{\displaystyle \frac{C\left({\overline{k}}^3{W}^2\left(\overline{k}\right)\right)}{\rho v_A}}\right],$$

(45)

where the interaction is

$$T_{\mathrm{NL}}^{-1}\left(\overline{k}\right)=C_{\mathrm{NL}}{\displaystyle \frac{{\overline{k}}^2{W}_k\left(\overline{k}\right)}{m_{\mathrm{i}}{n}_{\mathrm{i}}{v}_{\mathrm{A}}}},$$

(46)

and the constant $C_{\mathrm{NL}}\sim 1$.

If the magnetic field fluctuations are injected by an external source at the scale $L=1/{\overline{k}}_L$, Equation (43) can be simplified to

$${\displaystyle \frac{\partial }{\partial \overline{k}}}\left({\displaystyle \frac{\overline{k}{W}_k}{T_{\mathrm{NL}}}}\right)=2{\Gamma }_{\mathrm{CR}}{W}_k+\Phi \delta (\overline{k}-{\overline{k}}_L).$$

(47)

This result has been applied to the spectrum of MHD turbulence in the FB envelope for $\overline{k}>{\overline{k}}_0$ (see [71]). The solution of Equation (47) is given by

$$W_k\left(\overline{k}\right)={\left({\displaystyle \frac{{\overline{k}}_0}{\overline{k}}}\right)}^{3/2}{W}_k\left({\overline{k}}_0\right)-{\displaystyle \frac{Z^2{e}^2{B}^2{v}_A}{8C{c}^2{\overline{k}}^{3/2}}}{\int }_{{\overline{k}}_0}^{\overline{k}}{\overline{k}}_L^{-5/2}d\overline{k}{\int }_{p_{\mathrm{res}}\left({\overline{k}}_L\right)}^{\infty }{\displaystyle \frac{F\left(p\right)dp}{p}},$$

(48)

where $W_k\left({\overline{k}}_L\right)={\overline{k}}_L^{-3/2}\sqrt{\rho v_A\Phi /C}$, where $\Phi$ describes the source injection at ${\overline{k}}_L$.

The coefficient of momentum diffusion of CRs is (see [70])

$$D_p\left(p\right)=p^2\kappa \left(p\right),$$

(49)

where

$$\kappa \left(p\right)={\displaystyle \frac{3{\overline{k}}_{\mathrm{res}}^2{W}_k\left({\overline{k}}_{\mathrm{res}}\right)}{\rho v}}.$$

(50)

Here, ${\overline{k}}_{\mathrm{res}}=1/r_L=ZeB/pc$, $r_L$ is the particle Larmor radius, and B is the magnetic field strength.

The momentum diffusion coefficient $D_p$ for the bubble parameters is shown in Figure 9 (solid line). For comparison, the dash-dotted line is the diffusion coefficient for the Kraichnan spectrum of turbulence without CR absorption.

As shown in Figure 9, the wave damping by cosmic rays can terminate the cascade for relatively small CR momenta p. Brunetti et al. [55] showed that the time of damping was considerably shorter than the cascade time for large wavenumbers. They also concluded that the damping rate for protons largely dominates that for electrons. These protons can exhibit a resonance with the relativistic electrons, which may be important for their acceleration.

From Equation (48), the distribution function of CR electrons, $F\left(p\right)$, can be estimated from the observed gamma-ray [2,72] and microwave [3] emissions for the expected parameter of the FB envelope of hydrodynamic turbulence.

If the effect of CR damping is insignificant and $2{\Gamma }_{\mathrm{CR}}{W}_k$ can be ignored in Equation (43), then Equation (43) is simply the balance of wave cascading and MHD excitation by the hydrodynamic turbulence:

$${\displaystyle \frac{\partial }{\partial \overline{k}}}\left({\displaystyle \frac{\overline{k}{W}_k}{T_{\mathrm{NL}}}}\right)\approx I_A\left(\overline{k}\right).$$

(51)

With Equation (33), the spectrum of MHD fluctuations $W_k$ is

$$W_k^2\left(\overline{k}\right)\approx {\displaystyle \frac{\rho v_A}{C_{\mathrm{NL}}{\overline{k}}^3}}\left[{\mathcal{P}}_0-{\displaystyle \frac{I_0{k}_0}{(s_t-1)}}{\left({\displaystyle \frac{k_0}{\overline{k}}}\right)}^{(s_t-1)}\right],$$

(52)

where

$${\mathcal{P}}_0={\displaystyle \frac{I_0{k}_0}{(s_t-1)}}+{\displaystyle \frac{C_{\mathrm{NL}}}{\rho v_A}}{k}_0^3{W}_k^2\left(k_0\right).$$

(53)

A cutoff in the spectrum of the waves $W_k$ can be estimated from the balance between damping and the cascade at large wavenumbers.

Unlike the classical GALPROP code with accepted arbitrary parameters of the kinetic diffusion (e.g., [73,74,75]), we derive the coefficients of the kinetic equation for CRs in the FB envelope from Equation (47) for the MHD turbulence:

$${\displaystyle \frac{\partial F\left(p\right)}{\partial t}}+{\displaystyle \frac{F\left(p\right)}{{\tau }_{\mathrm{esc}}}}-Q(p,z)={\displaystyle \frac{1}{p^2}}{\displaystyle \frac{\partial }{\partial p}}{p}^2\left[D_p\left(p\right){\displaystyle \frac{\partial F\left(p\right)}{\partial p}}-\left({\displaystyle \frac{dp}{dt}}\right)F\left(p\right)\right]+{\displaystyle \frac{\partial }{\partial z}}{D}_{zz}\left(p\right){\displaystyle \frac{\partial F\left(p\right)}{\partial z}},$$

(54)

where $(dp/dt)$ (>0) is the rate of continuous energy losses, ${\tau }_{\mathrm{esc}}$ is catastrophic CR losses or the characteristic time of CR escape from the envelope, $Q(p,z)$ is the internal sources of CRs, $D_{zz}$ is the coefficient of spatial diffusion,

$$D_{zz}\left(p\right)={\displaystyle \frac{v{B}^2}{6{\pi }^2{k}^2{W}_k\left(\overline{k}\right)}}={\displaystyle \frac{2\rho v_A^{2}v}{3\pi {\overline{k}}^{2}W\left(\overline{k}\right)}},$$

(55)

and $D_p$ is the coefficient of momentum diffusion and is described by Equations (49) and (50). Here, $\overline{k}=1/r_L=ZeB/pc$ (cf. Equations (49) and (50)).

Our goal is to derive the spectrum of CRs from a combination of kinetic MHD/CR equations and to estimate the proper and correct coefficients of Equation (54). However, there is still a gap between the correct coefficients of the kinetic equations and some rough estimations of the spatial and momentum diffusion from the observed gamma-ray and microwave emissions from the FBs. At present we are unable to derive reliable numerical values for these coefficients, and try to estimate these parameters roughly from the data, ignoring the equation for the origin of MHD turbulence needed for CR scattering and propagation. These parameters of the spatial and momentum diffusion coefficients have been roughly or arbitrarily estimated, e.g., by weak random waves of a hydromagnetic turbulence (see [15]), by a supersonic turbulence (see [76]), or by simple estimations of electron acceleration from shocks of the FBs (see [13]), etc.

In the following, we describe how to roughly estimate the parameters of the spatial and momentum diffusion coefficients from the observed gamma-ray and microwave emissions from the FBs.

7. Leptonic and Hadronic Origins of the Radiation from the Fermi Bubbles

The origin of CRs in the envelope of the giant bubbles is still an open question. The structure of the bubbles is complicated. It is seen in thermal X-rays as an outer envelope with the following parameters:

$$\begin{array}{ccc}& & \mathrm{Power}\mathrm{of}\mathrm{hydrodynamic}\mathrm{turbulence}\text{:}\sim {10}^{39}\mathrm{erg}{\mathrm{s}}^{-1};\hfill \\ & & \mathrm{Thickness}\mathrm{of}\mathrm{the}\mathrm{envelope}\text{:}\sim 100\mathrm{pc};\hfill \\ & & \mathrm{Scale}\mathrm{of}\mathrm{eROSITA}\mathrm{bubbles};\sim 14\mathrm{kpc};\hfill \\ & & \mathrm{Magnetic}\mathrm{field};\sim 8\times {10}^{-6}\mathrm{G};\hfill \\ & & \mathrm{Gas}\mathrm{density}\text{:}\sim 4\times {10}^{-3}{\mathrm{cm}}^{-3};\hfill \\ & & \mathrm{Alfvén}\mathrm{velocity}\text{:}\sim 3\times {10}^7\mathrm{cm}{\mathrm{s}}^{-1}.\hfill \end{array}$$

An inner envelope of size ∼8 kpc is seen in the nonthermal gamma-ray and microwave emissions (see Figure 10). The microwave emission is evidently produced by the synchrotron losses of relativistic electrons, while the origin of the gamma rays is not clear.

The total gamma-ray luminosity of the bubbles between 100 MeV and 500 GeV is $F_{\gamma }\simeq 4.4\times {10}^{37}$ erg ${\mathrm{s}}^{-1}$ (see [72]). The spectrum can be described by a power law, $d{F}_{\gamma }/d{E}_{\gamma }\propto E_{\gamma }^{-1.87}$, with a cutoff $E_{\mathrm{cut}}\simeq 113$ GeV.

7.1. Origin of Gamma-Ray Emission from the Fermi Bubbles

The spectrum of gamma rays can be fitted by either a leptonic or hadronic model.

• Leptonic model: The rate of gamma-ray production by relativistic electrons interacting with low-energy interstellar photons is (see [72])

  $${\epsilon }_{\mathrm{IC}}\left(E_{\gamma }\right)=c\sum _i{n}_i\left({ϵ}_{\mathrm{ph}}\right)\int {\displaystyle \frac{d{\sigma }_{\mathrm{IC}}(E_{\gamma },E_e,{ϵ}_{\mathrm{ph}})}{d{E}_{\gamma }}}{N}_e\left(E_e\right)d{E}_e,$$

  (56)

  where ${\sigma }_{\mathrm{IC}}$ is the inverse Compton (IC) cross-section (see, e.g., [78]). The parameters of the CR electron spectrum were derived from the observed gamma-ray emission from the envelope: $N_e\propto E_e^{-2.17}$, $E_{\mathrm{cut}}\sim 1.25$ TeV. The required total energy in electrons above 1 GeV is ${\mathcal{E}}_e\sim {10}^{52}$ erg.
• Hadronic model: Gamma rays can be produced by proton–proton (p-p) collisions. For calculations of the emission from p-p collision, Ackermann et al. [72] used the p-p cross-section from Kamae et al. [79]. The rate of gamma-ray production by p-p collisions is

  $${\epsilon }_{\mathrm{pp}}\left(E_{\gamma }\right)=c\int {\displaystyle \frac{d{\sigma }_{\mathrm{pp}}(E_{\gamma },E_p)}{d{E}_{\gamma }}}{n}_H{N}_p\left(E_p\right)d{E}_p.$$

  (57)
  The required spectrum of CR protons is expressed as $d{N}_p/d{E}_p\propto E_p^{-2.13}exp(-E_p/E_{\mathrm{cut}})$, where $E_{\mathrm{cut}}\sim 14$ TeV. The total required energy in CR protons above 1 GeV is ${\mathcal{E}}_p\sim 3.5\times {10}^{55}$ erg for $n_H=0.01{\mathrm{cm}}^{-3}$.

In principle, we can interpret the gamma-ray emission (the right panel of Figure 10) by both leptonic (IC) and hadronic (p-p) models for the correspondingly derived parameters of CRs. The question is whether the observed microwave spectrum from the FBs is also compatible with the leptonic or hadronic model.

7.2. Microwave Origin in Cosmic-Ray Electron Model

The origin of the microwave radiation (see left panel of Figure 10) was analyzed for the leptonic and hadronic models (see [72]), with the goal of fitting the gamma-ray and microwave observations within the same model.

The electrons in the IC scenario should also produce the observed WMAP and Planck microwave spectrum and flux. Their properties can be derived from the observed density of gamma rays produced by relativistic electrons that interact with the low-energy interstellar photons in the IC scenario for a magnetic field in the FBs in the range $5$ μG to $20$ μG (see [72]). The best-fit magnetic field is about $8.4$ μG. The synchrotron flux of the FBs from these electrons can be estimated as (see [80,81])

$${\Phi }_{\nu }\simeq 4\pi {\displaystyle \frac{\sqrt{3}{e}^3}{m_e{c}^2}}{\int }_{r_{min}}^{r_0}B\left(r\right)r^{2}dr\int {\displaystyle \frac{\nu }{{\nu }_c}}{N}_e\left(E\right)dE{\int }_{\nu /{\nu }_c}^{\infty }{K}_{5/3}\left(\eta \right)d\eta ,$$

(58)

where $K_{\mu }\left(\eta \right)$ is the McDonald function, and

$${\nu }_c(r,E)={\displaystyle \frac{3eB\left(r\right)}{4\pi m_{e}c}}{\left({\displaystyle \frac{E}{m_e{c}^2}}\right)}^2.$$

(59)

${\Phi }_{\nu }\left(\nu \left(E_e\right)\right)$ can be estimated from the density of relativistic electrons $N_e\left(E_e\right)$ in the FB envelope as

$$N_e\left(E_e\right)\simeq {\displaystyle \frac{\chi }{V_0}}{\displaystyle \frac{{\Phi }_{\nu }\left(\nu \left(E_e\right)\right)}{E_e}},$$

(60)

where the parameter $\chi$ is

$$\chi ={\displaystyle \frac{3m_e^3{c}^5{\nu }_0}{e^4{B}_0^2}}.$$

(61)

Although the main contribution to the IC signal comes from electrons at energies $>100$ GeV, while the main contribution to the Planck frequencies comes from electrons between 10∼30 GeV, they match each other in the IC model.

7.3. Microwave Origin in Cosmic-Ray Proton Model

In the pure hadronic model (see [72,82]), the FB synchrotron emission is produced by secondary electrons from collisions of primary protons.

The approximated equations for the spectrum of secondary electrons produced by p-p and knock-on (KO) collisions ($pe$ and $ee$), can be written as (see, e.g., [83,84])

$$\begin{array}{ccc}& & N_{\mathrm{se}}(E_e,r)={\tau }_e{\int }_{{\displaystyle \frac{m_p}{m_e}}{E}_e}^{\infty }{N}_p(E_p,r)n_H{v}_{p}d{\sigma }_{\mathrm{pp}}\left(E_p\right),\hfill \end{array}$$

(62)

$$\begin{array}{ccc}& & N_{\mathrm{ss}}(E_e,r)={\tau }_e{\int }_{{\displaystyle \frac{m_p}{m_e}}{E}_e}^{\infty }{N}_p(E_p,r)v_p{n}_H{\displaystyle \frac{d{\sigma }_{\mathrm{KO}}(E_p,E^{\prime })}{d{E}^{\prime }}}d{E}_p,\hfill \end{array}$$

(63)

where ${\tau }_e$ is an integral over the rate of electron energy losses:

$${\tau }_e\sim {\int }_{E_e}^{E_{\max }}{\displaystyle \frac{d{E}_e}{{(d{E}_e/dt)}_i}}.$$

(64)

Here ${(d{E}_e/dt)}_i$ can be determined by bremsstrahlung, ionization, synchrotron, or particle escape.

The synchrotron emission of secondary electrons from the FB envelope can be calculated from Equation (58).

The pure hadronic model is unable to reproduce both gamma-ray and radio fluxes from the FBs at the same time. The problem is that the secondary electrons and positrons in the hadronic scenario produce synchrotron radiation with a spectrum that is too soft compared to the microwave haze spectrum, whereas the overall normalization of the synchrotron radiation from the secondary particles is at least a factor of three to four smaller than the microwave level that a hadronic model requires (see [82]).

Thus, we conclude that a purely hadronic origin of the nonthermal emission (gamma and radio) from the FBs is problematic.

8. Number of Relativistic Electrons in the Fermi Bubble Envelope

The origin of relativistic electrons in the FB envelope is an open question. It was assumed that the FB envelopes might be bounded by a shock with the velocity $v_{\mathrm{sh}}\sim {10}^8$ cm ${\mathrm{s}}^{-1}$ (see, e.g., [13,85,86]). It was proposed that these CRs were accelerated at the shock by the standard mechanism of shock acceleration, with the CR spectrum $Q\left(E_e\right)\propto E_e^{-2}$ (see, e.g., [87,88,89,90]). However, eROSITA (see [1]) found that the velocity of the shock with which X-ray giant bubbles propagate is about 340 km ${\mathrm{s}}^{-1}$ and its Mach number is only ≈1.5, which does not correspond to an effective CR shock acceleration. Thus, CRs are not accelerated by a shock near the outer shell of eROSITA. Therefore, CRs should be produced by in situ stochastic acceleration by MHD turbulence $W_k$ near the inner bubble surface (see Section 6 the function $W_k$).

In the following, we focus on the in situ stochastic (Fermi) acceleration of CRs by a hydromagnetic/supersonic turbulence.

8.1. Electrons Accelerated from Background Plasma in the Fermi Bubbles

In the model of Cheng et al. [71], CR electrons can be directly accelerated from a background plasma. They suggested that the acceleration from the background plasma is able to explain the origin of the nonthermal particles responsible for producing the observed fluxes of radio and gamma-ray emissions from the bubbles.

The kinetic equation for the distribution function of electrons, $F(p,t)$, in the case of in situ acceleration is described as

$${\displaystyle \frac{\partial F(p,t)}{\partial t}}+{\displaystyle \frac{F(p,t)}{{\tau }_{\mathrm{esc}}}}={\displaystyle \frac{1}{p^2}}{\displaystyle \frac{\partial }{\partial p}}{p}^2\left[-{\left({\displaystyle \frac{dp}{dt}}\right)}_{\mathrm{C}}F(p,t)+\left\{D_{\mathrm{C}}\left(p\right)+D_{\mathrm{F}}\left(p\right)\right\}{\displaystyle \frac{\partial F(p,t)}{\partial p}}\right].$$

(65)

The distribution function includes the thermal and nonthermal components of the particle distribution. The coefficient ${(dp/dt)}_{\mathrm{C}}$ describes the particle ionization/Coulomb energy losses. $D_{\mathrm{C}}\left(p\right)$ describes the diffusion in the momentum space due to Coulomb collisions (for details, see [91]). The stochastic (Fermi) acceleration is described as diffusion in the momentum space with the diffusion coefficient $D_{\mathrm{F}}\left(p\right)$. ${\tau }_{\mathrm{esc}}$ is the lifetime of particles in the region of acceleration, e.g., due to escape from the region.

In an ionized plasma, the equilibrium (Maxwellian) spectrum of background charged particles is formed by Coulomb collisions. For the case of CR acceleration, there is a boundary, $E=E_{\mathrm{inj}}$, between the equilibrium Maxwellian distribution and a power-law nonthermal spectrum of accelerated particles:

$${\displaystyle \frac{dE}{dt}}={\alpha }_{0}E-{\nu }_{0}E{\left({\displaystyle \frac{kT}{E}}\right)}^{3/2}.$$

(66)

The energy of injection is

$$E_{\mathrm{inj}}\sim kT{\left({\displaystyle \frac{{\nu }_0}{{\alpha }_0}}\right)}^{2/3}.$$

(67)

Here, the parameters ${\alpha }_0$ and ${\nu }_0$ are the acceleration and the ionization loss (by Coulomb collisions), respectively.

Making use of Equation (65), Gurevich [92] studied the process of the formal connection between the equilibrium Maxwellian distribution of background particles and a power-law non-equilibrium spectrum of accelerated particles.

For slow time variations, the runaway flux of particles from the region of thermal Maxwellian distribution into the acceleration range is (see, e.g., [93])

$$S_{\zeta }=S_0\left(t\right)\sqrt{{\displaystyle \frac{2}{\pi }}}{\int }_0^{\zeta }{\zeta }^{2}exp\left(-{\displaystyle \frac{{\zeta }^2}{2}}\right)d\zeta ,$$

(68)

where $\zeta =p/p_0$ is the normalized momentum, and

$$S_0\left(t\right)=\sqrt{{\displaystyle \frac{2}{\pi }}}n\left(t\right)exp\left(-{\int }_0^{\infty }{\displaystyle \frac{\zeta d\zeta }{(1+{\alpha }_0{\zeta }^5/{\nu }_0)}}\right),$$

(69)

and $n\left(t\right)$ is the slow variation in the gas density.

In a non-equilibrium case, this process forms an escape flux of runaway particles for energy $E>E_{\mathrm{inj}}$. On the other hand, this process also forms an excess density at energies $E\lesssim E_{\mathrm{inj}}$, which distorts the thermal Maxwellian distribution. This can be interpreted as a “second” effective temperature, higher than the gas equilibrium temperature.

This model of Gurevich [92] has been applied to the processes of particle acceleration from background plasma in galaxy clusters (see [93,94,95]), in the GC (see [71]), and in the galactic disk [96], in which excesses above thermal particles in the X-ray range is expected.

For example, the spectrum of X-ray emission from the galactic plane can be described as a multi-temperature emission. Dogiel et al. [96] interpreted this X-ray emission as the flux of runaway particles from background gas as a common effect of Coulomb collision (ionization losses) and stochastic acceleration (see Figure 11). When compared with the spectrum of a simple combination of thermal and nonthermal gas, the spectrum with runaway flux is larger, in particular, in the transition range between the thermal and nonthermal parts.

The model of particle excess from background gas proposed for CRs by Dogiel [93] was challenged by Petrosian [97], Wolfe & Melia [98], and Petrosian & East [99]. The problem was that the stochastic acceleration of the accelerated particles from a background plasma would (over) heat the plasma, because their energy would be quickly dumped into the thermal plasma by ionization losses. The energy gained by the particles is distributed to the whole plasma on a time scale much shorter than that of the acceleration process itself. As a result of the relative inefficiency of bremsstrahlung for cooling the accelerated electrons, this tail is quickly dumped into the thermal body of the background plasma (plasma overheating without a prominent tail of accelerated particles). This effect completely prevents the formation of nonthermal spectra from background plasma.

However, Chernyshov et al. [100] showed that the effect of overheating depends on the parameters of acceleration. It is insignificant if the stochastic acceleration is effective. This model depends on a value of $p_{\mathrm{inj}}$ and a free parameter of stochastic acceleration $p_0$ in the form of

$$D_{\mathrm{F}}\left(p\right)=D_0{p}^{\varsigma }\theta (p-p_0),$$

(70)

where $D_0$ and $\varsigma$ are constants. In general, $p_{\mathrm{inj}}$ is determined by $D_{\mathrm{F}}\left(p_{\mathrm{inj}}\right)={\left[p{(dp/dt)}_{\mathrm{C}}\right]}_{p_{\mathrm{inj}}}$. In this model, the injection momentum is given by

$$D_0{p}_{\mathrm{inj}}^{(\varsigma -1)}={\left[{\left({\displaystyle \frac{dp}{dt}}\right)}_{\mathrm{C}}\right]}_{p_{\mathrm{inj}}}.$$

(71)

For a high value of the acceleration momentum $p_0$, the runaway flux of thermal particles cools the plasma down from the very beginning. In spite of energy supply by external sources, the plasma temperature drops (analogous to Maxwell demon). The acceleration generates a prominent tail of accelerated particles but the “excess”, that was expected in Gurevich [92], is not produced in the range around $p_{\mathrm{inj}}$; see Figure 12. For $p_0>p_{\mathrm{inj}}$, plasma overheating is insignificant and stochastic acceleration works well. For $p_0<p_{\mathrm{inj}}$, plasma overheating is significant and stochastic acceleration is inhibited.

In any case, numerical calculations showed that the permitted parameters of the FB model are strongly restricted, and the model is unable to explain the observed fluxes of radio and gamma-ray emissions from the bubbles (see [71]).

8.2. Cosmic-Ray Electrons Re-Accelerated in the Fermi Bubbles

CR electrons can be generated by sources in the galactic disk (e.g., supernova remnant shocks, SNR shocks). Due to synchrotron, inverse Compton, and adiabatic losses in the halo, Cheng et al. [101] deemed that these CR electrons with energies above several GeV are unable to reach the height of the FB envelope (which is about 8∼10 kpc). With appropriate parameters in the FB envelope, these electrons can be re-accelerated in situ up to an energy of about ${10}^{12}$ eV, which is needed to reproduce the observed radio and gamma-ray emissions from the FBs and supply the required power.

The steady-state kinetic equation for the relativistic electrons in FBs can be described in the form (see [101])

$$Q\left(p\right)\delta \left(z\right)=-\nabla ·\left[D_s\nabla F-\mathbf{v}F\right]+{\displaystyle \frac{1}{p^2}}{\displaystyle \frac{\partial }{\partial p}}{p}^2\left[\left({\displaystyle \frac{dp}{dt}}-{\displaystyle \frac{\nabla ·\mathbf{v}}{3}}p\right)F-D_p{\displaystyle \frac{\partial F}{\partial p}}\right],$$

(72)

where $F(p,r,z)$ is the particle distribution function, $Q(p,r)\propto p^{-\gamma }$ is the source function of electrons in the galactic disk ($z=0$), $dp/dt=\mu E^2$ is the rate of synchrotron and inverse Compton energy losses, $D_s(p,r,z)$ is the coefficient of spatial diffusion, $D_p(p,r,z)=\kappa p^2$ is the coefficient of momentum diffusion (coefficient of the Fermi re-acceleration), and $\mathbf{v}(r,z)=v_z{\widehat{e}}_z=3v^{\prime }z{\widehat{e}}_z$ ($v^{\prime }=d{v}_z/dz$) is the wind velocity of advection in the halo in the z direction. The effect of the wind advection leads to adiabatic losses of CRs, $dp/dt=-p\nabla ·\mathbf{v}/3$. Consequently, the spectrum of electrons in the FBs in the acceleration region is harder than the case without advection.

Cheng et al. [101] showed that the gamma-ray and radio emissions of the re-accelerated electrons nicely reproduced the Fermi-LAT and Planck data points for the parameters: the spatial diffusion coefficient $D_s={10}^{29}{\mathrm{cm}}^2{\mathrm{s}}^{-1}$, the energy loss rate $\nu =2\times {10}^{-16}{\mathrm{s}}^{-1}{\mathrm{GeV}}^{-1}$, the velocity gradient of the advection in the halo $v^{\prime }={10}^{-15}{\mathrm{s}}^{-1}$, the magnetic field strength is $B=3$ μG, the thickness of the re-acceleration region (say, the FB envelope) is about $\Delta r_{\mathrm{FB}}=60$ pc, and the parameter of re-acceleration in the FBs is $\kappa =2\times {10}^{-14}{\mathrm{s}}^{-1}$.

The total power $\dot{\mathcal{E}}$ needed to be supplied by sources of Fermi re-acceleration in the FBs to produce high-energy electrons is

$$\dot{\mathcal{E}}=-{\int }_0^{\infty }4\pi E{\displaystyle \frac{\partial }{\partial p}}\left(p^2{D}_p{\displaystyle \frac{\partial F}{\partial p}}\right)dp,$$

(73)

where p and E are the particle momentum and particle kinetic energy, respectively.

To define the spectrum of accelerated electrons we estimated the number of GeV electrons that can reach an altitude of several kpc when the effect of advection $v_z\left(z\right)$ is essential in the galactic halo (see [102,103,104,105]). The spectrum of re-accelerated SNR electrons in the FBs is shown in Figure 13 (cf. [101]). In the figure, the thick solid line is the spectrum of CR electrons from their sources in the galactic disk (see, e.g., [70]). When re-acceleration (stochastic acceleration) and adiabatic losses are taken into account, the spectrum of electrons (thin dashed line) becomes harder than that of the spectrum emitted by sources (thick solid line), but softer than the spectrum of pure re-acceleration (thin dash-dotted line). The spectrum for a velocity gradient $v^{\prime }={10}^{-15}{\mathrm{s}}^{-1}$ (see, e.g., [103,105]) is consistent with the observed one (dashed line in Figure 13).

The model of re-acceleration within the envelope coincides nicely with the observed microwave and gamma-ray emissions, shown in the left and right panels of Figure 10, respectively.

In the phenomenological model, the power, $\dot{\mathcal{E}}$, is estimated numerically from the observed FB gamma-ray and microwave fluxes, and is about $\dot{\mathcal{E}}\sim 2\times {10}^{38}$ erg ${\mathrm{s}}^{-1}$. The density of high-energy electrons needed for the observed gamma-ray flux is shown by the shaded gray region in Figure 13.

With appropriate parameters, these electrons can be re-accelerated up to an energy of ${10}^{12}$ eV, which explains the origin of the observed gamma-ray and radio emissions from the FBs in this model. However, although the model gamma-ray spectrum is consistent with the Fermi results, the model radio spectrum in the pure diffusion model is steeper than that observed by WMAP and Planck.

If adiabatic losses due to plasma outflows from the galactic central regions are taken into account, we expect that the spectrum of electrons in the acceleration region will be harder than the one without advection. Our calculations with divergent outflows show that the gamma-ray and radio emissions of the re-accelerated electrons nicely reproduce the Fermi-LAT and Planck data (see Figure 10).

In essence, both gamma-ray and microwave observations can be explained by only one source of high-energy electrons. The basic idea is summarized as follows. CR electrons from SNRs in the galactic disk are re-accelerated in the FB via supersonic turbulence (or multiple shocks). The resulting spectrum is hard, and in the high-energy region the density of electrons exceeds that required to produce the observed gamma-ray emission. With adiabatic loss by the divergent flow, the density reduces but the spectrum is still hard enough to produce the observed microwave emission. Inevitably, some delicate balance or fine-tuning of parameters is needed.

9. Cosmic-Ray Protons Escaping from the Fermi Bubbles into the Galaxy

The fundamental question of the sources of CRs in the galaxy is still open. We present a number of models which may interpret the origin of CRs escaping from the FBs into the galactic disk. The effect can be observed from the spectrum of CRs near Earth.

As described in Berezinsky et al. [70], the classical model of CR origin in the galaxy is that CRs are generated by SNRs in the galactic disk with energies below ${10}^{15}$ eV. They escape into the galactic halo with an effective spatial diffusion coefficient about $D_{\mathrm{G}}=D_{*}{(E/4\mathrm{GeV})}^{0.6}$ ($D_{*}=6.2\times {10}^{28}{\mathrm{cm}}^2{\mathrm{s}}^{-1}$), which is estimated from the observed chemical composition of CRs (for a modern nonlinear model of CR, see [106] and references therein). Plane shock acceleration produces a CR spectrum $F\left(E\right)\propto E^{-2}$ (see, e.g., [87,88,89,90]). The actual CR spectrum observed outside the acceleration region is a result of the process of acceleration in the region together with the process of particle leakage or escape from the region. Furthermore, the maximum energy that can be attained by CR particles depends on the size and/or lifetime of the acceleration region. Suppose $v_{\mathrm{sh}}$ is the shock velocity and $D_{\mathrm{sh}}$ is the spatial diffusion coefficient in the shock vicinity. The minimum value of $D_{\mathrm{sh}}$ follows Bohm’s limit, $D_{\mathrm{sh}}\approx u{r}_L/3$ (where $u\approx c$ is the speed of CR particles). The maximum energy is constrained by the size of the acceleration region becoming smaller than the diffusion length scale of the particles $l_D\sim D_{\mathrm{sh}}/v_{\mathrm{sh}}$, and/or the lifetime of the region being smaller than the acceleration time scale of the particles ${\tau }_{\mathrm{acc}}\sim D_{\mathrm{sh}}/v_{\mathrm{sh}}^2$.

The CR spectrum (above 1 GeV) observed from Earth can be described by a broken power law; see the classic spectrum in Swordy [107] (or Figure 14, which shows the part above ${10}^{13}$ eV). The power-law index is $-2.7$ between ${10}^9$ and ${10}^{15}$ eV, and $-3.1$ between ${10}^{15}$ and several ${10}^{18}$ eV. An energy around ${10}^{15}$ eV is called the ‘knee’, where the spectrum changes from a harder one to a softer one. The spectrum beyond several ${10}^{18}$ eV becomes harder again, and this region is called the ‘ankle’. The apparent cutoff at somewhat less than ${10}^{20}$ eV is commonly attributed to the Greisen–Zatsepin–Kuzmin limit due to the interaction of ultrahigh-energy CRs with the cosmic microwave background [108,109]. At energies smaller than 1 GeV, the spectrum is heavily affected by solar modulation and activity of the Sun. Although the spectrum is a broken power law, the ‘joints’ (say the ‘knee’ and the ‘ankle’) are smooth. The origins of different parts of the spectrum should be somehow related (e.g., [110]).

SNR shocks are believed to be the source of CRs (e.g., [70]). They can produce the spectral index $-2.7$ reasonably well. With the magnetic field around the SNR comparable to the general interstellar field, Lagage & Cesarsky [111] and Berezhko & Völk [112] estimated that the maximum energy of protons from SNRs is about ${10}^{13}$∼${10}^{14}$ eV. However, with instabilities caused by cosmic ray streaming (e.g., non-resonant hybrid instability), Bell [113] and Bykov et al. [114] found that the fluctuated magnetic field can be orders of magnitude larger at the SNR shock (and so $D_{\mathrm{sh}}\ll D_{\mathrm{G}}$). The confinement time of CRs is longer and they can be accelerated up to ${10}^{15}$ eV (see also [115]).

Figure 14. CR spectrum at the Earth as a combination of the contributions from the SNRs in the galactic disk and the stochastic acceleration in the FBs. Figure reproduced from Cheng et al. [116] with permission.

Figure 14. CR spectrum at the Earth as a combination of the contributions from the SNRs in the galactic disk and the stochastic acceleration in the FBs. Figure reproduced from Cheng et al. [116] with permission.

The SNR framework is able to describe the spectrum in the range $E<{10}^{15}$ eV (see Figure 14). However, for $E>{10}^{15}$ eV, we need some other ideas (see, e.g., [76]). Larger and longer lifetime shocks are required, such as superbubbles (like FBs) and strong galactic winds (e.g., [85,86,116,117,118,119]). In the following, we analyze the spectrum of CRs accelerated in the FBs and discuss whether the bubble’s contribution may explain the ‘knee’ steepening. The idea is based on acceleration by multiple shocks. In a multiple-shock system, two length scales are important: (1) the average separation between two shocks $l_{\mathrm{sh}}\sim v_{\mathrm{sh}}{\tau }_{\mathrm{c}}$ (where ${\tau }_{\mathrm{c}}$ is the average time between the creation of two consecutive shocks), and (2) the diffusion length scale at the shock $l_D$ (∼$D_{\mathrm{sh}}/v_{\mathrm{sh}}$). We focus on the regime $l_{\mathrm{sh}}\ll l_D$ (the supersonic turbulence regime), for energies beyond the ‘knee’ ($E>{10}^{15}$ eV).

An illustration of a possible multiple-shock structure in the FBs is shown in Figure 15.

Escape of Cosmic-Ray Protons Re-Accelerated by Supersonic Turbulence Inside the Fermi Bubbles

Cheng et al. [116] suggested an alternative model of CRs in the FBs (see also [11,12,120,121,122]). They assumed that up to several hundred TDEs might have occurred in the past 10 Myr. This would have generated a series of shocks propagating through the central part of the galactic halo, which would produce relativistic CRs via multiple-shock acceleration. The average separation between two shocks is then

$$\begin{array}{ccc}& l_{\mathrm{sh}}& =v_{\mathrm{sh}}{\tau }_{\mathrm{cap}}\hfill \\ & & \simeq 30\left({\displaystyle \frac{{\tau }_{\mathrm{cap}}}{3\times {10}^4\mathrm{yr}}}\right)\left({\displaystyle \frac{v_{\mathrm{sh}}}{{10}^8{\mathrm{cm}\mathrm{s}}^{-1}}}\right)\mathrm{pc},\hfill \end{array}$$

(74)

where ${\tau }_{\mathrm{cap}}$ is the average time between two stellar captures by the SMBH.

We applied the model of Bykov & Fleishman [123] and Bykov & Toptygin [76] for CR acceleration by multiple shocks in the FBs. Under the conditions of supersonic turbulence (multiple-shock structure) the regime of acceleration is characterized by $l_{\mathrm{sh}}\ll l_{\mathrm{D}}$.

The steady state kinetic equation in axisymmetric geometry is

$${\displaystyle \frac{\partial }{\partial z}}\left[D_s(p,r){\displaystyle \frac{\partial F}{\partial z}}\right]+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial }{\partial r}}\left[D_s(p,r)r{\displaystyle \frac{\partial F}{\partial r}}\right]+{\displaystyle \frac{1}{p^2}}{\displaystyle \frac{\partial }{\partial p}}\left[D_p(p,r)p^2{\displaystyle \frac{\partial F}{\partial p}}\right]=-Q(p,r,z).$$

(75)

The spatial diffusion coefficient inside and outside the bubble is

$$D_s(p,r)=D_B\theta (r_{\mathrm{B}}-r)+D_{\mathrm{G}}\theta (r-r_{\mathrm{B}}),$$

(76)

where $r_{\mathrm{B}}$ is the radius of the bubbles. Inside the bubbles, as a result of interaction with supersonic turbulence, $D_s=D_{\mathrm{B}}\approx u{l}_{\mathrm{sh}}/3\approx c{l}_{\mathrm{sh}}/3$. Outside the bubbles, $D_s=D_{\mathrm{G}}$, which is the average diffusion coefficient in the galaxy, i.e., the one described in Berezinsky et al. [70]. The momentum diffusion coefficient is nonzero inside the bubbles only:

$$D_p(p,r)={\kappa }_{\mathrm{B}}{p}^2\theta (r_{\mathrm{B}}-r),$$

(77)

and ${\kappa }_{\mathrm{B}}\sim v_{\mathrm{sh}}^2/D_{\mathrm{B}}$. The parameters in the FBs can be estimated numerically from the observed CR spectrum.

The method of separation variables as in Bulanov et al. [124] and Bulanov & Dogel [125] can be used to solve Equation (75) for the spectrum of CRs generated by SNRs with the standard model of CR propagation and escape in the galactic halo (see, e.g., [70,80,126]). This model describes the observed CR spectrum near Earth in the range below the ‘knee’, ${10}^{15}$ eV (see Figure 14).

Cheng et al. [116] interpreted the CR spectrum near Earth in the energy range from ${10}^{15}$ eV to a few ${10}^{18}$ eV (from the ‘knee’ to the ‘ankle’) as a combined result of acceleration in the FB and escape. The acceleration is provided by the supersonic turbulence in the FB (see [76]). Under the set of parameters in Cheng et al. [116], the ratio of the escape time to the acceleration time is about 1.9.

In order to derive the spectrum near Earth, Cheng et al. [116] matched the solutions for the spectra inside the FB and outside (in the halo) at $r=r_{\mathrm{B}}$ (see Equations (75) and (76)). Basically, the model can reproduce the CR spectrum from ${10}^{13}$ eV to several ${10}^{18}$ eV. However, we should point out that the numerical result described here [116] has free parameters and some physics have been ignored (e.g., adiabatic loss), and the result might not be very solid. It is imperative to perform further investigations.

A brief summary of the idea is as follows. CRs from SNRs (<${10}^{15}$ eV) are re-accelerated in the FB by multiple shocks or supersonic turbulence. As the shocks are larger and live longer in the FB, CRs can be accelerated to much higher energies (up to ${10}^{19}$). Acceleration by multiple shocks inside the bubbles gives a harder spectrum, hence the contribution of FB to CRs within the ‘knee’ (<${10}^{15}$ eV is subordinate to SNR). The high-energy CRs escaping from the bubbles constitute the sole source of CRs in the range between the ‘knee’ and the ‘ankle’ (${10}^{15}$ to several ${10}^{18}$ eV) observed at Earth. In the model, some fine-tuning of parameters is inevitable.

10. Summary

Here, we present a brief summary of our perspective of the Fermi Bubbles at the Galactic Center.

• The key point of the bubbles is a huge energy release of ${10}^{55}$∼${10}^{56}$ erg in the Galactic Center, whose origin is still unknown.
• We assume that the energy source of the bubbles could be a routine tidal disruption of stars near the central supermassive black hole. Each disruption of a star releases a total energy about ${10}^{52}$∼${10}^{53}$ erg. For a typical rate of stellar capture (once every ${10}^4$ years), this can provide a luminosity ≳${10}^{41}$ erg ${\mathrm{s}}^{-1}$ from the Galactic Center. These processes of stellar tidal disruption events can be directly observed in some external galaxies.
• Hydrodynamic models can describe the envelope of bubble propagation in the galactic halo where the gas distribution is nonuniform. The distribution is commonly characterized by an exponential or a power-law function. The observed shape of the Fermi bubbles seems to suggest an exponential halo. If the velocity of the top of the envelope prevails over the sound velocity in the halo, then the envelope may reach the size of 10 kpc.
• The surface of the top of the envelope propagates with acceleration in the halo. As a result, Rayleigh–Taylor instabilities are developed and they will destroy the bubble envelope at the top. We expect excitation of hydrodynamic instabilities and generation of hydrodynamic turbulence there.
• Turbulent motions act as a source of waves, which are manifested as a hierarchy of eddies, and act as a direct source of energy to the MHD waves (via the Lighthill mechanism). For small Mach numbers, a small fraction of the power radiated by the turbulent motion is Afvén waves.
• The coefficients of the spatial and momentum diffusion of the system of nonlinear kinetic equations of the cosmic-ray distribution function are derived from the spectrum of MHD waves. These coefficients were calculated analytically, but we were unable to estimate the numerical values for the bubbles because of the lack of available observations on the wave spectrum.
• We roughly estimated the spatial and momentum diffusion of cosmic rays in the envelope from the data of gamma-ray and microwave radiations from the Fermi bubbles.
• We concluded that the observed gamma-ray and microwave radiations from the envelope of the Fermi Bubbles are generated by cosmic-ray electrons only. The contribution of cosmic-ray protons can be neglected.
• We prefer the model that GeV cosmic-ray electrons from supernova remnants in the galactic disk are re-accelerated in situ in the bubbles to TeV energy range. With the help of a divergent flow, this can reproduce the data of both gamma-ray and microwave observations.
• On the other hand, high-energy cosmic ray protons can escape the bubbles and reach the Earth. Cosmic ray protons from supernova remnants can be accelerated in the bubbles by supersonic turbulence to higher energies. We found that the escaped high-energy protons that arrive the Earth can reproduce the spectrum and flux of cosmic rays in the range ${10}^{15}$∼${10}^{18}$ eV (from the ‘knee’ to the ‘ankle’), observed near Earth.