A Supernova-driven, Magnetically Collimated Outflow as the Origin of the Galactic Center Radio Bubbles

The simulation is based on a 3D MHD cartesian frame with a grid of 512³, equivalent to a physical volume of 200³ pc³ and a linear resolution of 0.39 pc. We set the z-axis to be perpendicular to the Galactic disk (north as positive), the y-axis to be parallel to the line of sight (the observer at the negative side), and the x-axis to run along decreasing Galactic longitude. Because the radio bubbles are roughly symmetric about the Galactic plane, we only simulate the z > 0 volume, sufficient to enclose the northern bubble, which exhibits a size of ∼120 pc (width) × 190 pc (height). We adopt an outflow boundary condition.

The simulation is governed by the ideal MHD conservation equations,

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}})=0,\\ \displaystyle \frac{\partial (\rho {\boldsymbol{v}})}{\partial t}+{\rm{\nabla }}\cdot {\left[\rho {\boldsymbol{v}}{\boldsymbol{v}}-\displaystyle \frac{{\boldsymbol{B}}{\boldsymbol{B}}}{4\pi }+1\left(p+\displaystyle \frac{{{\boldsymbol{B}}}^{2}}{8\pi }\right)\right]}^{T}=-\rho {\rm{\nabla }}{\rm{\Phi }},\\ \displaystyle \frac{\partial {E}_{t}}{\partial t}+{\rm{\nabla }}\cdot \left[\left(\displaystyle \frac{\rho {{\boldsymbol{v}}}^{2}}{2}+\rho e+p+\rho {\rm{\Phi }}\right)\displaystyle \frac{{\boldsymbol{v}}-{\boldsymbol{v}}\times {\boldsymbol{B}}\times {\boldsymbol{B}}}{4\pi }\right]\\ =\,-\displaystyle \frac{\partial \left(\rho {\rm{\Phi }}\right)}{\partial t},\\ \displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}-{\rm{\nabla }}\times ({\boldsymbol{v}}\times {\boldsymbol{B}})=0,\end{array}\end{eqnarray} \tag{ 1 }$$

where ρ is the mass density, p the thermal pressure, v the velocity, B the magnetic field, 1 the dyadic tensor, Φ the gravitational potential, and E_t the total energy density, defined as:

$$\begin{eqnarray}&&{E}_{t}=\rho \epsilon +\displaystyle \frac{{\left(\rho {\boldsymbol{v}}\right)}^{2}}{2\rho }+\displaystyle \frac{{{\boldsymbol{B}}}^{2}}{8\pi },\end{eqnarray} \tag{ 2 }$$

where is the internal energy. We use an ideal equation of state, i.e., = p/(Γ − 1), in which the ratio of specific heats Γ = 5/3.

As mentioned in Section 1, the orientation of the NTFs indicates that a vertical magnetic field is prevalent in the Galactic center. We adopt a dipole magnetic field structure generated by a current loop with a diameter of 300 pc, which can be expressed analytically (Simpson et al. 2001). With this large diameter, the magnetic field lines remain approximately vertical to the disk within our simulation volume. There is ample evidence that the Galactic center has an average magnetic field strength substantially higher than in the disk (Ferrière 2009). Crocker et al. (2010) derived a lower limit of 50 μG for the central 400 pc, based on an upper limit in the detected diffuse γ-ray flux. Given the observed radio spectral energy distribution of the Galactic center, a weaker magnetic field would lead to more relativistic electrons and consequently a higher γ-ray flux due to inverse Compton emission. In fact, energy equipartition between the magnetic field, X-ray-emitting hot plasma, and turbulent gas implies a magnetic field strength of ∼100 μG (Crocker et al. 2010). On the other hand, Thomas et al. (2020) suggested a stronger magnetic strength of 200 μG in the NTFs. In our fiducial run of simulation, the initial magnetic field strength at the origin (x = y = z = 0) is set as B₀ = 80μG. Values of 50 μG and 200 μG are also tested to examine the effect of a weaker/stronger magnetic field (see Section 2.4).

The simulation neglects viscosity and thermal conduction, but takes into account radiative cooling. We adopt the TABULATED cooling function implemented in PLUTO, which is generated with Cloudy for an optically thin plasma and solar abundances (Ferland et al. 2017). We neglect the synchrotron cooling of relativistic electrons, which are presumably produced by the SN shocks (see Section 2.4).

The gravitation in the Galactic center mainly originates from two components, namely, the nuclear star cluster (NSC), which dominates the innermost ∼20 pc, and the nuclear stellar disk (ND) that occupies the inner few hundred parsecs. We neglect larger-scale structures such as the bar and the Galactic disk. The SMBH, which has four million solar masses and a sphere of influence of a few parsecs in radius, can also be ignored given the scales of interest here. The NSC/ND will not evolve significantly on the timescale involved in our simulations, hence we adopt a fixed gravitational potential, which, following Stolte et al. (2008), can be approximated by a logarithmic form,

$$\begin{eqnarray}&&{\rm{\Phi }}=0.5{v}_{0}^{2}\mathrm{log}\left({R}_{c}^{2}+\displaystyle \frac{{x}^{2}}{{a}^{2}}+\displaystyle \frac{{y}^{2}}{{b}^{2}}+\displaystyle \frac{{z}^{2}}{{c}^{2}}\right),\end{eqnarray} \tag{ 3 }$$

