Reconstruction of Fermi and eROSITA Bubbles from Magnetized Jet Eruption with Simulations

Abstract

The Fermi bubbles and the eROSITA bubbles around the Milky Way Galaxy are speculated to be the aftermaths of past jet eruptions from a supermassive black hole in the galactic center. In this work, a 2.5D axisymmetric relativistic magnetohydrodynamic (RMHD) model is applied to simulate a jet eruption from our galactic center and to reconstruct the observed Fermi bubbles and eROSITA bubbles. High-energy non-thermal electrons are excited around forward shock and discontinuity transition regions in the simulated plasma distributions. The $\gamma$-ray and X-ray emissions from these electrons manifest patterns on the skymap that match the observed Fermi bubbles and eROSITA bubbles, respectively, in shape, size and radiation intensity. The influence of the background magnetic field, initial mass distribution in the Galaxy, and the jet parameters on the plasma distributions and hence these bubbles is analyzed. Subtle effects on the evolution of plasma distributions attributed to the adoption of a galactic disk model versus a spiral-arm model are also studied.

1. Introduction

The Fermi and the eROSITA bubbles recently discovered were speculated to be the aftermaths of past jet eruptions from the center of our Milky Way Galaxy [1,2,3] millions of years (Myrs) ago. A supermassive black hole (SMBH) is commonly observed in the center of a galaxy, which spews powerful kpc-scale jets intermittently [4,5].

Two $\gamma$-ray Fermi bubbles were manifested on the skymap in the range of $\pm {50}^{\circ }$ in latitude, extending about ${40}^{\circ }$ in longitude [2]. The plasma dynamics behind the Fermi bubbles has been studied by applying various hydrodynamic (HD) or magnetohydrodynamic (MHD) models and simulations. However, the time scales and jet parameters in different models may deviate significantly since they cannot be directly observed [6,7,8]. In [6], Fermi bubbles were reconstructed by simulating an active galactic nucleus (AGN) jet that erupted from the SMBH in the galactic center 1.2 Myrs ago with a 3D MHD model. In [7], the morphology of Fermi bubbles was revealed by simulating an AGN jet that erupted about 5 Myrs ago with an HD model, neglecting the effects of the magnetic field. In [8], Fermi bubbles were studied by applying an MHD model, with a jet that erupted about 2.6 Myrs ago.

The jet morphology and the spectrum of Fermi bubbles could be explained in terms of cosmic rays [8], generated via leptonic, hadronic and in situ acceleration processes [9]. A leptonic process may last for a few Myrs and a hadronic process may last for more than 10 Myrs. A hadronic process alone could not account for the observed emission spectrum; hence, a hybrid model including cosmic-ray electrons from other sources was proposed [8]. Leptons, hadrons or their combination are accelerated by shocks to emit $\gamma$-rays, manifesting sharp edges of the bubbles. However, this model predicted constant volume emissivity, which could not explain the flat surface brightness projected on the skymap [9]. An RMHD model may have a better chance to simulate the plasma distributions and subsequently reconstruct the bubbles that match the observation in shape, size and radiation intensity.

In [10,11], the emission from Fermi bubbles was simulated by applying an MHD model coupled with a diffusion–advection equation of cosmic rays. The emission spectra of hadronic, leptonic and hybrid cosmic rays manifest similar features, with stronger emission in X-ray and $\gamma$-ray from hadronic cosmic rays than from leptonic cosmic rays [10]. The X-ray and $\gamma$-ray were mainly contributed via inverse Compton scattering, weakly contributed by bremsstrahlung and negligibly by synchrotron radiation.

Tidal disruption events (TDEs) and winds from hot accretion flow have also been proposed to explain the formation of Fermi bubbles [12,13]. In [12], the energy released by a TDE of stars, triggered by an SMBH, was studied. Routine TDEs near a galactic center could accumulate energy to form large-scale structures like Fermi bubbles. However, the mass of a black hole at our galactic center is not sufficient to drive TDEs out of the surrounding stars. In the simulation of [13], Fermi bubbles were inflated by the winds from past hot accretion flow in Sgr A*, which lasted for about 10 Myrs and quenched 0.2 Myrs ago.

The eROSITA bubbles emitting soft X-ray were observed to manifest a roughly spherical shape, extending in the range of $\pm {80}^{\circ }$ in latitude and about ${80}^{\circ }$ in longitude [3]. The edges of eROSITA bubbles appear brighter and hotter than the ambient, consistent with the features of a forward shock released from a strong energy release event at the galactic center.

The front edge of the eROSITA bubbles and the boundary between eROSITA and Fermi bubbles are clues of past jet eruptions [3]. Four possible mechanisms were proposed to explain the boundary between eROSITA bubbles and Fermi bubbles, including a forward shock from a sequence of energy release events, a reverse shock, a wind-termination shock and a contact discontinuity. If both the eROSITA bubbles and the Fermi bubbles were driven by the same energy release event at the galactic center, the boundary between them would be a contact discontinuity, and the front-edge of the eROSITA bubbles would be a shock plowing through the halo gas.

Bubble-like structures were also observed around external galaxies [14,15]. In [14], a diffuse $\gamma$-ray halo, with a luminosity of $(3.2\pm 0.6)\times {10}^{38}$ (erg/s) in 0.3–100 GeV, was discovered around the Andromeda Galaxy M31 by analyzing the Fermi LAT data over almost 7 years. Two 6–7.5 kpc bubbles, symmetrical about the M31 galactic disc, were revealed by using a best-fit halo template. The bubbles looked similar to the Fermi bubbles around our galactic center. Similar eROSITA bubbles around M31 were also conjectured. In [15], a nuclear superbubble emitting hard X-rays was observed on the southwest side of the galactic center of NGC 3079. No hard X-ray emission was detected on the northeast side, which could not be explained with inverse Compton scattering or stellar X-ray sources. However, the combination of synchrotron emission and thermal bremsstrahlung at about 1 keV could account for the broadband radio/hard X-ray spectra.

Different sets of parameters have been proposed to account for similar observation features of Fermi bubbles and eROSITA bubbles, possibly attributed to the complicated and nonlinear physical processes underpinning these phenomena. In [8], both Fermi bubbles and eROSITA bubbles were simulated under a single jet eruption from the galactic center 2.6 Myrs ago. The morphology and multi-wavelength spectra of the observed bubbles were reconstructed, and different sets of parameters might lead to similar morphology. The reported age of bubbles ranges from 1.2 Myrs to 9 Myrs [6,7,8]. The Fermi bubbles and the eROSITA bubbles would take shape more quickly if the initial mass distribution in the Galaxy was not included. A more general study on the interactions among the background magnetic field, jet parameters and the observation data will improve understanding about the physical mechanisms behind the formation of Fermi bubbles and eROSITA bubbles.

A typical jet speed is about $0.1c$ and a typical plasma flow speed is about 0.01–0.1c [6,7,8]. The plasma flow near the dilute spine region of a jet is sensitive to the gas pressure and can be accelerated closer to the speed of light. Thus, a relativistic magnetohydrodynamic (RMHD) model will be more robust for simulating the Fermi bubbles and the eROSITA bubbles.

In this work, a 2.5D axisymmetric RMHD model is applied to simulate a jet eruption from the center of our Milky Way Galaxy, arousing a forward shock and a tangential discontinuity. The RMHD model is validated by comparing with the results in the literature. The effects of background magnetic field, initial mass distribution in the Galaxy and jet parameters on the evolution of plasma distributions are simulated. Subtle effects on the evolution of plasma distributions attributed to the adoption of a galactic disk model versus a spiral-arm model are also studied. The emissivities of leptonic and hadronic processes from the plasma distribution are investigated to successfully reconstruct the morphology of Fermi bubbles in $\gamma$-ray and eROSITA bubbles in X-ray on the skymap. The simulated bubbles are compared with the observation data in shape, size and radiation intensity.

The rest of this work is organized as follows. The RMHD model and emission mechanisms are presented in Section 2, simulations on magnetized jets in the Milky Way Galaxy are discussed in Section 3, simulations on the emission from Fermi bubbles and eROSITA bubbles are analyzed in Section 4, and some conclusions are drawn in Section 5.

2. RMHD Model and Emission Mechanisms

Relativistic magnetohydrodynamic (RMHD) equations are chosen to describe the evolution of plasma distribution in the Galaxy, where the plasma velocity in some regions can be sufficiently high, while the effect of general relativity is ignored since the neighborhood of SMBH is not of our main concern. The RMHD equations are given by [16,17,18]

$$\begin{array}{ccc}& & {\partial }_{\mu }\left(\rho u^{\mu }\right)=0\hfill \end{array}$$

(1)

$$\begin{array}{ccc}& & {\partial }_{\nu }{T}^{1\nu }=\rho g_r+\rho {\Phi }_{,1}\hfill \end{array}$$

(2)

$$\begin{array}{ccc}& & {\partial }_{\nu }{T}^{2\nu }=0\hfill \end{array}$$

(3)

$$\begin{array}{ccc}& & {\partial }_{\nu }{T}^{3\nu }=\rho {\Phi }_{,3}\hfill \end{array}$$

(4)

$$\begin{array}{ccc}& & {\partial }_{\nu }{T}^{0\nu }=\frac{1}{c}(\rho v_r{g}_r+\rho v^i{\Phi }_{,i}),i=1,2,3\hfill \end{array}$$

(5)

$$\begin{array}{ccc}& & {\partial }_{\nu }{G}^{\mu \nu }=0\hfill \end{array}$$

(6)

where