where v₀ is the asymptotic velocity of a flat rotation curve, R_c is the core radius, and a, b, and c are stretching parameters. We adopt v₀ = 98.6 km s⁻¹, R_c = 2 pc, a = b = c = 1 for the NSC, and v₀ = 190 km s⁻¹, R_c = 90 pc, a = b = 1, c = 0.71 for the ND, from Table 1 of Stolte et al. (2008). The combined NSC+ND potential has been found to provide a good match to the observed stellar mass distribution in the Galactic center (Launhardt et al. 2002).

At the beginning of the simulation, the interstellar medium (ISM) is assumed to be isothermal and in hydrostatic equilibrium with the gravitational potential,

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\nabla }}P}{\rho }=-{\rm{\nabla }}{\rm{\Phi }},\end{eqnarray} \tag{ 4 }$$

where P = n_t kT is the thermal pressure, and n_t is the total number density of gas particles including protons, electrons, and heavy elements. As usual we define ρ = μ m_p n_t , where m_p is the proton mass and μ ≈ 0.6 is the mean molecular weight for solar abundance. The initial temperature is set to be 10⁶ K, which is roughly the virial temperature given the enclosed gravitational mass of 1 × 10⁹ M_⊙ within 100 pc. The prevalence of hot gas (with temperatures ≳10⁶ K) in the Galactic center has been established observationally (e.g., Baganoff et al. 2003; Ponti et al. 2015). While cooler gas (with temperatures ≲10⁴ K) is also known to exist in the Galactic center, it tends to concentrate in dense filaments and clouds near the midplane and is not expected to play a significant role in the bubble formation. We discuss possible effects of a multi-phase ISM on the observed properties of the bubble in Section 4.4. The initial density distribution can then be derived by solving Equation (4), as shown in Figure 1 along with the initial magnetic field distribution. From the adopted initial conditions, it can be shown that the thermal pressure of the ISM (n_t kT ∼ 10⁻¹² − 10⁻¹⁰ dyn cm⁻²) is everywhere significantly lower than the magnetic pressure (${B}_{0}^{2}/8\pi \sim 2.5\times {10}^{-10}\,\mathrm{dyn}\,{\mathrm{cm}}^{-2}$), perhaps except in the innermost few parsecs. In the meantime, the Alfvén speed, ${V}_{{\rm{A}}}={\left({B}_{0}^{2}/4\pi \rho \right)}^{\tfrac{1}{2}}\lesssim {10}^{3}\,\mathrm{km}\,{{\rm{s}}}^{-1}$, is much lower than the typical expansion velocity of the SN. Therefore, the present case of the Galactic center satisfies the moderately strong field condition defined by Insertis & Rees (1991).

Zoom In Zoom Out Reset image size

In the simulations, SNe are set to explode within a predefined cylindrical volume. The cylinder has a diameter of 50 pc in the x–y plane and a thickness of 10 pc along the z-axis, to mimic the concentration of massive stars near the Galactic plane (Kruijssen et al. 2015). We have tested a wider explosion area in the x–y plane (e.g., 100 pc in diameter, closer to that of the CMZ), finding that the resultant bubble would become significantly fatter, inconsistent with the observed morphology. In reality, the CMZ may provide a horizontal confinement to the bubble. However, a self-consistent implementation of the CMZ would necessarily introduce more free parameters, and is beyond the scope of the present work. The base of the radio bubbles shows a small but appreciable offset to the west of Sgr A* (Heywood et al. 2019). Thus we place the center of the cylinder at x = 5 pc to mimic this behavior. Due to the otherwise axisymmetry in the simulation, this appears to be the most viable way to reproduce the observed offset.

The fiducial SN birth rate is set to be 1 kyr⁻¹ (Di Teodoro et al. 2018), which is estimated by assuming an SFR of 0.1 M_⊙ yr⁻¹, a Kroupa (2001) initial mass function (IMF), and a minimum mass of 8 M_⊙ for the progenitor star of a core-collapse SN. Barnes et al. (2017) and Sormani et al. (2020) estimated a current SFR of 0.1 M_⊙ yr⁻¹ inside the CMZ, while Nogueras-Lara et al. (2020) found that star formation in the ND (which has a similar radial extent as the CMZ) has been relatively active in the past 30 Myr, with an SFR of 0.2–0.8 M_⊙ yr⁻¹. Our assumed SFR of 0.1 M_⊙ yr⁻¹ is compatible with the smaller radial extent of our adopted exploding region, which may be the case if SN events have been episodic and clustering on an approximately megayear or less timescale. We also test the effect of a lower SN birth rate of 0.5 kyr⁻¹ (see below). We have neglected Type Ia SNe, which have a birth rate of ≲0.05 kyr⁻¹ according to the enclosed stellar mass in the ND/NSC (Mannucci et al. 2005), though a recent study by Zhou et al. (2021) found evidence that Sgr A East, one of the few currently known SNRs in the Galactic center, was created by a Type Iax SN.