$$\begin{array}{ccc}& & u^{\mu }=\frac{d{x}^{\mu }}{d\tau }=\left(\frac{c}{\sqrt{1-v^2/c^2}},\frac{\overline{v}}{\sqrt{1-v^2/c^2}}\right)\hfill \\ & & T^{\mu \nu }=\left(\rho c^2+\epsilon +p+\frac{{\left|b\right|}^2}{4\pi }\right)\frac{u^{\mu }}{c}\frac{u^{\nu }}{c}+\left(p+\frac{{\left|b\right|}^2}{8\pi }\right)g^{\mu \nu }-\frac{b^{\mu }{b}^{\nu }}{4\pi }\hfill \\ & & G^{\mu \nu }=b^{\mu }\frac{u^{\nu }}{c}-b^{\nu }\frac{u^{\mu }}{c}\hfill \\ & & g_r=-\frac{G{M}_{*}}{r^2}\hfill \end{array}$$

are the four-velocity, stress–energy tensor, dual of Faraday tensor and gravitational acceleration, respectively; G is the gravitational constant, $\overline{v}$ is the three-velocity, $d\tau =dt/\gamma$, $\gamma$ is the Lorentz factor, $x^{\mu }=(x^0,x^1,x^2,x^3)=(ct,r,\theta ,z)$ in cylindrical coordinates, $\rho$ is the rest-mass density, $\epsilon$ is the internal energy density, p is the gas pressure, $b^{\mu }$ is the magnetic four-vector, with [18]

$$\begin{array}{ccc}& & b^i=\frac{c{B}^i+u^i{b}^0}{u^0},i=1,2,3\hfill \end{array}$$

(7)

$$\begin{array}{ccc}& & b^0=g_{i\mu }{B}^i\frac{u^{\mu }}{c}\hfill \end{array}$$

(8)

$B^i$ is the magnetic three-vector, and ${\left|b\right|}^2=b^{\mu }{b}^{\nu }{g}_{\mu \nu }=-{\left(b^0\right)}^2+{\left(b^1\right)}^2+{\left(b^2\right)}^2+{\left(b^3\right)}^2$. The gravitation acceleration $g_r$ is attributed to the Sgr ${\mathrm{A}}^{*}$ with mass $M_{*}=4.4\times {10}^6{M}_{\odot }$ [19].

The gravitational field attributed to the mass in the Milky Way Galaxy is given by $\nabla \Phi$, with $\Phi$ the potential function that satisfies the Poisson’s equation, ${\nabla }^2\Phi =4\pi G{\rho }_{\mathrm{bulge}}$, and

$$\begin{array}{ccc}& & {\rho }_{\mathrm{bulge}}={\rho }_b+{\rho }_{\mathrm{thin}}+{\rho }_{\mathrm{thick}}+{\rho }_{\mathrm{HI}}+{\rho }_m+{\rho }_d\hfill \end{array}$$

(9)

where ${\rho }_b$, ${\rho }_{\mathrm{thin}}$, ${\rho }_{\mathrm{thick}}$, ${\rho }_{\mathrm{HI}}$, ${\rho }_m$ and ${\rho }_d$ are the mass densities in the galactic bulge, thin stellar disk, thick stellar disk, $H_I$ (neutral atomic hydrogen) gas disk, molecular disk and dark-matter halo, respectively [7].

The initial temperature of the plasma is assumed to be $T=2.32\times {10}^6$ K [7]. The gas pressure is governed by the ideal gas law, with a mean particle weight of $\mu =0.61$ [7]. The gas pressure is related to the internal energy density as $p=(\Gamma -1)\epsilon$, where $\Gamma =c_p/c_v$ is the ratio of specific heats, which is approximated as $\Gamma =5/3$ in the regions where the plasma flow is mild ($\gamma \simeq 1$).

The initial magnetic-field distribution is assumed to be [20]

$$\begin{array}{ccc}& & {\overline{B}}_p=\frac{1}{r}\nabla \psi \times \widehat{\varphi }\hfill \end{array}$$

(10)

where

$$\begin{array}{ccc}& & \psi =r^2{B}_0{e}^{-r/r_0-z/z_0}\hfill \end{array}$$

(11)

is a flux function modified from the default model in GALPROP, with $|\overline{B}|$ converging at $r=0$ and $\nabla ·\overline{B}=0$ [8,21]; $B_0$ is the field strength, $r_0=2$ kpc and $z_0=2$ kpc.

To ensure numerical stability, all the physical variables are normalized with respect to $L_0=1$ kpc $=3.0856\times {10}^{21}$ cm in length, c in speed, $t_0=c/L_0$ in time, ${\rho }_0={10}^7{M}_{\odot }$/kpc³$=6.8\times {10}^{-25}$ g/cm³ in mass density, $p_0={\rho }_0{c}^2=6.1\times {10}^{-4}$ dyn/cm² in pressure, $T_0=6.6\times {10}^{12}$ K in temperature and $b_0=\sqrt{4\pi {\rho }_0{c}^2}=8.5\times {10}^{-2}$ G in magnetic field. Thus, the four-velocity, stress–energy tensor, dual of Faraday tensor and magnetic four-vector are normalized as

$$\begin{array}{ccc}& & u^{\prime \mu }=\gamma (1,{\overline{v}}^{\prime })\hfill \end{array}$$

(12)

$$\begin{array}{ccc}& & T^{\prime \mu \nu }=({\rho }^{\prime }+{\epsilon }^{\prime }+p^{\prime }+|b^{\prime }{|}^2)u^{\prime \mu }{u}^{\prime \nu }+\left(p^{\prime }+\frac{|b^{\prime }{|}^2}{2}\right)g^{\mu \nu }-b^{\prime \mu }{b}^{\prime \nu }\hfill \end{array}$$

(13)

$$\begin{array}{ccc}& & G^{\prime \mu \nu }=b^{\prime \mu }{u}^{\prime \nu }-b^{\prime \nu }{u}^{\prime \mu }\hfill \end{array}$$

(14)

$$\begin{array}{ccc}& & b^{\prime \mu }=\gamma \left[{\overline{v}}^{\prime }·{\overline{B}}^{\prime },\frac{{\overline{B}}^{\prime }}{{\gamma }^2}+{\overline{v}}^{\prime }({\overline{v}}^{\prime }·{\overline{B}}^{\prime })\right]\hfill \end{array}$$

(15)

$$\begin{array}{ccc}& & |b^{\prime }{|}^2=b^{\prime \mu }{b}_{\mu }^{\prime }=\frac{|{\overline{B}}^{\prime }{|}^2}{{\gamma }^2}+{({\overline{v}}^{\prime }·{\overline{B}}^{\prime })}^2\hfill \end{array}$$

(16)

In the subsequent discussions, the ′ notation is omitted for clarity.

Next, MATLAB codes are developed to solve the RMHD equations by applying a Godunov-type finite-volume algorithm with HLL-type Riemann solvers and slope limiters [22]. A constrained transport method is applied to maintain the divergence-free condition of the magnetic field [23]. The codes are validated in Orszag–Tang and blast wave tests, simulating shock tubes, magnetic flux emergence and relativistic magnetized jets [24,25].

In the simulations, a jet of radius $r_j$ erupts from the center of our Milky Way Galaxy, in both $\widehat{z}$ and $-\widehat{z}$ directions, and lasts for $t_j$. The gas density, thermal energy density, temperature and injection speed are ${\rho }_j$, ${\epsilon }_j$T and $v_j$, respectively. A reflective boundary condition is imposed along the z axis ($r=0$) and the region $r>r_j$ at $z=0$. On the jet base ($r\le r_j$ at $z=0$), an injection boundary condition is imposed during $t<t_j$, and a reflective boundary condition is imposed during $t\ge t_j$. An open boundary condition is imposed on the outer boundaries in the r and z directions. The simulation results with fine grid sizes of $\Delta r=0.01$ and $\Delta z=0.1$ are almost the same as those obtained with coarse grid sizes of $\Delta r=0.02$ and $\Delta z=0.18$.

2.1. Emission Mechanisms of Leptons

The radiation (emission) from the Fermi bubbles and the eROSITA bubbles is governed by the time-independent radiative transfer equation [26]

$$\begin{array}{ccc}& & \frac{\partial I_{\nu }(\overline{r},\overline{\Omega })}{\partial s}=j_{\nu }(\overline{r},\overline{\Omega })-{\kappa }_{\nu }(\overline{r},\overline{\Omega })I_{\nu }(\overline{r},\overline{\Omega })\hfill \end{array}$$

(17)

where $I_{\nu }$ is the radiation intensity, $j_{\nu }$ is the emission coefficient, ${\kappa }_{\nu }$ is the absorption coefficient, s is the radiation path, $\overline{\Omega }=(\partial x/\partial s,\partial y/\partial s,\partial z/\partial s)$ is the direction cosines of radiation intensity and $dt=ds/c$. The radiation is computed by integrating (17) along the line-of-sight (LoS) path, starting from the farther boundary of the computational domain [25], with the result projected onto the skymap. The light traverse effect is neglected because the time scale of plasma evolution is about 10 Myrs, much longer than the light traverse time of several kyrs [25].

The leptonic emission is contributed by synchrotron radiation, inverse Compton scattering and bremsstrahlung [11]. Synchrotron radiation originates from relativistic charged particles accelerated by the magnetic field [27]. The inverse Compton scattering results in net energy transfer from fast moving electrons to photons. The bremsstrahlung is caused by the acceleration of charged particles in the Coulomb field of other charged particles [28]. The interaction of fast-moving electrons with rest protons and ions induces non-thermal bremsstrahlung, while the electron–ion interaction at nonrelativistic temperatures induces thermal bremsstrahlung.