Individual SNe are thus injected at random positions inside the cylindrical volume, one after another, with a fixed interval according to the assumed birth rate. Each SN has an ejecta mass of M_ej = 10 M_⊙ and a kinetic energy of E_ej = 1 × 10⁵¹ erg (Poznanski 2013). This energy is deposited into a sphere with a radius of ${R}_{\mathrm{SN}}=4$ pc, ignoring any intrinsic anisotropy. The analytic solution within ${R}_{\mathrm{SN}}$ is derived from Truelove & McKee (1999), in which the newly born SN is divided into two parts, the inner uniform density core region and the outer power-law density envelope region. The radius of the former is 10 times that of the latter, and the power-law index is set as zero.

In this work, we perform four runs of simulation, each with a unique combination of magnetic field strength and SN explosion interval. Our fiducial simulation is represented by run B80I1, where B and I indicate the magnetic field and explosion interval, respectively. The fiducial run has B₀ = 80 μG and I = 1 kyr. The other three runs have either one of the two parameters varied. B50I1 has B₀ = 50 μG and B200I1 has B₀ = 200 μG, covering the empirical lower and upper limits inferred for the Galactic center (Section 2.1). Finally, B80I2 has an explosion interval of 2 kyr. The total elapsed time is set to be 330 kyr for all four runs. In the fiducial simulation, this is about the time when the top of the bubble approaches the edge of the simulation box. The time step is adaptive and ranges between 1 and 40 yr. The simulation parameters are summarized in Table 1.

To facilitate comparison with the observations, we generate synthetic radio and X-ray maps for the final snapshot (i.e., t = 330 kyr) of the simulation. We include synchrotron radiation and free–free emission in the radio band (default at 1284 MHz, to be consistent with the MeerKAT observation), while, for the X-ray band, only thermal emission from a collisionally ionized, optically thin plasma is considered.

First we need to distinguish regions inside and outside the evolving bubble. This is realized by adding a tracer parameter, Q, evaluated at every pixel in the simulation, which obeys a simple conservation equation:

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho Q)}{\partial t}+{\rm{\nabla }}\cdot (\rho Q{\boldsymbol{v}})=0.\end{eqnarray} \tag{ 5 }$$

Q has a value of 1 for pure SN ejecta and 0 for the unpolluted ISM and a value between 0 and 1 for pixels with mixed ejecta and ISM. We further calculate the Mach number for every pixel. The synthetic maps only take into account pixels with a non-zero tracer parameter or a Mach number greater than 2. The latter condition is employed to ensure that pixels with a high Mach number but a zero tracer parameter, such as those at or immediately behind the shock, are included.

Synchrotron emissivity depends on the magnetic field strength and the density of relativistic electrons. However, the latter cannot be directly obtained from our simulation and thus requires some working assumption. Here, we assume that the relativistic electron density at a given pixel of interest is proportional to the local gas density (Orlando et al. 2007; Zhang et al. 2017), normalized to have a mean energy density of 0.1 eV cm⁻³ across the bubble volume. This is compatible with the estimated mean cosmic-ray energy density of 10 eV cm⁻³ in the bubble (Heywood et al. 2019) and the empirical fact that relativistic electrons account for ∼1% of the total cosmic-ray energy density in the GeV band (Blasi 2013). We calculate the synchrotron emissivity in each pixel and integrate along the line of sight (i.e., the y-axis) to derive the synchrotron intensity map. In this calculation, the y-component of the magnetic field is neglected due to the nature of synchrotron radiation.

Radio free–free emission is calculated following the standard formula of Longair (2011), which, at a give pixel, scales with density squared and is a function of temperature. A temperature threshold of 10⁴ K is adopted when calculating the free–free emission. We find that only a tiny fraction of all pixels in any of our simulations has a temperature below 10⁵ K. The X-ray emissivity of an optically thin thermal plasma in collisional ionization equilibrium (Smith et al. 2001), also scaling with density squared, is extracted from ATOMDB, ⁵ version 3.0.9, for which we adopt a solar abundance. The free–free and X-ray intensity maps are again derived by integrating along the y-axis. We find that self-absorption is negligible in both the radio and X-ray bands, thanks to the relatively low column density involved.

In Figures 2 and 3, we show the gas density and temperature maps of run B80I1. In each figure, the density or temperature distribution is shown for a slice through the z = 0 (left columns), y = 0 (middle columns), and x = 0 (right columns) plane, after a simulation time of t = 30 (top rows), 180 (middle rows), and 330 (bottom rows) kyr.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By design, 30 SNe have exploded by the time of 30 kyr. The forward shock front of several of the youngest SNe are clearly revealed in the density map, as well as by the overlaid projected velocity vectors. A high-density region forms and persists around the origin (x = y = z = 0), because of the steep gravitational potential even in the presumed absence of an SMBH. As the shocks propagate, they compress and heat the ambient gas and also frequently collide with each other, eventually forming an expanding complex of post-shock gas with temperatures of 10^7–8 K.