Most electrons in a plasma follow the Maxwell–Boltzmann distribution, while a small fraction of non-thermal electrons follow a power-law distribution in terms of the internal energy as [27,29,30,31]

$$\begin{array}{ccc}& & d{n}_e\left(E\right)=N_0{E}^{-p_e}dE\hfill \end{array}$$

(18)

where $N_0$ is a normalization factor and $p_e=2.2$ is a power-law index. Under a relativistic condition of ${\gamma }_e\gg 1$, $E\simeq {\gamma }_e{m}_e{c}^2$ and (18) is reduced to

$$\begin{array}{ccc}& & d{n}_e\left({\gamma }_e\right)=N\left({\gamma }_e\right)d{\gamma }_e=N_0{\left(m_e{c}^2\right)}^{-p_e+1}{\gamma }_e^{-p_e}d{\gamma }_e\hfill \end{array}$$

(19)

which is integrated to obtain the number density of non-thermal electrons

$$\begin{array}{ccc}& & N_e={\int }_{E_{min}}^{E_{max}}d{n}_e\left(E\right)=\frac{N_0}{p_e-1}(E_{min}^{-p_e+1}-E_{max}^{-p_e+1})=a_{1}n\hfill \end{array}$$

(20)

where $E_{min}$ and $E_{max}$ are the minimum and the maximum internal energies, respectively, where $E_{max}$ is determined by cooling mechanism, n is the total electron number density and $a_1$ is an empirical fraction.

Similarly, the internal energy density of non-thermal electrons is given by

$$\begin{array}{ccc}& & U_e={\int }_{E_{min}}^{E_{max}}Ed{n}_e\left(E\right)=\frac{N_0}{p_e-2}(E_{min}^{-p_e+2}-E_{max}^{-p_e+2})=a_2\epsilon \hfill \end{array}$$

(21)

where $\epsilon =p/(\Gamma -1)$ is the total internal energy density and $a_2$ is another empirical fraction.

In [32], the time evolution of bubbles was simulated by using an MHD model, along with the cosmic-ray (CR) spectrum. The electron energy of cosmic rays achieved maximum near the front edge of the bubble, at $t=1.2$ Myrs, which was about hundreds of GeV at low galactic latitudes and about 1 TeV at high latitudes. Initially, the maximum energy was $E_{max}\simeq {10}^4$ GeV; then, it decayed rapidly within ∼0.3 Myrs due to synchrotron and inverse-Compton cooling, and finally settled to several hundreds of GeV, compatible with the observed value of $E_{max}\simeq 250$ GeV in the present day.

In this work, both non-thermal electrons and protons are assumed to follow the same power-law distribution over the energy internal of $[E_{min},E_{max}]$. Based on the simulation and observation data in [32], we set $E_{max}=250$ GeV. Then, $E_{min}$ is determined by dividing (21) with (20) to obtain

$$\begin{array}{ccc}& & \frac{U_e}{N_e}=\frac{E_{min}^{-p_e+2}-E_{max}^{-p_e+2}}{E_{min}^{-p_e+1}-E_{max}^{-p_e+1}}\frac{p_e-1}{p_e-2}\hfill \end{array}$$

(22)

which is solved for $E_{min}$ by using a secant method [33]. The normalization factor $N_0$ is determined by substituting $E_{min}$ into (20).

The spectral emissivity of leptons is given by [10,11]

$$\begin{array}{ccc}& & j_{\nu }\left(ϵ\right)=j_{\mathrm{sy}}\left(ϵ\right)+j_{\mathrm{ic}}\left(ϵ\right)+j_{\mathrm{ff}}^{\mathrm{nt}}\left(ϵ\right)+j_{\mathrm{ff}}^{\mathrm{t}}\left(ϵ\right)\hfill \end{array}$$

(23)

where $ϵ=h\nu$ is the photon (radiation) energy, and $j_{\mathrm{sy}}\left(ϵ\right)$, $j_{\mathrm{ic}}\left(ϵ\right)$, $j_{\mathrm{ff}}^{\mathrm{nt}}\left(ϵ\right)$ and $j_{\mathrm{ff}}^{\mathrm{t}}\left(ϵ\right)$ are the spectral emissivities of synchrotron radiation, inverse Compton scattering, non-thermal bremsstrahlung and thermal bremsstrahlung, respectively.

The emissivity of synchrotron radiation is given by [11,34]

$$\begin{array}{ccc}& & j_{\mathrm{sy}}\left(ϵ\right)=\frac{4\sqrt{3}\pi r_e}{9}{m}_{e}c\frac{{\nu }^2}{{\nu }_B}{\int }_{{\gamma }_{min}}^{{\gamma }_{max}}\frac{d{\gamma }_e}{{\gamma }_e^4}\frac{d{n}_e}{d{\gamma }_e}F\left(x\right)\hfill \end{array}$$

(24)

where $x=\nu /\left(3{\nu }_B{\gamma }_e^2\right)$, ${\nu }_B=q_{e}B/\left(2\pi m_{e}c\right)$ and

$$\begin{array}{ccc}& & F\left(x\right)=K_{4/3}\left(x\right)K_{1/3}\left(x\right)-\frac{3}{5}x\left[K_{4/3}^2\left(x\right)-K_{1/3}^2\left(x\right)\right]\hfill \end{array}$$

(25)

$K_m\left(x\right)$ is the mth-order modified Bessel function of the second kind.

The spectral emissivity of inverse Compton scattering is given by [28]

$$\begin{array}{ccc}& & j_{\mathrm{ic}}\left(ϵ\right)=3c{\sigma }_Tϵ{\int }_{{\gamma }_{emin}}^{{\gamma }_{emax}}\int \frac{1}{4{\gamma }_e^2{ϵ}_0}N\left({\gamma }_e\right)n_{\mathrm{ph}}\left(\frac{{ϵ}_0}{m_e{c}^2}\right)f(x_{\mathrm{ic}},{ϵ}_0)d{ϵ}_{0}d{\gamma }_e\hfill \end{array}$$

(26)

where ${ϵ}_0$ is the energy of incident photon, $n_{\mathrm{ph}}$ is the number density of the cosmological microwave background (CMB) and interstellar radiation field (ISRF), $x_{\mathrm{ic}}=ϵ/\left(4{\gamma }_e^2{ϵ}_0\right)$, ${\sigma }_T=8\pi r_e^2/3$ is the Thomson cross-section, $N\left({\gamma }_e\right)$ is the distribution of electrons and

$$\begin{array}{ccc}& & f(x_{\mathrm{ic}},{ϵ}_0)=2x_{\mathrm{ic}}ln{x}_{\mathrm{ic}}+x_{\mathrm{ic}}+1-2x_{\mathrm{ic}}^2\hfill \end{array}$$

Since non-thermal electrons follow the power-law distribution, (26) is transformed to

$$\begin{array}{ccc}& & j_{\mathrm{ic}}\left(ϵ\right)=\frac{3}{4\pi }c{\sigma }_{T}C{2}^{p_e-2}{ϵ}^{-(p_e-1)/2}{\int }_{x_{\mathrm{ic}min}}^{x_{\mathrm{ic}max}}\int {ϵ}_0^{(p_e-1)/2}\hfill \\ & & n_{\mathrm{ph}}\left(\frac{{ϵ}_0}{m_e{c}^2}\right)x_{\mathrm{ic}}^{(p_e-1)/2}f(x_{\mathrm{ic}},{ϵ}_0)d{ϵ}_{0}d{x}_{\mathrm{ic}}\hfill \end{array}$$

(27)

where $C=N_0{\left(m_e{c}^2\right)}^{-p_e+1}$, $x_{\mathrm{ic}min}=ϵ/\left(4{\gamma }_{max}^2{ϵ}_0\right)$ and $x_{\mathrm{ic}max}=ϵ/\left(4{\gamma }_{min}^2{ϵ}_0\right)$. The spectral emissivities of non-thermal and thermal bremsstrahlung can be found in [11,34], respectively.

The absorption coefficient $\kappa \left(\nu \right)$ in the observer coordinates is transformed from ${\kappa }^{\prime }\left({\nu }^{\prime }\right)$ in the comoving coordinates of plasma as follows [35]:

$$\begin{array}{ccc}& & \kappa \left(\nu \right)={\delta }^{-1}{\kappa }^{\prime }\left({\nu }^{\prime }\right)\hfill \end{array}$$

(28)

where

$$\begin{array}{ccc}& & {\kappa }^{\prime }\left({\nu }^{\prime }\right)=-\frac{3\sqrt{3}{\sigma }_{T}c{U}_B{c}_1}{16{\pi }^2{m}_e{\nu }^{\prime 2}{\nu }_B}[h({\nu }^{\prime },x_1,p_e,K)-h({\nu }^{\prime },x_2,p_e,K)]\hfill \end{array}$$

(29)

is the electron self-absorption coefficient, $\delta =1/\left[\gamma (1-vcos\theta /c)\right]=\nu /{\nu }^{\prime }$ is the Doppler factor, $\theta$ is the angle between the plasma flow direction and the line-of-sight direction,

$$\begin{array}{ccc}& & h(\nu ,x,p_e,K)=(p_e+2)3^{p_e/2}{\left(\frac{{\nu }^{\prime }}{{\nu }_B}\right)}^{-p_e/2}{x}^{c_2+p_e/2}{\left(c_{3}x\right)}^{-c_2-p_e/2}\Gamma (c_2+p_e/2,c_{3}x)\hfill \end{array}$$