By 180 kyr, this hot gas complex has developed into a bubble structure with a common dense shell, most clearly seen in the x–z and y–z planes. Inside the bubble, the density is low as a result of expansion, while the temperature remains high due to repeated shock heating. Numerous arc-like features are evident in the temperature map, especially in the x–z and y–z planes, which are the relic of individual SN shocks. At this stage, the bubble looks fat, with a similar extent (∼100 pc) along the three dimensions. However, the overall expansion starts to show a preference along the vertical (positive z) direction, with the vertical expansion velocity of the shell now being ∼690 km s⁻¹, substantially larger than the average expansion velocity of ∼120 km s⁻¹ in the x–y plane. This is primarily due to the collimation effect by an ordered magnetic field (Insertis & Rees 1991; Stone & Norman 1992; Rozyczka & Tenorio-Tagle 1995; Wu & Zhang 2019). Specifically, the SN shocks tend to push the semi-vertical magnetic field to the sides, greatly suppressing the magnetic field inside the bubble and, in the meantime, amplifying the magnetic field near the bubble shell. In turn, the latter decelerates and even halts the horizontal expansion of the bubble. The vertical expansion, on the other hand, feels no such magnetic confinement, thus a high velocity along this direction remains. The relatively strong gravitational potential in the x–y plane also contributes to retarding the horizontal expansion and facilitates the bubble collimation along the z-axis.

As a result, by the time of 330 kyr, the bubble becomes much more elongated. The top of the bubble almost reaches the edge of the simulation box (z = 200 pc), with a vertical expansion velocity still as high as ∼600 km s⁻¹, whereas its horizontal extent has not grown significantly since t = 180 yr. The width of the bubble at its base is about 120 pc, with a small but appreciable offset toward the positive x-axis, both in agreement with the observed bubble. Arc-like features tracing the sequential SN shocks remain prominent throughout the bubble interior. Near some of these arcs, locally enhanced magnetic fields are evident, which is the result of shock compression, as illustrated in Figure 4. The magnetic field strength takes a highest value of 175 μG across the bubble. Our simulation ends at this point.

Zoom In Zoom Out Reset image size

Similarly, we show the snapshots of density and temperature maps of simulation runs B80I2, B50I1, and B200I1 in Figures 5 and 6, all at t = 330 kyr. These three simulations share some common features with the fiducial simulation. In particular, a vertically collimated, bubble-like structure is formed in all these simulations. The bubble is delineated by a dense outer shell with compressed magnetic field and has a low-density, high-temperature interior with vertically oriented velocities and generally weak magnetic field. The bubble interior is not smooth, rather, it is filled with chaotic small-scale structures, again due to the sequential SN shocks and mutual interactions between them. Below, we shall describe the more unique features in the individual simulations.

Zoom In Zoom Out Reset image size

Simulation B80I2 (top row in Figures 5 and 6) adopts a lower explosion frequency than the fiducial case. This leads to a smaller energy injection rate, but is still sufficient to form a bubble. The bubble evolves more slowly, reaching a height of only 140 pc by the time of 330 kyr. The width of the bubble is also somewhat smaller than in the fiducial case (thus also narrower than the observed bubble), which remains the case even if we followed the bubble growth to a height of 200 pc. This occurs because, given a weaker SN energy injection but the same magnetic confinement, reduction in the horizontal expansion is greater than in the vertical expansion. We note that, at a further reduced explosion frequency, much of the SN ejecta would not be able to escape from the strong gravity near the midplane and a bubble would never form.

Zoom In Zoom Out Reset image size

In B50I1 (middle row in Figures 5 and 6), which has a weaker magnetic field compared to the fiducial run, the resultant vertical collimation is less effective and thus the bubble appears fatter. We note that a thinner bubble could still be achieved, should a lower explosion frequency be adopted in combination with the weaker magnetic field, for the reason explained above. However, in this case, it would take a much longer time for the bubble to grow to the observed height of 190 pc.

In contrast, B200I1 results in a significantly thinner structure. The magnetic field in this run is so strong that it can resist the compression of the SN shocks and consequently there is little sweeping of magnetic field inside the bubble. With the strong magnetic collimation, some SN ejecta are able to rapidly propagate along the field lines, forming vertical protrusions (several of these are captured in the x − z and y − z slides). The overall morphology is obviously inconsistent with the observed bubble.

Here we shall provide a more quantitative comparison with the radio and X-ray observations, with a focus on the fiducial simulation, which has the best morphological agreement with the observed bubble.

The synthetic synchrotron and free–free intensity maps of B80I1, after an evolution time t = 330 kyr, are shown in the left and right panels of Figure 7. The overall morphology is quite similar between the synchrotron and free–free emission, which is partially owing to our assumption that the density of relativistic electrons scales with the local gas density. However, the synchrotron intensity is everywhere orders of magnitude higher than the free–free counterpart in the synthetic maps. This holds true even considering the uncertainties in the energy density of the relativistic electrons and the magnetic strength. Consequently, synchrotron dominates the total flux density at 1284 MHz, consistent with the MeerKAT observation (Heywood et al. 2019). It is noteworthy that both the hydrogen recombination line, H90α, at 8309 MHz and the 8.4 GHz continuum are found to trace the GCL (Nagoshi et al. 2019), which exhibits a loop-like structure spatially coincident with the northern radio bubble. This suggests that the thermal component may have an increasingly larger contribution toward higher frequencies, which can be due to a combined effect of substantial synchrotron cooling at higher frequencies and the presence of ambient cooler gas not taken into account in our simulation.

Zoom In Zoom Out Reset image size