$\Gamma (x,y)$ is the incomplete gamma function, $U_B={\left|\overline{B}\right|}^2/\left(8\pi \right)$ is the magnetic energy density, $p_e=2.2$, $c_1=0.78$, $c_2=0.25$, $c_3=2.175$, $K=N_0{\left(m_e{c}^2\right)}^{-p_e+1}$, $x_1={\nu }^{\prime }/\left(3{\gamma }_{max}^2{\nu }_B\right)$, $x_2={\nu }^{\prime }/\left(3{\gamma }_{min}^2{\nu }_B\right)$ and ${\nu }_B=q_{e}B/\left(2\pi m_{e}c\right)$. The absorption of photons by electrons can be neglected if the photon frequency is sufficiently high.

In short, the distribution $n_e\left(E\right)$ of non-thermal electrons and the associated parameters $N_0$ and $E_{min}$ are determined first. Then, the emissivities and absorption coefficients are computed by using (23) and (28), respectively. Finally, the radiation intensity on the skymap is computed by solving (17).

2.2. Emission Mechanisms of Hadrons

The emission from Fermi bubbles and eROSITA bubbles may be contributed by an inelastic proton–proton collision that produces $\gamma$-rays and secondary particles [10,36,37]. A proton–proton collision may lead to $p+p\to p+{\Delta }^{+}$ or $p+p\to n+{\Delta }^{++}$ [10,36], where ${\Delta }^{+}$ ($uud$) and ${\Delta }^{++}$ ($uuu$) stand for Delta baryons, which decay rapidly on a time scale of $5.63\times {10}^{-24}$ s to produce charged and neutral pions. In other words, pions are produced alternatively via $p+p\to p+p+{\pi }^0$, $p+p\to p+p+2{\pi }^0$, $p+p\to p+n+{\pi }^{+}+{\pi }^0$, $p+p\to D+{\pi }^{+}+{\pi }^0$ or $p+p\to p+p+{\pi }^{+}+{\pi }^{-}$[36], where D stands for deuterium. The fractions of pions thus produced were $({\pi }^{+},{\pi }^{-},{\pi }^0)=(0.6,0.1,0.3)$ at 1 GeV [10]. The charged pions decay to produce secondary electrons, and the neutral pions decay to emit $\gamma$-rays via ${\pi }^0\to \gamma +\gamma$ [37].

The spectral emissivity of $\gamma$-rays attributed to the decay of neutral pions is given by [10,36]

$$\begin{array}{ccc}& & j_{\pi \gamma }\left(ϵ\right)=ϵ{\lambda }_c{m}_e{c}^2{n}_H{\int }_{{\gamma }_{pmin}}^{{\gamma }_{pmax}}\frac{d\sigma }{dϵ}{N}_p\left({\gamma }_p\right)d{\gamma }_p\hfill \end{array}$$

(30)

with the differential cross-section of $\gamma$-rays via the process ${\pi }^0\to \gamma +\gamma$ given by [36]

$$\begin{array}{ccc}& & \frac{d\sigma ({\gamma }_p,ϵ)}{dϵ}=P\left({\gamma }_p\right)F({\gamma }_p,ϵ)\hfill \end{array}$$

(31)

where $P\left({\gamma }_p\right)$ is a peak function and $F({\gamma }_p,ϵ)$ is a spectral shape function. The number density $N_p\left({\gamma }_p\right)$ of high-energy protons at Lorentz factor ${\gamma }_p$ is related to the plasma number density n as $N_p=a_{1p}n$, and the high-energy proton internal energy density $U_p$ is related to the total internal energy density $\epsilon$ as $U_p=a_{2p}\epsilon$. In this work, we make the approximation that $a_{1p}=a_1$ and $a_{2p}=a_2$.

The steady-state distribution of secondary electrons is approximated as [10,37]

$$\begin{array}{ccc}& & N_e\left(E_e\right)=t_{\mathrm{cool}}\left(E_e\right){\int }_{E_e}^{E_e^{max}}{Q}_e\left(E_e^{\prime }\right)d{E}_e^{\prime }\hfill \end{array}$$

(32)

where

$$\begin{array}{ccc}& & t_{\mathrm{cool}}\left(E_e\right)={\left(\frac{1}{{\tau }_{\mathrm{Comb}}}+\frac{1}{{\tau }_{\mathrm{brem}}}+\frac{1}{{\tau }_{\mathrm{sync}+\mathrm{IC}}}\right)}^{-1}\left(\mathrm{yr}\right)\hfill \end{array}$$

(33)

is the effective electron cooling time scale [38], with the contributions from Coulomb, bremsstrahlung, and synchrotron plus inverse Compton scattering given by [38]

$$\begin{array}{ccc}& & {\tau }_{\mathrm{Comb}}=2.53\times {10}^9\left(\frac{|\overline{p}|}{100m_{e}c}\right){\left(\frac{n}{{10}^{-3}}\right)}^{-1}\left(\mathrm{yr}\right)\hfill \\ & & {\tau }_{\mathrm{brem}}=6\times {10}^{12}{\left(\frac{n}{{10}^{-3}}\right)}^{-1}\left(\mathrm{yr}\right)\hfill \\ & & {\tau }_{\mathrm{sync}+\mathrm{IC}}=2.3\times {10}^8{\left(\frac{|\overline{p}|}{{10}^4{m}_{e}c}\right)}^{-1}{\left(1+\frac{u_B}{u_{\mathrm{cmb}}}\right)}^{-1}\left(\mathrm{yr}\right)\hfill \end{array}$$

respectively, where $|\overline{p}|$ is the electron momentum, and $u_B$ and $u_{\mathrm{cmb}}$ are the energy densities of magnetic field and cosmic microwave background, respectively. The Coulomb interaction dominates low-energy electrons, synchrotron radiation and inverse Compton scattering mainly affecting high-energy electrons, while the bremsstrahlung radiation is negligible.

The energy spectrum of secondary electrons is given by [37]

$$\begin{array}{ccc}& & Q_e\left(E_e\right)=2{\int }_0^1{f}_e\left(x_e\right)J_{\pi }(E_e/x_e)\frac{d{x}_e}{x_e}\hfill \end{array}$$

(34)

where $f_e\left(x_e\right)$ is the distribution function of electrons, $x_e=E_e/E_{\pi }$ and $J_{\pi }(E_e/x_e)$ is the production rate of pions, which is given by [37]

$$\begin{array}{ccc}& & J_{\pi }\left(E_e/x_e\right)\frac{d{x}_e}{x_e}=-J_{\pi }\left(E_{\pi }\right)\frac{d{E}_{\pi }}{E_{\pi }}=-F_{\pi }(x,E_p)\frac{dx}{E_{\pi }}\hfill \end{array}$$

(35)

where $E_p$ is the proton energy, $x=E_{\pi }/E_p$ and $F_{\pi }(x,E_p)$ is the energy spectrum of pions.

By substituting (35) into (34), the energy spectrum of secondary electrons is reduced to

$$\begin{array}{ccc}& & Q_e\left(E_e\right)=2{\int }_{E_e}^{E_{pmax}}{f}_e\left(E_e/E_{\pi }\right)F_{\pi }(E_{\pi }/E_p,E_p)\frac{d{E}_{\pi }}{E_{\pi }{E}_p}\hfill \end{array}$$

(36)

which is weighted by the number density distribution of protons, $N_p\left(E_p\right)$, to become

$$\begin{array}{ccc}& & Q_e\left(E_e\right)=2{\int }_{E_{pmin}}^{E_{pmax}}{\int }_{E_e}^{E_{pmax}}{f}_e\left(E_e/E_{\pi }\right)F_{\pi }(E_{\pi }/E_p,E_p)\frac{N_p\left(E_p\right)}{E_{\pi }{E}_p}d{E}_{\pi }d{E}_p\hfill \end{array}$$

(37)

In short, the emission from hadrons is obtained by first computing the spectral emissivity of $\gamma$-rays attributed to the decay of neutral pions with (30). Then, the steady-state distribution of secondary electrons is computed with (32), and the radiation intensity is obtained by solving (17).

3. Simulations on Magnetized Jets in Milky Way Galaxy

In this section, we will simulate jet eruption from the center of our Milky Way Galaxy with the RMHD model.

Table 1 lists the simulation cases and relevant parameters in this work and the literature. Only Fermi bubbles were simulated in [6,7], while both Fermi and eROSITA bubbles were simulated in [8]. The evolution of Fermi and eROSITA bubbles were simulated with the HD or MHD models [6,7,8]. The magnetic field plays an important role in the jet evolution, but was not included in an HD model [25].

Debates on jet eruptions and their aftermaths on large-scale structures like bubbles have not been settled. The evolution time of bubbles after jet eruption was conjectured to be from 1.2 Myrs to 9 Myrs, accompanied by different sets of simulation parameters. For example, the inclusion of initial mass distribution in the Galaxy may prolong the formulation of Fermi and eROSITA bubbles.

The jet parameters, $t_j$, $r_j$, ${\rho }_j$, ${\epsilon }_j$ and $v_j$, cannot be directly acquired by observation. Only the radiation intensity from the bubbles and the galactic disk are observed. Typically, the jet parameters were selected and justified by matching the simulation results with the observation data. Different sets of parameters have been proposed to account for similar observation features of Fermi and eROSITA bubbles, possibly due to the complicated and nonlinear physical processes underpinning these phenomena.

Table 2 lists the simulation parameters adopted in this work. Run A is designed to validate the RMHD model by comparing the simulation results with the literature [8]. The parameters in run B and run C are the same as in run A, except that the initial mass distribution in the Galaxy is considered in run B, and the magnetic field is absent in run C.

The parameters in run D are selected to form another benchmark for comparison. The initial mass distribution in the Galaxy and the magnetic field remain the same as in run B. Compared with runs A, B and C, the duration of jet injection ($t_j$) is significantly extended, the jet radius ($r_j$) is slightly reduced, and the mass density (${\rho }_j$) and the internal energy density (${ϵ}_j$) of the jet are also reduced.

The parameters in runs E, F and G are the same as in run D, except the mass density is increased to ${\rho }_j=1.95\times {10}^{-26}$ (g/cm³) in run E and run G, while the internal energy density is decreased to ${ϵ}_j=1.29\times {10}^{-11}$ (erg/cm³) in run F and run G.

Table 3 lists the composition of ejected power and total ejected energy, based on the jet parameters listed in Table 2. The ejected power is the product of the energy density, cross-section and speed of the jet. The total ejected energy is the product of the ejected power and duration of the jet injection.

Figure 1 shows the initial and evolutional distributions of mass density, gas pressure/ velocity and temperature in our Milky Way Galaxy, with the same simulation parameters as in [8]. Figure 1a manifests a disk-like distribution of mass density around the Galactic Equator. The mass density distribution, as given in (9), is composed of a galactic bulge, thin stellar disk, thick stellar disk, $H_I$ gas disk, molecular disk and dark-matter halo. The galactic bulge extends from the galactic center to $r\simeq 5$ and $z\simeq 2.5$. The mass density of thin and thick stellar disks decay in the z direction. The thicknesses of the $H_I$ gas disk and molecular gas disk are about 0.5 and 0.2, respectively. The dark-matter halo permeates the whole computational domain, and is denser around the galactic center.

Figure 1b,c show the initial distributions of gas pressure and temperature, respectively. The gas pressure is proportional to $\rho$, under the ideal gas law. The initial temperature is uniform with $T=2.32\times {10}^6$ (K) [7].

Figure 1d shows that the plasma jet induces a forward shock to compress the ambient medium, forging a shell that reaches $z\simeq 11$ and $r\simeq 7$, respectively. The jet induces a high-mass region centered at $z\simeq 7.5$, with $\rho \simeq {10}^{-27}$ (g/cm³). A low-mass region in $2<z<4$ near the z axis follows the high-mass region, with $\rho \simeq {10}^{-29}$ (g/cm³).

Figure 1e shows that the plasma flows outwards, with a maximum speed of about $3\times {10}^8$ (cm/s) around the forward shock, and decreases to about ${10}^7$ (cm/s) behind the forward shock. Figure 1f shows that the temperature around the forward shock reaches about ${10}^8$ K. A low-temperature region is observed near the z axis around $z\simeq 5$, with $T\simeq {10}^5$ K. A similar low-temperature region was also observed in [8].

In short, the simulated bubbles with the proposed RMHD model in run A match well with the literature in shape, size, large-scale structure, velocity distribution and evolution time.

Figure 2 shows the distributions of mass density and temperature in run B and run C, respectively, at $t=2.6$ Myrs. With run A as the benchmark, run B includes the initial mass distribution in the Galaxy, and run C excludes the magnetic field.

Figure 2a shows that the mass density within the bubble ranges from ${10}^{-28}$ to ${10}^{-27}$ (g/cm³), which is higher than that in run A. The mass density around the forward shock is higher than that in the central part of the bubble. The forward shock travels to $z\simeq 9$ and $r\simeq 5$, respectively. The forward shock expands slightly slower than that in run A, impeded by the initial mass distribution in the Galaxy. In other words, the bubble takes a longer time to evolve to the same size as in run A when the initial mass distribution is present.

Figure 2b manifests a discontinuity in temperature distribution around the forward shock. The spread of the high-temperature region near the Galactic Equator is obviously impeded as compared with run A.

Figure 2c manifests a high-mass region dragged by the jet, followed by a low-mass region. The features in $z<4$ and $r<4$ become simpler as compared with run A. Figure 2d manifests a discontinuity in temperature distribution around the forward shock. No complex features are observed in $z<4$ and $r<4$, as in Figure 2c. The spread of the high-temperature region near the Galactic Equator is obviously impeded as compared with run A, but moves farther than in run B.

In order to reconstruct both Fermi bubbles and eROSITA bubbles from a single jet eruption, the jet parameters in run D are adjusted to match the observed morphology, as shown in Figure 3. The duration of jet injection ($t_j$) is significantly extended to induce a hot region in the central part, reconstructing the Fermi bubbles that emit $\gamma$-rays. The jet radius $r_j$ is slightly reduced to control the size of the bubbles.

With run D as benchmark, the mass density is increased to ${\rho }_j=1.95\times {10}^{-26}$ (g/cm³) in run E and run G, while the internal energy density is decreased to ${ϵ}_j=1.29\times {10}^{-11}$ (erg/cm³) in run F and run G, to investigate their effects on bubble formation.

Figure 3 shows the distributions of mass density, velocity and temperature of the plasma in run D, at $t=5.15$ Myrs, which roughly reveal the morphology of Fermi bubbles and eROSITA bubbles. All the three distributions manifest an outer shell enclosed by a forward shock (white curve) and a central lobe enclosed by a discontinuity transition (red curve). The inner part of the central lobe, featuring low mass density and low temperature, is enclosed by a yellow curve, which will be used to computed the $\gamma$-ray residual map later.

The outer shell expands to $r\simeq 7$ and $z\simeq 13$. Its shape and size are comparable to those of eROSITA bubbles [3]. Compared with the central lobe, the outer shell has a higher density of $\rho \simeq {10}^{-28}$–${10}^{-27}$ (g/cm³) and a lower temperature of ${10}^7$–$7\times {10}^7$ (K).

The central lobe extends to $r\simeq 4$ and $z\simeq 10$, respectively. Its shape and size are comparable to those of Fermi bubbles [3]. The mass density is about $\rho ={10}^{-29}$ (g/cm³). The outer edge of the central lobe, with a thickness of about 2.5 (kpc), manifests a high temperature of about ${10}^8$–$2\times {10}^8$ (K). The temperature off the outer edge decreases to about ${10}^7$ (K).

Figure 3b shows a boundary segment (indicated by a bold black arrow) between the central lobe and the outer shell, across which the plasma motion is parallel to the boundary; the mass density, speed and temperature manifest significant discontinuity. The discontinuity across the other boundary segment is not so obvious, as the plasma flows across the boundary. The tangential discontinuity and the forward shock mark the edges of the Fermi bubbles and eROSITA bubbles, respectively, brought by the same jet eruption.

The plasma is composed of both thermal and non-thermal particles, with the latter much fewer in number than the former. In the outer edge of the central lobe, the rising temperature increases the number of relativistic particles, which are responsible for $\gamma$-ray emission via inverse Compton scattering and non-thermal bremsstrahlung.

Figure 4 shows the distributions of mass density and temperature in runs E, F and G, respectively, at $t=5.15$ Myrs, which are compared with their counterparts in run D.

Referring to Table 2 and Table 3, the jet mass density and the ejected kinetic power in run E are ten times larger than their counterparts in run D, and the ejected thermal power is the same in both run D and run E.

In run E, the mass density in the central lobe increases to about ${10}^{-27}$, and the temperature in the central lobe decreases to about ${10}^6$–${10}^7$ K. A low-temperature region appears near the z axis behind the discontinuity transition, similar to Figure 1f and Figure 2d, because larger kinetic energy implies lower thermal energy, hence lower temperature.

In run F, the jet mass density and the ejected kinetic power are the same as those in run D, and the ejected thermal power is ten times smaller than that in run D. The forward shock reaches $r\simeq 5$ and $z\simeq 8$, respectively, and the bubbles are significantly smaller in size than those in run D.

The resulting mass density distribution manifests an outer shell with higher density and a central lobe with lower density. The temperature in the front edge of the central lobe is higher than that behind it. In short, the distributions appear similar to their counterparts in run D, but it takes much longer evolution time in run F to achieve a similar size to that in run D because of the lower ejected thermal energy in run F.

In run G, the jet mass density and the ejected kinetic power are the same as those in run E, but the ejected thermal power is ten times smaller than that in run E. The resulting mass distribution appears very similar to that in run E, but the size of bubbles is slightly smaller. The temperature in the central lobe is lower than that in run E because of the lower ejected thermal energy in run G.

In summary, the initial mass distribution in the Galaxy tends to impede the expansion of bubbles. The features internal to the bubbles become more complicated if the magnetic field is present. If the jet mass density and the ejected kinetic power increase while the ejected thermal power remains the same, the resulting mass density increases and the temperature decreases in the central lobe. If the ejected thermal power decreases while the jet mass density and the ejected kinetic power remain the same, the bubbles evolve more slowly. A jet carrying large kinetic energy tends to induce a depletion region with low temperature behind the discontinuity transition.

Spiral-Arm Model

So far, the simulation results have manifested salient features of Fermi bubbles and eROSITA bubbles. A spiral-arm model of mass density is adopted to substitute the disk model near the Galactic Equator, including more details of mass distribution in the Galaxy. In our azimuthal symmetric model, the stellar disk is approximated as concentric tubes, with the cross-sectional sizes compatible with the observed pattern in the spiral arms [39]. Explicitly, the mass density in the spiral arms is given by

$$\begin{array}{ccc}& & {\rho }_{\mathrm{thin}}=\frac{{\Sigma }_{\mathrm{thin}}^{\prime }}{2z_{\mathrm{thin}}}{e}^{-\left|z\right|/z_{\mathrm{thin}}-r/r_{\mathrm{thin}}}\hfill \end{array}$$

$$\begin{array}{ccc}& & \left[e^{-{(r/2.2)}^2}+e^{-{((r-5.2)/0.8)}^2}+e^{-{((r-6.9)/0.4)}^2}+e^{-{((r-8.85)/0.55)}^2}\right]\hfill \end{array}$$

(38)