The overall extent of the synthetic synchrotron emission highly resembles that of the northern radio bubble (delineated by the red dotted line in Figure 7), which, has a width of 120 pc at its base and a height of 190 pc. Another interesting feature in the simulation is the presence of numerous filaments both at the edge of and inside the bubble, which closely resemble the NTFs (Yusef-Zadeh et al. 1984), although the ones in the simulation appear thicker and fuzzier in general, which may be partly owing to our moderate resolution. In the simulation, these filaments originate from the sequential SN shocks and their mutual interactions, and are associated with locally amplified magnetic field (Figure 4). Their possible relation with the NTFs will be further addressed in Section 4.

The 1284 MHz synchrotron flux density of the simulated bubble is found to be 5801 Jy, which is to be contrasted with our rough estimate of the observed flux density in the MeerKAT image, 970 Jy, obtained by assuming a mean flux density of 3 mJy beam⁻¹ across the projected area of the bubble. We caution that the MeerKAT mosaic image presented in Heywood et al. (2019) was not corrected for the primary beam attenuation and that the extended emission from the bubble suffers from potential flux loss in the interferometric image (I. Heywood, private communication), thus our estimate should be treated as a lower limit of the true flux density. On the other hand, the simulated flux density depends heavily on the assumed energy density of relativistic electrons. Therefore, the apparently large discrepancy between the observed and simulated radio flux densities should be taken as a point for future improvement rather than a failure of the simulation.

The synthetic 0.5–1.5 keV and 1.5–10 keV X-ray intensity maps are shown in Figure 8. Compared to its radio morphology, the simulated bubble appears smoother in the X-rays. The expanding shell of the bubble (Figure 2) leaves no significant sign of limb brightening in the 1.5–10 keV map, which is roughly consistent with the X-ray observations. This might be due to the fact that the shell is, on average, cooler than the bubble interior (Figure 3). Indeed, in the 0.5–1.5 keV map, which is more sensitive to gas temperatures below ∼1 keV, limb brightening is more evident especially at the northwestern side of the bubble, although this energy band is not directly observable due to the large foreground absorption column density (a few 10²² cm⁻²; Ponti et al. 2019). The 1.5–10 keV map also exhibits much fewer small-scale structures in the bubble interior, except near the x–y plane where the gas density is high and the most recent SNe freshly deposit a fraction of their kinetic energy. In particular, remnants of two newly exploded SNe are evident near the center (marked in the right panel of Figure 8), although they are not clearly seen in the synthetic radio map. An SNR evolving near Sgr A* will be heavily shaped by the strong gravity, with a large part of the ejecta pulled to the midplane, resulting in an appearance resembling the bipolar X-ray lobes detected in the innermost 15 parsecs of the Galactic center (Ponti et al. 2015, 2019).

Zoom In Zoom Out Reset image size

The thermal, kinetic, and magnetic energy of the bubble is calculated by summing over all "bubble pixels" (Section 2.4), which is found to be 1.9, 1.2, and 0.1 × 10⁵² erg, respectively. The initial thermal and magnetic energy within the bubble volume are 0.7 and 1.1 × 10⁵² erg. A net decrease of the magnetic energy underscores the sweep-up of the magnetic field. Ponti et al. (2019) estimated a thermal energy of 4 × 10⁵² erg for the X-ray chimneys (sum of the northern and southern halves), which is well matched by the simulated value of 1.9 × 10⁵² erg for the northern chimney.

Ponti et al. (2019) also measured density and temperature profiles along selected Galactic longitude, l = 0°, and Galactic latitude, b = 07. For a direct comparison, we construct density and temperature profiles at l = 0° and b = 07 from the simulation, as shown in Figure 9. Precisely speaking, Sgr A* is located at l = 005579, b = −004608, but here we neglect this small difference and simply take the x = 0 plane and z = 98 pc plane for comparison. We calculate the density-weighted mean density along the line of sight as

$$\begin{eqnarray}&&\langle n\rangle =\sqrt{\displaystyle \frac{\int {n}_{{\rm{t}}}^{2}{dV}}{\int {dV}}},\end{eqnarray} \tag{ 6 }$$

where dV = Adl, A is the projected area, and the line-of-sight integration (dl) is from the farthest side to the nearest side of the bubble. The projected area varies across the profiles to approximate the rather irregular spectral extraction regions used in Ponti et al. (2019). At l = 0°, the width is 12.5 pc, and the lengths are 20 pc and 70 pc, respectively, for b < 026 and b > 026. At b = 07, the width is 12.5 pc, and the length is always 70 pc. The emissivity-weighted mean temperature is calculated as

$$\begin{eqnarray}&&\langle T\rangle =\displaystyle \frac{\int {{Tn}}_{{\rm{t}}}^{2}{\rm{\Lambda }}(T,Z){dV}}{\int {n}_{{\rm{t}}}^{2}{\rm{\Lambda }}(T,Z){dV}},\end{eqnarray} \tag{ 7 }$$

where Λ is the tabulated X-ray emissivity as a function of temperature and metallicity extracted from ATOMDB. By examining the distribution of the SN ejecta through the tracer parameter, we have verified that the assumption of a uniform metallicity is a reasonable approximation.

Zoom In Zoom Out Reset image size