$$\begin{array}{ccc}& & {\rho }_{\mathrm{thick}}=\frac{{\Sigma }_{\mathrm{thick}}^{\prime }}{2z_{\mathrm{thick}}}{e}^{-\left|z\right|/z_{\mathrm{thick}}-r/r_{\mathrm{thick}}}\hfill \end{array}$$

$$\begin{array}{ccc}& & \left[e^{-{(r/2.2)}^2}+e^{-{((r-5.2)/0.8)}^2}+e^{-{((r-6.9)/0.4)}^2}+e^{-{((r-8.85)/0.55)}^2}\right]\hfill \end{array}$$

(39)

where ${\Sigma }_{\mathrm{thin}}^{\prime }$ and ${\Sigma }_{\mathrm{thick}}^{\prime }$ are adjusted to make the total mass in the spiral-arm model the same as that in the disk model.

Figure 5 shows the distributions of mass density in the thin and thick spiral arms, respectively, as well as their sum in the spiral-arm model. The mass density distribution is composed of a central disk with radius 2.2 in the r direction, a first tube of width 0.8 centered at $r=5.2$, a second tube of width 0.4 centered at $r=6.9$ and a third tube of width 0.55 centered at $r=8.85$. The central disk and the three tubes are concentric about the z axis.

Figure 6 shows the distributions of mass density and temperature in run D, with the spiral-arm model, at $t=5$ Myrs. The central lobe extends to $r\simeq 4$ and $z\simeq 10$, and the outer shell extends to $r\simeq 7$ and $z\simeq 12$. The outer shell has a higher mass density and the central lobe has a lower mass density. The temperature in the outer shell is about ${10}^7$–$8\times {10}^7$ (K), and that in the outer edge of the central lobe, with a thickness of about $2.5$ (kpc), is about ${10}^8$–$2\times {10}^8$ (K), which decreases to about ${10}^7$ (K) in the inner part.

In general, the shape, size and internal features are very similar to those shown in Figure 3, simulated with the galactic disk model. The mass density in the outer shell with the spiral-arm model is a little lower than that with the galactic disk model, and the temperature in the outer shell with the spiral-arm model is a little higher than that with the galactic disk model.

4. Simulations on Emission from Fermi Bubbles and eROSITA Bubbles

Figure 7 shows the simulation scenario and the coordinate system used to compute $\gamma$-ray and X-ray emission by solving the radiative transfer Equation (17). The computational domain is a cylindrical region which is axisymmetric about the z axis. The Sgr ${\mathrm{A}}^{*}$ is located at the origin and our Sun lies on the x axis, at a distance of 8.32 kpc from the origin [40]. The spectral radiation intensity $I_{\nu }$ observed on the Earth is computed by solving (17) and integrating along the line-of-sight coordinate s, marked by red line, shown in Figure 7a. Figure 7b,c show the side view and top view, respectively, of the scenario. The bubble is bounded in elevation by galactic latitude $b=b_0$ and in azimuth by galactic longitude $\ell ={\ell }_0$.

Figure 8a shows the probability density functions of thermal electrons and non-thermal electrons. The thermal electrons follow the Maxwell–Boltzmann distribution at $T={10}^7$ K, a typical temperature in Fermi bubbles [41]. The energy of most thermal electrons falls in $E<{10}^4$ eV ($\gamma =1.019$). Non-thermal electrons follow a power-law distribution

$$\begin{array}{ccc}& & f_n\left(E\right)=\frac{N_0}{N_e}{E}^{-p_e}\hfill \end{array}$$

(40)

The probability density at $E={10}^5$ eV is about $f\left(E\right)={10}^{-5}$. The minimum energy $E_{min}$ of non-thermal electrons, as marked in Figure 8a, is determined by solving (22).

Figure 8b shows the spectrum of CMB plus ISRF, which are the photons induced by inverse Compton scattering in the Milky Way Galaxy. The spectra marked by the red curve and black curve are computed at $(r,z)=(0,0)$ and $(8,0)$, respectively, via GALPROP [42] and a radiation transfer (RT) model [43], respectively, which were well calibrated for use in our Milky Way Galaxy [42,43].

The ISRF spectral profile in the computational domain is derived by interpolating in the r direction from the available data at $(r,z)=(0,0),(4,0),(12,0)$ in [42] and those at $(1,0)$ and $(8,0)$ in [43]. Since these five profiles are similar to one another, only those at $(r,z)=(0,0)$ and $(8,0)$ are shown in Figure 8b. Note that the contribution of emission from the region $r>8$ is negligible.

The ISRF spectrum varies from place to place in the Milky Way Galaxy. In general, the ISRF spectrum is stronger near the galactic center and decreases outwards [42]. The spectrum at $(r,z)=(8,0)$ (kpc) is roughly comparable to that at $(r,z)=(0,0)$ (kpc) in $3\times {10}^{-2}<{ϵ}_0<{10}^{-1}$ (eV). In ${10}^{-1}<{ϵ}_0<{10}^1$ (eV), the spectrum at $(r,z)=(8,0)$ (kpc) is about one order smaller than that at $(r,z)=(0,0)$ (kpc). At $(r,z)=(8,0)$, the maximum number density of photons is $n_{\mathrm{ph}}\simeq 5\times {10}^3$ (1/cm³/eV) at ${ϵ}_0=6\times {10}^{-3}$ eV, and decreases to $n_{\mathrm{ph}}\simeq 2\times {10}^{-4}$ (1/cm³/eV) at ${ϵ}_0=10$ eV.

Similarly, the ISRF spectral profile in the z direction is interpolated from the available data at $\lambda =24,62,100\mathsf{\mu }$m in [43] and those at $\lambda =2.2,340\mathsf{\mu }$m in [44], as shown in Figure 8c. It is observed that the ISRF ($n_{\mathrm{ph}}$) approximately follows the power-law distribution with respect to z as $n_{\mathrm{ph}}(\lambda ,z)=n_{\mathrm{ph}}(\lambda ,0)\times z^{-\left|L\right|}$. By regressing with the data in [43,44], the power-law indices at $\lambda =2.2,24,62,100,340\mathsf{\mu }$m are estimated as $L=-1.2,-0.89,-0.79,-0.90$ and $-0.47$, respectively. The profiles of $n_{\mathrm{ph}}(\lambda ,z)$ at other wavelengths are interpolated in terms of these five power-law profiles.

Figure 8d,e show the distributions of $j_{\mathrm{ic}}$ attributed to inverse Compton scattering, computed with the ISRF spectrum at $(r,z)=(8,0)$ via the RT model [43] and the interpolated ISRF spectrum, respectively. Each distribution manifests an outer shell with a higher magnitude and a central lobe with a lower magnitude. The emissivity around the origin and the outer shell near the Galactic Equator in Figure 8d is weaker than that in Figure 8e. The emissivity around the forward shock in Figure 8d is slightly stronger than that in Figure 8e. In this work, the interpolated ISRF spectrum is adopted in the subsequent simulations.

Figure 8f shows the spectral emissivity of $j_{\mathrm{ic}}$, with $N_0=1$. The magnitude of $j_{\mathrm{ic}}\left(\nu \right)$ decays exponentially with the photon frequency. Its value is $j_{\mathrm{ic}}\simeq {10}^{-7}$ (erg/cm³/sr/s/erg) at $\nu ={10}^{15}$ Hz (4.1 eV), and $j_{\mathrm{ic}}\simeq {10}^{-13}$ (erg/cm³/sr/s/erg) at $\nu ={10}^{25}$ Hz (41 GeV).

Figure 9a shows the distribution of minimum internal energy $E_{min}$ (eV) in run D, which is about ${10}^4$ eV in the outer edge of the central lobe, and about ${10}^3$ behind the forward shock. The value of $E_{min}$ is determined by solving the self-consistent Equations in (20)–(22), given the empirical values of $a_1$, $a_2$ and $E_{max}$. These three equations are based on the power-law distribution of non-thermal electrons, which may lose their energy through Coulomb interaction, bremsstrahlung, synchrotron radiation and inverse Compton scattering. The Coulomb interaction dominates low-energy electrons, and the cooling time scale at low energy is typically shorter than the time scale of plasma evolution. The cooling of non-thermal electrons has a negligible effect on the thermal plasma due to the tiny fraction of non-thermal electrons.

Figure 9b shows the spectral emissivity $d{j}_{\mathrm{ic}}/dE$ attributed to inverse Compton scattering, with $N_0=1$. It is observed that the electrons with energy of about ${10}^6$ (eV) are responsible for X-ray emission in 0.6–1 (keV), and those with energy higher than ${10}^9$ (eV) are responsible for $\gamma$-ray emission in 2–5 GeV. Non-thermal electrons with power-law energy spectrum are the dominant source of emission in this work. In situ acceleration and cooling effect may affect the electron energy spectrum [9,38], leading to subtle differences in radiation.

The time scales of cooling due to Coulomb loss on non-thermal electrons at $E={10}^6$ and ${10}^9$ (eV) are estimated as ${10}^4$ and ${10}^7$ years, respectively. The evolution time scale of Fermi bubbles and eROSITA bubbles in the simulations is about ${10}^6$–${10}^7$ yrs. By considering the cooling effect, the power-law spectrum of non-thermal electron number density in (18) is modified as [45]

$$\begin{array}{ccc}& & d{n}_e\left(E\right)=N\left(E\right)dE=N_0{\left(E^{3/2}+1.16\times {10}^{-5}{\lambda }_{\mathrm{ee}}nt\right)}^{-(2p_e+1)/3}{E}^{1/2}dE\hfill \end{array}$$

(41)

where t (in sec) is the time since jet eruption,