At l = 0°, the simulated density profile peaks at the midplane and decreases untill z ≈ 40 pc, beyond which it flattens. This general trend is in reasonable agreement with the observed density profile. Notably, the observed density profile has a significantly higher peak at low z. This may be due partly to the smaller line-of-sight depth adopted by Ponti et al. (2019) for the two inner data points, and partly to contamination from unresolved stellar objects and nonthermal extended features to the apparently diffuse X-ray emission near the midplane (Zhu et al. 2018). The simulated temperature profile appears bumpy around a mean value of ∼1.0 keV. The "bumps" are most likely due to consecutive SN shocks propagating upward. Near the top of the expanding shell, the temperature quickly drops to ∼0.6 keV. The observed temperature profile, on the other hand, appears flatter and has a lower value of 0.7–0.8 keV between 20–150 pc. We note that the observed temperature was derived using a single-temperature spectral model to the underlying plasma having a range of temperatures (Ponti et al. 2019). The Galactic center hot ISM is expected to have a somewhat lower temperature than in the whole bubble interior. Inclusion of the hot ISM in the observed spectrum could have led to a lower observed temperature.

At b = 0.7°, the simulated density profile peaks at the eastern and western edges of the bubble shell, which is consistent with Figure 1. However, there is no clear sign of limb brightening in the observed density profile; an enhanced density is only weakly seen near the eastern edge (X ≈ −60 pc) but is absent near the western edge (x ≈ 70 pc). One possibility is that the soft X-ray emission from the denser and cooler western shell has largely dropped out of the observation band (but could have been seen in the 0.5–1.5 keV band, as shown in the left panel of Figure 8). The simulated temperature profile shows a roughly inverse "U" shape, with values peaking at ∼1.25 keV at x = 20 pc. It is noteworthy that the outermost few points in the simulated profile are actually outside the bubble volume, whose values only reflect the unperturbed ISM. The observed temperature profile, again derived from a spectral fit using a single-temperature model, appears flat around a mean value of 0.8 keV.

Ponti et al. (2019) did not provide an explicit total X-ray luminosity of the chimneys. A rough estimate of this value can be made by adopting a cylinder of 150 pc in both diameter and height, as assumed by Ponti et al. (2019), a mean density of 0.1 cm⁻³ and a mean temperature of 1 × 10⁷ K (0.86 keV), which are representative of the X-ray chimneys. This leads to an estimated 1.5–10 keV luminosity of ∼2.8 × 10³⁶ erg s⁻¹ for the northern chimney, again well matched by the simulated value of 2.0 × 10³⁶ erg s⁻¹.

For completeness, the synthetic radio and X-ray maps of runs B80I2, B50I1, and B200I1 are shown in Figure 10. While these maps exhibit some interesting features, it is immediately clear that none of them matches the observed bubble morphology (again approximated by the red dotted line).

Zoom In Zoom Out Reset image size

The simulations presented in the previous section show that an outflow driven by sequential SN explosions and collimated by a vertical magnetic field can provide a reasonable explanation for the observed radio bubbles/X-ray chimneys in the Galactic center. In particular, the simulations can well reproduce the overall morphology, X-ray luminosity, and thermal energy of the northern bubble.

This scenario relies on two key ingredients: SN explosions clustering in the nuclear disk to provide a semicontinuous energy input, and a vertical, moderately strong, magnetic field to provide the collimation. Both ingredients are very likely available in the Galactic center. Indeed, direct evidence for contemporary SN explosions in the Galactic center was provided by at least a few SNRs clearly visible in radio or X-ray images (e.g., Ponti et al. 2015). Moreover, about 200 emission-line objects have been detected in the Galactic center, most of which are likely evolved massive stars (Dong et al. 2012). These stars may belong to the same population that gave rise to the SNe responsible for launching the bubbles. As for the magnetic field, it is widely thought that it is predominantly poloidal in the Galactic center, at least in regions outside the giant molecular clouds (Ferrière 2009). In this regard, an SNe-driven, magnetically collimated outflow should naturally develop in the Galactic center, provided the correctness of our simulations.

As mentioned in Section 1, a competing driver of a large-scale outflow is the kinetic power from the central SMBH, even though Sgr A* is by no means comparable with a classical AGN. While our simulations cannot automatically rule out an AGN-driven outflow, they share useful insight on the latter case. Compared to the distributed SN explosions, energy input from the SMBH is highly concentrated. Thus an AGN-driven outflow on the hundred-parsec scale may either acquire a highly elongated shape in the case of a canonical jet-driven outflow (e.g., Zhang & Guo 2020), or inflate a fat bubble in the case of a more isotropic wind symbiotic with the hot accretion flow onto a weakly accreting SMBH (Yuan et al. 2015). Magnetic collimation may also shape the wind-blown bubble, but one expects that the resultant structure is again a highly elongated one. Thus matching the morphology of the radio bubbles with an AGN wind-driven outflow may require some fine-tuning, which awaits a detailed investigation.