$$\begin{array}{ccc}& & {\lambda }_{\mathrm{ee}}=30.9-ln\left[n^{1/2}\left(\frac{1}{k_B{T}_e}\right)\right]\hfill \end{array}$$

is an empirical formula of Coulomb logarithm and $k_B{T}_e$ is in unit of keV. Note that (41) reduces to (18) at $t=0$.

Figure 9c shows the spectra of non-thermal electron number density in (41), under cooling effect, with $N_0=1$, $n=6\times {10}^{-5}$ (1/cm³) and $T_e={10}^8$ (K). Due to the cooling effect, the number density $N\left(E\right)$ decreases with time if $E<{10}^5$ (eV), but changes slightly if $E>{10}^5$ (eV). The emission of X-rays and $\gamma$-rays via inverse Compton scattering is dominated by non-thermal electrons with $E>{10}^5$ (eV), of which the number density barely changes by the cooling effect in $t\le 5$ Myrs. Thus, reconstructions of the bubbles in X-rays and $\gamma$-rays, respectively, on the skymap with the spectrum of non-thermal electron number density in (41) are indiscernible from their counterparts reconstructed with the power-law spectrum in (18).

The image of Fermi bubbles was acquired from the received photons with energy of 1–50 GeV [2], and that of eROSITA bubbles was acquired from the received photons with energy of 0.3–2.3 keV [3]. In this work, the RMHD model is applied to simulate the plasma distributions, as shown in Figure 3. Then, the spectral emissivities in $\gamma$-rays and X-rays, respectively, are computed by using the simulated plasma distributions and are compared with the observed images [2,3].

It is assumed that the fraction $a_1$ of non-thermal electron number density is the same as the fraction $a_2$ of non-thermal internal energy density, and the fraction $a_{1p}$ of high-energy proton number density is the same as the fraction $a_{2p}$ of high-energy proton internal energy density. In other words, $a_1=a_2$ and $a_{1p}=a_{2p}$. The plasma distribution in the RMHD model is composed of thermal and non-thermal particles. The fraction of non-thermal particles is quite small; hence, the plasma velocity is in the order of 100 km/h, much slower than the speed of light.

A set of parameters are adjusted to successfully reconstruct both images in $\gamma$-rays and X-rays, respectively, on the skymap. The fractions $a_1$ and $a_{1p}$ are estimated by comparing the reconstructed image and the observed one in $\gamma$-rays, as well as the relative contributions between leptonic and hadronic processes. As will be shown in Figure 10 and Figure 11, the hadronic process is much more efficient than the leptonic process in $\gamma$-ray emission, leading to $a_1\gg a_{1p}$. The simulation results also indicate that the hadronic process alone cannot reconstruct the observed image in X-rays. Hence, both leptons and hadrons are considered in the simulations, as will be shown in Figure 12.

Figure 10 shows the distributions of spectral emissivities (erg/cm³/sr/s/eV) at $ϵ=2$ GeV ($4.8\times {10}^{23}$ Hz) and the skymap of radiation intensity at 2–5 GeV, with $a_1=1\times {10}^{-7}$ and $a_{1p}={10}^{-18}$. The total internal energy of non-thermal electrons is about $2.1\times {10}^{48}$ (erg), of which about ${10}^{37}$ (erg) is contributed by hadrons. In [8], the energy of cosmic-ray electrons exerted in forming the observed Fermi bubbles was estimated to be about ${10}^{51}$ (erg).

Figure 10a shows that the spectral emissivity of inverse Compton scattering at $\nu =4.8\times {10}^{23}$ Hz extends to $r\simeq 7$ and $z\simeq 13$, respectively. The magnitude is about ${10}^{-36}$ erg/cm³/sr/s/eV. The emissivity in the outer shell is higher than that in the central lobe. Figure 10b shows the spectral emissivity of non-thermal bremsstrahlung. The highest emissivity appears behind the forward shock, near the Galactic Equator. Its magnitude is about ${10}^{-45}$ erg/cm³/sr/s/eV, much smaller than that of inverse Compton scattering.

Figure 10c shows the spectral emissivity of $\gamma$-ray produced via the decay of neutral pions. The highest emissivity appears behind the forward shock, near the Galactic Equator, and its magnitude is about ${10}^{-53}$ erg/cm³/sr/s/eV. Figure 10d shows the spectral emissivity of $\gamma$-rays produced by secondary electrons. The emissivity reaches a maximum of about ${10}^{-36}$ erg/cm³/sr/s/eV near the edges of the central lobe, and decreases to about ${10}^{-39}$ erg/cm³/sr/s/eV in the outer shell.

Figure 10e shows the skymap of radiation intensity at 2–5 GeV, which is drawn with an Aitoff projection [46]. The magnitude at low galactic latitudes ($\left|b\right|<{15}^{\circ }$) is about ${10}^4$ keV/cm²/s/sr, and decays towards high galactic latitudes. This skymap is contributed by cosmic rays only, excluding the isotropic background $\gamma$-rays in the Milky Way Galaxy [14].

In our simulations, non-thermal electrons are produced by the jet and the accompanying shock waves. The forward shock reaches $r\simeq 7$ at the present day, while the Earth is located at $r=8.32$. Since non-thermal electrons have not traversed beyond the forward shock, the observed radiation at $|\ell |>{80}^{\circ }$ is contributed by other mechanisms not considered in this work. The excess $\gamma$-rays observed by the Fermi LAT were attributed to pulsars or dark-matter annihilation in $\left|b\right|<{30}^{\circ }$ and $|\ell |<{10}^{\circ }$ [47], which are not considered in our simulations.

Figure 11 shows the spectral emissivities (erg/cm³/sr/s/eV) at $ϵ=2$ GeV ($4.8\times {10}^{23}$ Hz) in the inner part of the central lobe, which is enclosed by a yellow curve as shown in Figure 3. The total internal energy of non-thermal electrons is about $8.4\times {10}^{45}$ (erg).

Figure 11a shows the spectral emissivity of inverse Compton scattering, $j_{\mathrm{ic}}\left(ϵ\right)$, which is stronger near the bottom half area above the Galactic Equator. Figure 11b shows the spectral emissivity of non-thermal bremsstrahlung, $j_{\mathrm{ff}}^{\mathrm{nt}}\left(ϵ\right)$, with a magnitude of about ${10}^{-46}$ erg/cm³/sr/s/eV, much smaller than that of inverse Compton scattering.

Figure 11c shows the spectral emissivity of $\gamma$-rays, $j_{\pi \gamma }\left(ϵ\right)$, produced via the decay of neutral pions, with a magnitude of about ${10}^{-54}$ erg/cm³/sr/s/eV. Figure 11d shows the spectral emissivity of $\gamma$-rays, $j_{\sec }\left(ϵ\right)$, produced by secondary electrons. It is computed by using the emissivity of inverse Compton scattering in (26), with $N\left({\gamma }_e\right)$ replaced by the energy spectrum of secondary electrons, $Q_e\left(E_e\right)$, in (34). The magnitude is about ${10}^{-36}$ erg/cm³/sr/s/eV.

Summarized from Figure 11a–d, the emissivity at $ϵ=2$ GeV in the central lobe is mainly contributed by the inverse Compton scattering of non-thermal electrons and secondary electrons.

Figure 11e shows the skymap of radiation intensity at 2–5 GeV, contributed by non-thermal electrons in the central lobe. Two bubbles appear in the north and the south hemispheres, respectively, extending in the range $\pm {50}^{\circ }$ in latitude and ${30}^{\circ }$ in longitude. The magnitude is about ${10}^{1.5}$ keV/cm²/s/sr, which is much smaller than that from the outer shell, as shown in Figure 10e.

Figure 11f shows the residual map transformed from [2], which was obtained by removing the radiation from the galactic disk to highlight the contribution of the jet. The skymap of radiation intensity shown in Figure 11e roughly matches the residual map in Figure 11f in terms of size, shape and magnitude. It is reconfirmed that the Fermi bubbles in $\gamma$-ray emission are mainly contributed by inverse Compton scattering in the central lobe.

Figure 12 shows the spectral emissivities (erg/cm³/sr/s/eV) at $ϵ=0.6$ keV ($1.46\times {10}^{17}$ Hz) and the skymap of radiation intensity at 0.6–1 keV. Figure 12a shows that the spectral emissivity of inverse Compton scattering manifests a central lobe and an outer shell, which roughly delineate the shapes of Fermi bubbles and eROSITA bubbles, respectively. The magnitude of the emissivity is about ${10}^{-33}$ erg/cm³/sr/s/eV.

Figure 12b shows the spectral emissivity of non-thermal bremsstrahlung. The morphology is similar to that in Figure 12a, while the magnitude decreases to about ${10}^{-44}$ erg/cm³/sr/s/eV. Figure 12c shows the spectral emissivity of thermal bremsstrahlung, with a magnitude about ${10}^{-38}$ erg/cm³/sr/s/eV. Figure 12d shows the spectral emissivity attributed to secondary electrons, with a magnitude that is about ${10}^{-41}$ erg/cm³/sr/s/eV in the outer shell and decreases to about ${10}^{-43}$ erg/cm³/sr/s/eV beyond the bubbles. The results in Figure 12a–d indicate that the X-ray emission at 0.6 keV is dominated by the inverse Compton scattering.

Figure 12e shows the skymap of X-ray radiation intensity (keV/cm²/s/sr) at 0.6–1.0 keV. Two bubbles appear in the north and the south hemispheres, extending to about $\pm {75}^{\circ }$ in latitude and covering about ${100}^{\circ }$ in longitude. These bubbles match the observed eROSITA bubbles in size and shape.