We now turn to consider the fate of the radio bubbles. In the framework of our simulations, the SNe-driven outflow is necessarily an evolving structure. In fact, at the end of our fiducial simulation, the top of the bubble still expands at a speed of ∼600 km s⁻¹ (Section 3.1). Provided a continuous energy injection from future SNe, which is quite likely given the evolved massive stars near the disk plane (Dong et al. 2012), the bubbles should continue to grow and gradually evolve into a more "chimney"-like structure, as long as a moderately strong magnetic field persists to greater heights. Conversely, if SNe were temporarily shut off, one expects that the bubble/chimney would ultimately disperse and collapse within a time not much greater than the sound-crossing time (a few hundred kiloyears). We have run a test simulation to examine such a case. Specifically, we adopt the same setting as the fiducial simulation, except that SN explosions cease after a time of 200 kyr. It is found that the upper edge of the bubble can still climb to a height of ∼190 pc with its accumulated momentum. However, the interior of the bubble, especially its lower portion, begins to collapse soon after the shutoff of the SNe, due to the loss of energy injection against the strong central gravity. In addition, the mean gas temperature inside the bubble gradually declines. Such an effect might bring the simulated temperature profile into better agreement with the observed temperature profile (Figure 9), although we have no evidence that the Galactic center is currently experiencing a substantial drop in the SN birth rate.

It is interesting to ask whether the radio bubbles/X-ray chimneys have a causal relation with the Fermi bubbles (Su et al. 2010) and eROSITA bubbles (Predehl et al. 2020) found on much larger scales. We note that the age of the radio bubbles inferred from our simulations is only a few hundred kiloyears, much shorter than the dynamical timescale of a few megayears originally suggested by Heywood et al. (2019). However, Heywood et al. (2019)'s estimate was based on the assumption of a constant expansion velocity of the bubbles, which is implausible, hence a shorter timescale is expected. The estimated age of the Fermi bubbles, on the other hand, ranges from 1 Myr (Yang et al. 2013) to 1 Gyr (Crocker & Aharonian 2011). Thus, in the context of our supernova-based model for the origin of the radio bubbles/chimneys, the radio bubbles would be a dynamically younger and independent structure simply evolving in the interior of the Fermi/eROSITA bubbles, which themselves were formed by older activities in the Galactic center.

Alternatively, as suggested by Ponti et al. (2019), the X-ray chimney may be a channel that transports energy from the Galactic center to the high-latitude region currently occupied by the Fermi bubbles. In this case, the channel should have existed for tens of megayears, so that star formation in the Galactic center can be sufficient to supply the total energy content of the Fermi bubbles, ∼10⁵⁶ erg (Carretti et al. 2013). However, such a picture contradicts with the capped morphology of the radio bubbles (the southern bubble is not obviously capped in X-rays; Ponti et al. 2021), which, according to our simulations, is naturally explained as the expanding shell of a newly born outflow. This picture may be reconciled if star formation in the Galactic center has been episodic on a timescale of ∼10 Myr (Krumholz & Kruijssen 2015). In this case, the "chimney" is (re)established by consecutive generations of mini-starbursts and collapses inbetween. Of course, over such a long interval, the activity of Sgr A* can also play an important role in contributing to the inflation of the chimneys, especially in view of the fact it was likely much more active in the recent past (Ponti et al. 2010, 2013; Camilo et al. 2018). In a hybrid scenario, Sgr A*, with SNe and even stellar winds, can simultaneously sustain the "chimney" and transport energy to larger scales, implying X-ray emission beyond the edge of the radio bubbles, which is also suggested by Ponti et al. (2021).

The origin of the NTFs has been extensively debated since their discovery nearly four decades ago. Proposed models for the NTFs include expanding magnetic loops (Heyvaerts et al. 1988), induced electric fields (Benford 1988; Morris & Yusef-Zadeh 1989), thermal instability in relativistic gas (Rosso & Pelletier 1993), cosmic strings (Chudnovsky et al. 1986), magnetic reconnection (Lesch & Reich 1992; Serabyn & Morris 1994; Morris 1996; Banda-Barragán et al. 2016, 2018), analogs of cometary plasma tails (Shore & LaRosa 1999), a turbulent magnetic field (Boldyrev & Yusef-Zadeh 2006), stellar winds or SNe of the young star cluster (Yusef-Zadeh 2003; Yusef-Zadeh & Wardle 2019), pulsar wind nebulae (Barkov & Lyutikov 2019), and the tidal destruction of gas clouds (Coughlin et al. 2021). Of course, a multi-SNe hypothesis has also been suggested (Sofue 2020).

In our simulations, filamentary features resembling the observed NTFs trigger and form primarily at the interface of colliding shocks of individual SNe (Figure 4). Magnetic fields are compressed and amplified in these filaments, where particle acceleration (e.g., due to diffusive shock acceleration) is expected to take place. Also the Radio Arc finds its possible counterpart in the simulations, which arises from the piling of consecutive SN shocks at the sides of the bubble (Figure 2). Comparing Figures 7 and 10, it occurs that an SN-driven outflow evolving in a weaker magnetic field produces more filaments. This is because a strong magnetic field can more easily confine an SN shock and reduce its chance of encountering other shocks. We note that, in the simulation, many filaments are indeed one-dimensional structures, i.e., they have a distinct long axis roughly oriented vertically, but some others arise from a projection effect, i.e., a two-dimensional surface viewed edge on. Such a surface is also the result of colliding shock fronts. We stress that the moderate resolution of our simulation would smear the appearance of the shock fronts, so we anticipate that additional apparent filaments would show up with higher resolution. The viability of this formation mechanism for the NTFs could be assessed by direct comparison of the cross-sectional profiles of the filaments appearing in the simulations with those of observed NTFs, but a higher resolution simulation is needed for such a comparison.

We note that there are NTFs found outside the radio bubbles (Heywood et al. 2019). These might have been formed in a past generation of clustering SN explosions, and they exist for a longer time than the associated outflow. Of course, we cannot rule out the aforementioned alternative models for all NTFs. In reality, the NTFs can have a mixed origin, i.e., different processes, including SN shocks, stellar winds, and pulsar winds can produce seeds of NTFs which are further shaped by the compressed magnetic field or other mechanisms.

The magnetic field is a crucial component of the Galactic center environment. At present, the average field strength is still quite uncertain. The assumption of energy equipartition between the magnetic field and relativistic particles leads to estimates up to ∼1 mG in the brightest NTFs and as low as 10 μG in the more diffuse background. Crocker et al. (2010) derived a lower limit of ∼50 μG based on the diffuse γ-ray flux and suggested a typical value of ∼100 μG in the central 400 pc region.

In our simulation B50I1, which adopts a field strength of 50 μG, an outflow can be developed, although the resultant bubble appears fatter due to the reduced magnetic confinement compared to the fiducial simulation (Section 3.2). This lends some support to the above lower limit.

On the other hand, simulation B200I1, which assumes a field strength of 200 μG, is obviously inconsistent with the observation (Figure 10). This conclusion holds even if the other parameter, the SN birth rate, were adjusted within a reasonable range. Qualitatively, at a lower SN birth rate, the shock and ejecta of individual SNe would be less resistant to the magnetic pressure, thus they are less likely to evolve into a mutual network. The resultant outflow hardly takes a bubble shape, rather it would consist of many barrel-like structures, through which individual SN ejecta propagate. Only a much higher SN birth rate can counteract the magnetic pressure, but this would be inconsistent with the currently accepted star formation rate in the Galactic center (∼0.1 M_⊙ yr⁻¹). Therefore, our simulations provide a meaningful constraint on the average magnetic field on 100 pc scales in the Galactic center, 50 μG ≲ B₀ ≲ 200 μ G.

Our fiducial run B80I1 demonstrates localized magnetic field amplification across the bubble, reaching a maximum field strength of 175 μG. It is expected that the global magnetic field would gradually restore to the initial configuration after the termination of clustering SN explosion and the dispersion/collapse of the outflow.

Despite the satisfactory reproduction of the major observed properties of the radio bubbles/X-ray chimneys, some notable discrepancies exist between our simulation results and the observations, which warrant the following remarks.

The observed edge-brightened radio bubbles have a low-surface-brightness interior, while in our simulation, the edge–interior contrast is less significant. A possible cause is that we have ignored synchrotron cooling. Using a magnetic field of 20 μG, Heywood et al. (2019) derived a synchrotron cooling time of 1–2 Myr by assuming that the electron energy density distribution has a power-law index of 2. Based on the same method, we estimate a cooling time of 250 kyr for 80 μG, which is comparable to the evolution time of the bubble in our simulation. Hence the relativistic electrons produced at the early stage and now filling the bubble interior should be subject to radiative cooling, an effect that is not taken into account but otherwise would enhance the edge–interior contrast.

An alternative and more likely cause is the absence of a cool gas shell in our simulation. The presence of cool gas (with a temperature of ∼10⁴ K) in the outer part of the GCL has been known for some time (Law 2010; Nakashima et al. 2019). This cool gas is not found in our simulations, owing to the very moderate radiative cooling even in the dense shell of post-shock gas. This is also the reason why the free–free emission predicted by our simulation is negligible compared to the synchrotron (Section 3.3). Hence the detected cool gas probably has an external origin that is missing in the framework of our simulation. Indeed a substantial amount of both cool and cold gas exist in the NSD/CMZ (Ferrière et al. 2007), and part of this gas may be swept into the bubble shell and/or entrained into the bubble interior. For example, Ponti et al. (2021) argued that a gas cloud associated with the bright 25 μm source AFGL5376 has been accelerated and is now defining part of the wall of the bubble. An additional source of cool gas is the stellar wind of the massive stars distributed in the nuclear disk.

In principle, the Galactic center outflow may also be driven by stellar winds (Chevalier 1992). Stellar winds as an additional energy and momentum source have not been included in our simulation. We can give a rough estimate of the collective energy input from the massive stars in the Galactic center. The stellar winds should be dominated by the Wolf-Rayet stars, which have a typical mass loss rate of 10⁻⁵ M_⊙ yr⁻¹ and a wind velocity of 2000 km s⁻¹. Thus the ∼200 evolved massive stars found by Dong et al. (2012) in the nuclear disk have a total kinetic power of 2.5 × 10³⁹ erg s⁻¹ and would release a kinetic energy of 2.6 × 10⁵² erg in 330 kyr. The massive stars in the central parsec provide an additional kinetic energy of 3 × 10⁵¹ erg in 330 kyr, assuming a collective mass loss rate of 10⁻³ M_⊙ yr⁻¹ and a wind velocity of 1000 km s⁻¹ (Najarro et al. 1997; Quataert 2004). Therefore, the energy input from the massive stars is about one order of magnitude smaller than that of the SNe in our simulation. Nevertheless, massive stars may start launching strong winds a few Myr before their core collapse, significantly shaping the ambient gas into which the bubbles expand. A self-consistent implementation of the stellar winds requires a reliable stellar evolution model and a much higher resolution, thus awaits future work.