In summary, a set of parameters are adjusted to successfully reconstruct bubbles in $\gamma$-rays and X-rays, respectively, on the skymap. The contributions of leptons and hadrons are comparable in reconstructing $\gamma$-ray bubbles. The relation $a_1\gg a_{1p}$ implies that the hadronic process is much more efficient than the leptonic process in $\gamma$-ray emission. The X-ray emission is mainly contributed by leptons. In the leptonic process, the inverse Compton scattering dominates both X-ray and $\gamma$-ray emissions. In the hadronic process, the secondary electrons contribute most of the $\gamma$-ray emission.

Figure 13 shows the longitudinal profiles of surface brightness I (counts/s/deg²) at 0.6–1.0 keV, at galactic latitudes of ${40}^{\circ }$, ${50}^{\circ }$ and ${60}^{\circ }$, respectively [3,8], where count means the number of photons. Figure 13a shows that, at $\left|b\right|={40}^{\circ }$ and $|\ell |<{50}^{\circ }$, the observed surface brightness is about 4–17 counts/s/deg² in the northern hemisphere, and about 6–10 counts/s/deg² in the southern hemisphere. The simulated surface brightness is about 8–9 counts/s/deg², which is comparable to the observed data. The maximum value occurs at $|\ell |\simeq {25}^{\circ }$, and decays quickly at $|\ell |>{25}^{\circ }$.

The observed surface brightness at $|\ell |>{50}^{\circ }$ decays to about 4–5 counts/s/deg². The simulated surface brightness at $|\ell |>{70}^{\circ }$ decays to $I<3$ counts/s/deg². In [3], the surface brightness estimated with an analytical model decays to $I<3$ counts/s/deg² at $|\ell |\simeq {60}^{\circ }$. The value observed at high $|\ell |$ is higher than simulation, possibly appended by detector noise, background cosmic X-ray and particles [3,48]. In [48], hydrodynamic simulations were performed to investigate the asymmetric eROSITA bubbles under the influence of background cosmic X-rays of about 2 counts/s/deg².

In short, the emission of eROSITA bubbles concentrates in $|\ell |<{50}^{\circ }$, and the observed surface brightness at $|\ell |>{50}^{\circ }$ is not emitted from the eROSITA bubbles.

Figure 13b shows that, at $\left|b\right|={50}^{\circ }$ and $|\ell |<{50}^{\circ }$, the observed surface brightness is about 4–15 counts/s/deg² in the northern hemisphere and about 5–8 counts/s/deg² in the southern hemisphere. The simulated surface brightness is about 11–12 counts/s/deg², comparable to the observed data in both hemispheres. The maximum value of $I\simeq 12$ counts/s/deg² appears at galactic longitude $|\ell |\simeq {15}^{\circ }$, and decays quickly at $|\ell |>{15}^{\circ }$.

Figure 13c shows that, at $\left|b\right|={60}^{\circ }$ and $|\ell |<{50}^{\circ }$, the observed surface brightness is about 4–13 counts/s/deg² in the northern hemisphere and about 5–7 counts/s/deg² in the southern hemisphere. The simulated surface brightness reaches a maximum of about 12 counts/s/deg² at $|\ell |=0^{\circ }$, and decreases below 3 counts/s/deg² at $|\ell |\simeq {40}^{\circ }$.

The observed X-ray skymap at 0.6–1 keV [3] was processed to enhance large-scale features of the bubbles; hence, it cannot be directly compared with the simulation data [3].

Figure 14 shows the spectral emissivity attributed to synchrotron radiation and the skymap of simulated synchrotron radiation at 23 GHz, transformed to the same unit as that of the skymaps in $\gamma$-rays and X-rays. In [49], spectral and morphological characteristics of diffuse galactic emission in the WMAP data were analyzed. Microwave haze or cosmic microwave background (CMB), with a magnitude about 0.2–6 (kJy/sr) at 23 GHz, was observed in the inner Galaxy on all the residual maps. The spectrum of microwave haze matched that of synchrotron emission from cosmic-ray electrons.

Figure 14a shows the differential spectral emissivity attributed to synchrotron radiation, $d{j}_{\mathrm{sy}}/dE$ (cm³/sr/s/eV), with $N_0=1$ and $|\overline{B}|=7$ ($\mathsf{\mu }$G). It is observed that the radiation at $\nu =23$ GHz is mainly contributed by the electrons with energy of $E\simeq {10}^9$–${10}^{10}$ (eV).

Figure 14b shows the skymap of simulated synchrotron radiation I at 23 GHz, contributed by non-thermal electrons. The maximum radiation intensity is about $3\times {10}^{-16}$ (keV/cm²/s/Hz/sr) $=5\times {10}^{-2}$ (Jy/sr), which is about four orders weaker than 0.2–6 (kJy/sr) in the observed microwave haze [49].

Various mechanisms have been proposed to explain microwave haze, including galactic wind, starburst, second-order acceleration, AGN jet and dark-matter annihilation [50]. In [51], a semi-analytical model was proposed to study microwave haze and its relevance to the Fermi bubbles. A reverse shock was predicted to re-accelerate cosmic-ray electrons and generate synchrotron radiation, explaining the spectrum and morphology of microwave haze.

The reverse shock in [51] was predicted at $r\simeq 1$ (kpc), which fell within the tangential discontinuity around $r\simeq 3$ (kpc). As shown in Figure 3a, the reverse shock was located in a very dilute region; hence, it could barely be observed in the simulations. In this work, possible mechanisms of galactic wind, starburst, second-order acceleration and dark-matter annihilation are not included; hence, the observed microwave haze cannot be well reconstructed.

Highlights and Prospects

An RMHD model, a leptonic emission model and a hadronic emission model have been integrated to reconstruct Fermi bubbles and eROSITA bubbles on the skymap, including shape, size and radiation intensity. The plasma distributions after a jet eruption from the galactic center are simulated with the RMHD model, manifesting forward shock and discontinuity transition that match the size and shape of the observed eROSITA bubbles and Fermi bubbles, respectively. A set of parameters are adjusted to successfully reconstruct skymaps of emission in $\gamma$-rays and X-rays, respectively, which match well with the observed Fermi bubbles and eROSITA bubbles, respectively, in shape, size and radiation intensity.

The effects of initial mass distribution in the Galaxy and the magnetic field on the evolution of plasma distributions have been investigated. The initial mass in the Galaxy tends to impede the expansion of bubbles, and the magnetic field induces more complicated features internal to the bubbles.

The effects of jet mass density (${\rho }_j$) and jet internal energy density (${ϵ}_j$) have been analyzed. If the jet mass density and the ejected kinetic power increase, while the ejected thermal power is fixed, the resulting mass density increases while the temperature decreases within the Fermi bubbles. If the ejected thermal power decreases while the jet mass density and the ejected kinetic power are fixed, the bubbles will expand slower. If the ejected kinetic energy is large, a depletion region with low temperature will be induced behind the forward shock.

Subtle differences between the spiral-arm model and the galactic disk model on the evolution of plasma distributions are compared. The shape, size and internal features in the resulting distributions of these two models are very similar. The mass density in the outer shell with the spiral-arm model is a little lower than that with the galactic disk model, and the temperature in the outer shell with the spiral-arm model is a little higher than that with the galactic disk model.

The energy spectrum of non-thermal electrons, the spectra of cosmological microwave background and interstellar radiation fields, as well as the spectrum of inverse Compton scattering haven been revisited.

Different emission mechanisms of leptonic and hadronic processes have been analyzed. The contributions to $\gamma$-rays from leptonic and hadronic processes are comparable, under the condition of $a_1\gg a_{1p}$, which implies the hadronic process is much more efficient than the leptonic process in emitting $\gamma$-rays. The X-ray emission is mainly contributed by the leptons. The emission of Fermi bubbles and eROSITA bubbles are attributed to both leptons and hadrons. In the leptonic process, inverse Compton scattering dominates in both X-ray and $\gamma$-ray emission. In the hadronic process, secondary electrons contribute most of the $\gamma$-ray emission.

5. Conclusions

An RMHD model, a leptonic emission model and a hadronic emission model have been integrated to reconstruct the observed Fermi bubbles and eROSITA bubbles. The plasma distributions induced by a jet eruption from the center of our Milky Way Galaxy are simulated, manifesting an outer shell and a central lobe which match the eROSITA bubbles and the Fermi bubbles, respectively, in shape and size. The skymap of radiation intensity computed with the simulated plasma distributions match well with the Fermi bubbles in $\gamma$-rays and the eROSITA bubbles in X-rays. The emission from Fermi bubbles and eROSITA bubbles is contributed by both leptons and hadrons.

The emission mechanisms of leptonic and hadronic processes have been revisited. Among the leptonic processes, inverse Compton scattering dominates in both X-ray and $\gamma$-ray emission. Among the hadronic processes, secondary electrons contribute most of the $\gamma$-ray emission. The number density of non-thermal leptons is much higher than that of non-thermal hadrons. The contributions from leptons and hadrons are comparable in $\gamma$-ray emission, implying that the hadronic process is much more efficient than the leptonic process in emitting $\gamma$-rays. The X-ray emission is mainly contributed by the leptons.

The initial mass distribution in the Galaxy tends to impede the expansion of bubbles, and the magnetic field tends to induce complicated features internal to the bubbles. Subtle differences attributed to different models or factors have also been investigated, including the spiral-arm model versus the galactic disk model on the plasma distributions, the energy spectrum of non-thermal electrons, the spectra of cosmological microwave background and interstellar radiation fields, and the spectrum of inverse Compton scattering. Synchrotron radiation from plasma distribution is computed, and possible mechanisms of microwave haze have been briefly discussed.