A Supernova-driven, Magnetically-collimated Outflow as the Origin of the Galactic Center Radio Bubbles
aa r X i v : . [ a s t r o - ph . GA ] J a n Draft version January 28, 2021
Typeset using L A TEX twocolumn style in AASTeX61
A SUPERNOVA-DRIVEN, MAGNETICALLY-COLLIMATED OUTFLOW AS THE ORIGIN OF THEGALACTIC CENTER RADIO BUBBLES
Mengfei Zhang,
Zhiyuan Li, and Mark R. Morris School of Astronomy and Space Science, Nanjing University, Nanjing 210023, China Key Laboratory of Modern Astronomy and Astrophysics (Nanjing University), Ministry of Education, Nanjing 210023, China Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA
ABSTRACTA pair of non-thermal radio bubbles recently discovered in the inner few hundred parsecs of the Galactic center bearsa close spatial association with elongated, thermal X-ray features called the X-ray chimneys. While the morphology,position, and orientation vividly point to an outflow from the Galactic center, the physical processes responsible forthe outflow remain to be understood. We use three-dimensional magnetohydrodynamic (MHD) simulations to test thehypothesis that the radio bubbles/X-ray chimneys are the manifestation of an energetic outflow driven by multiplecore-collapsed supernovae in the nuclear stellar disk, where numerous massive stars are known to be present. Oursimulations are run with different combinations of two main parameters, the supernova birth rate and the strengthof a global magnetic field orientated vertically with respect to the disk. The simulation results show that a hot gasoutflow can naturally form and acquire a vertically elongated shape due to collimation by the magnetic pressure. Inparticular, the simulation with an initial magnetic field strength of 80 µ G and a supernova rate of 1 kyr − can wellreproduce the observed morphology, internal energy and X-ray luminosity of the bubbles after an evolutionary timeof 330 kyr. On the other hand, a magnetic field strength of 200 µ G gives rise to an overly elongated outflow that isinconsistent with the observed bubbles. The simulations also reveal that, inside the bubbles, mutual collisions betweenthe shock waves of individual supernova remnants produce dense filaments of locally amplified magnetic field. Suchfilaments may account for a fraction of the synchrotron-emitting radio filaments known to exist in the Galactic center.
Keywords:
Galactic center (565), Superbubbles (1656), Magnetic fields (994), Magnetohydrodynamicalsimulations (1966) [email protected]@nju.edu.cn
M.F.Zhang et al. INTRODUCTIONGalactic outflows driven by energy and momentumof an active galactic nucleus (AGN) and/or super-novae (SNe) are now understood to be an indispens-able component of the galactic ecosystem (Fabian 2012;Heckman & Best 2014; Heckman & Thompson 2017;Zhang 2018). Multi-wavelength observations over thepast decades have established an ever-growing inven-tory of galactic outflows, leading to the recognition thatthese outflows typically involve multi-scales and multi-phases. However, our physical understanding of galacticoutflows, in particular their mass budget, energetics andlife cycle, is still far from complete.The Galactic center, loosely defined here as the in-nermost few hundred parsec region of our Galaxy, pro-vides the closest and perhaps the best laboratory forstudying the formation and early evolution of a galacticoutflow. Observational evidence has accumulated overrecent years for a multi-phase outflow from the Galac-tic center (Bland-Hawthorn & Cohen 2003; Law 2010;Nakashima et al. 2019), collectively known as the Galac-tic Center Lobe (GCL; Sofue & Handa 1984), a loop-like feature extending vertically out to & ◦ corresponds to 140pc) north of the disk mid-plane. Compelling evidencealso exists for outflows at still larger (kiloparsec-) scales(Su et al. 2010; Carretti et al. 2013; Di Teodoro et al.2018; di Teodoro et al. 2020; Predehl et al. 2020), butthe physical relation between the outflows on differentscales, e.g., whether they were produced by the samemechanism, remains an open question.More recently, our view of the Galactic center outfloware further sharpened. Based on high-resolution radiocontinuum observations afforded by the MeerKAT ra-dio telescope, Heywood et al. (2019) found evidence fora pair of radio bubbles in the Galactic center, whichare roughly symmetric about the disk mid-plane witha width of 140 pc and a full length of 430 pc. Thenorthern bubble is spatially coincident with the GCL,but it is more clearly limb-brightened. In particular,the eastern side of the radio bubbles is delineated bythe famous Radio Arc (Yusef-Zadeh et al. 1984) andits northern and southern extension toward higher lat-itudes; the western side is also bounded by prominentnon-thermal filaments (NTFs; Yusef-Zadeh et al. 1984).Non-thermal emission is predominant in the radio bub-bles at the observed frequency of 1284 MHz, althoughthe GCL is known to show substantial thermal emis-sion at different wavebands (Bland-Hawthorn & Cohen2003; Law 2010; Nagoshi et al. 2019). Strikingly, theshells of the radio bubbles delineate the so-called “X-ray chimneys” recently discovered by X-ray observations (Ponti et al. 2019), which is a pair of diffuse, thermal X-ray features extending above and below the mid-plane.This strongly suggests a physical relation between thetwo features, reminiscent of a collimated hot gas outflowwith an expanding shell (Ponti et al. 2021).Proposed origins for the Galactic center outflow aswell as for the outflows on larger scales (i.e., the Fermibubbles and the recently discovered eROSITA bubbles;Su et al. 2010; Predehl et al. 2020) fall in two categories(Cheng et al. 2011; Zubovas et al. 2011; Fujita et al.2014; Heywood et al. 2019; Zhang & Guo 2020): (i)past activity from the central super-massive black hole(SMBH), commonly known as Sgr A*, which is currentlyin a quiescent state; or (ii) episodic or continuous nuclearstar formation (Genzel et al. 2010). In principle, bothprocesses can drive an energetic outflow and produce thebubble-like structures observed at multi-wavelengthsand multi-scales. Therefore a quantitative modelingand close comparison with the observations are crucialto distinguish between the two scenarios. In the litera-ture, there have been a number of numerical simulationsof a large-scale outflow from the Galactic center, whichfocuses on the formation of the Fermi bubbles by AGNjets or AGN winds (Guo & Mathews 2012; Sarkar et al.2015; Cheng et al. 2015; Zhang & Guo 2020).In this work, we investigate the specific scenario thatthe radio bubbles/X-ray chimneys are the manifestationof an outflow driven by sequential SN explosions concen-trated in the Galactic center, using three-dimensionalmagnetohydrodynamic (MHD) simulations, which is thefirst attempt of this kind to our knowledge. Recently,Li et al. (2017) and Li & Bryan (2020) have performedadvanced numerical simulations to study SNe-drivenoutflows on a similar physical scale, but these simula-tions were run with physical conditions typical of galac-tic disks. The Galactic center, on the other hand, is aunique environment characterized by a strong gravity, aconcentration of massive stars, and a strong and orderedmagnetic field. In particular, the presence of the NTFs,which have a strong tendency vertically oriented withrespect to the disk, points to a vertical magnetic field inthe Galactic center (see review by Ferri`ere 2009). The-oretical studies have demonstrated that a strong exter-nal magnetic field can significantly affect the evolutionof a supernova remnant (SNR; Insertis & Rees 1991;Rozyczka & Tenorio-Tagle 1995; Wu & Zhang 2019), asthe magnetic pressure confines the expansion of the SNejecta in such a way that they preferentially propagatealong the direction of the magnetic field.We are thus motivated to perform numerical simula-tions to test the scenario of an SN-driven, magnetically-collimated outflow for the radio bubbles/X-ray chim- alactic Center Radio bubbles SIMULATIONWe use the publicly available MHD code
PLUTO (Mignone et al. 2007, 2012) to simulate sequential SNeexplosions in the Galactic center and the formation ofan SN-driven bubble. The global dynamical evolutionand fine structures of the bubble necessarily depend onmany physical processes and physical quantities of theGalactic center, some of which are not well constrained.Rather than pursuing a full degree of realism or a thor-ough exploration of the parameter space, our main aimhere is to test a simplified but well-motivated model forthe bubble formation.2.1. Basic MHD Equations and Magnetic FieldConfiguration
The simulation is based on a three-dimensional (3D)MHD cartesian frame with a grid of 512 , equivalent toa physical volume of 200 pc and a linear resolution of0.39 pc. We set the z -axis to be perpendicular to theGalactic disk (north as positive), the y -axis to be paral-lel to the line-of-sight (the observer at the negative side),and the x -axis to run along decreasing Galactic longi-tude. Because the radio bubbles are roughly symmetricabout the Galactic plane, we only simulate the z > ∼
120 pc (width) ×
190 pc (height).We adopt an outflow boundary condition. http://plutocode.ph.unito.it/ The simulation is governed by the ideal MHD conser-vation equations, ∂ρ∂t + ∇ · ( ρ v ) = 0 ,∂ ( ρ v ) ∂t + ∇ · (cid:20) ρ vv − BB π + (cid:18) p + B π (cid:19)(cid:21) T = − ρ ∇ Φ ,∂E t ∂t + ∇ · (cid:20)(cid:18) ρ v ρe + p + ρ Φ (cid:19) v − v × B × B π (cid:21) = − ∂ ( ρ Φ) ∂t ,∂ B ∂t − ∇ × ( v × B ) = 0 , (1)where ρ is the mass density, p the thermal pressure, v the velocity, B the magnetic field, the dyadic tensor,Φ the gravitational potential, and E t the total energydensity, defined as: E t = ρǫ + ( ρ v ) ρ + B π , (2)where ǫ is the internal energy. We use an ideal equationof state, i.e., ǫ = p/ (Γ − µ G for the central 400 pc, based on anupper limit of the detected diffuse γ -ray flux. Given theobserved radio spectral energy distribution of the Galac-tic center, a weaker magnetic field would lead to morerelativistic electrons and consequently a higher γ -rayflux due to inverse Compton emission. In fact, energyequipartition between the magnetic field, X-ray-emittinghot plasma and turbulent gas implies a magnetic fieldstrength of ∼ µ G (Crocker et al. 2010). On the otherhand, Thomas et al. (2020) suggested a stronger mag-netic strength of 200 µ G in the NTFs. In our fiducialrun of simulation, the initial magnetic field strength atthe origin ( x = y = z = 0) is set as B = 80 µ G. Valuesof 50 µ G and 200 µ G are also tested to examine the effectof a weaker/stronger magnetic field (see Section 2.4).
M.F.Zhang et al.
The simulation neglects viscosity and thermal conduc-tion, but takes into account radiative cooling which isapproximated by a power-law cooling function imple-mented in
PLUTO : Λ = a r ρ T α , where a r is a con-stant and α is set as 0.5. Although a more realisticcooling function may be implemented in a tabulatedform, our preliminary test shows that this is rathertime-consuming while having little effect on the dynam-ical evolution of the bubble, since the involved radiativecooling is only moderate. Therefore, we adopt the sim-ple power-law cooling function. We neglect synchrotroncooling of the relativistic electrons, which are presum-ably produced by the SN shocks (see Section 2.4).2.2. Gravitational Potential and Initial ISMConditions
The gravitation in the Galactic center mainly origi-nates from two components, namely, the nuclear starcluster (NSC), which dominates the innermost ∼
20 pc,and the nuclear stellar disk (ND) that occupies the innerfew hundred parsecs. We neglect larger-scale structuressuch as the bar and the Galactic disk. The SMBH, whichhas four million solar masses and a sphere of influenceof a few parsecs in radius, can also be ignored given thescales of interest here. The NSC/ND will not evolve sig-nificantly on the timescale involved in our simulations,hence we adopt a fixed gravitational potential, which,following Stolte et al. (2008), can be approximated by alogarithmic form,Φ = 0 . v log( R c + x a + y b + z c ) , (3)where v is the asymptotic velocity of a flat rotationcurve, R c is the core radius, and a , b and c are stretchingparameters. We adopt v = 98 . − , 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, fromTable 1 of Stolte et al. (2008). The combined NSC+NDpotential has been found to provide a good match to theobserved stellar mass distribution in the Galactic center(Launhardt et al. 2002).At the beginning of the simulation, the ISM is as-sumed to be isothermal and in hydrostatic equilibriumwith the gravitational potential, ∇ Pρ = −∇ Φ , (4)where P = n t kT is the thermal pressure, and n t isthe total number density of gas particles including pro-tons, electrons and heavy elements. As usual we de-fine ρ = µm p n t , where m p is the proton mass and µ ≈ . K, which is roughly the virial temperature given the en-closed gravitational mass of 1 × M ⊙ within 100 pc.The initial density distribution can then be derived bysolving Eqn. 4, as shown in Figure 1 along with the ini-tial magnetic field distribution. From the adopted initialconditions, it can be shown that the thermal pressureof the ISM ( n t kT ∼ − − − dyn cm − ) is ev-erywhere significantly lower than the magnetic pressure( B / π ∼ . × − dyn cm − ), perhaps except inthe innermost few parsecs. In the meantime, the Alfv´enspeed, V A = ( B / πρ ) . km s − , is much lowerthan the typical expansion velocity of the SN. Therefore,the present case of the Galactic center satisfies the mod-erately strong field condition defined by Insertis & Rees(1991). 2.3. Supernova Input
In the simulations, SNe are set to explode within a pre-defined cylindrical volume. The cylinder has a diameterof 50 pc in the x − y plane and a thickness of 10 pc alongthe z -axis, to mimic the concentration of massive starsnear the Galactic plane (Kruijssen et al. 2015). We havetested a wider explosion area in the x − y plane (e.g., 100pc in diameter, closer to that of the CMZ), finding thatthe resultant bubble would become significantly fatter,inconsistent with the observed morphology. In reality,the CMZ may provide a horizontal confinement to thebubble. However, a self-consistent implementation ofthe CMZ would necessarily introduce more free param-eters, and is beyond the scope of the present work. Thebase of the radio bubbles shows a small but appreciableoffset to the west of Sgr A* (Heywood et al. 2019). Thuswe place the center of the cylinder at x = 5 pc to mimicthis behavior. Due to the otherwise axisymmetry in thesimulation, this appears to be the most viable way toreproduce the observed offset.The fiducial SN birth rate is set to be 1 kyr − (Di Teodoro et al. 2018), which is estimated by as-suming an SFR of 0.1 M ⊙ yr − , a Kroupa (2001)initial mass function (IMF) and a minimum mass of8 M ⊙ for the progenitor star of a core-collapse SN.Barnes et al. (2017) and Sormani et al. (2020) estimateda current SFR of 0.1 M ⊙ yr − inside the CMZ, whileNogueras-Lara et al. (2020) found that star formationin the ND (which has a similar radial extent as theCMZ) has been relatively active in the past 30 Myr,with an SFR of 0 . − . ⊙ yr − . Our assumed SFRof 0.1 M ⊙ yr − is compatible with the smaller radialextent of our adopted exploding region, which may bethe case if SN events have been episodic and cluster-ing on a . Myr timescale. We also test the effect ofa lower SN birth rate of 0 . − (see below). We alactic Center Radio bubbles Table 1.
Simulation Parameters for the Radio BubblesFiducial Parameters ValueSN Ejecta Mass 10 M ⊙ SN Kinetic Energy 1 × ergInjection Radius 4 pcAmbient Temperature 1 × KDiameter of Explosion Region 50 pcHeight of Explosion Region 10 pcSimulation Runs B80I1 B80I2 B50I1 B200I1Magnetic Field Strength 80 µ G 80 µ G 50 µ G 200 µ GExplosion Interval 1 kyr 2 kyr 1 kyr 1 kyr have neglected Type Ia SNe, which have a birth rate of . .
05 kyr − according to the enclosed stellar mass inthe ND/NSC (Mannucci et al. 2005), though a recentstudy by Zhou et al. (2020) found evidence that SgrA East, one of the few currently known SNRs in theGalactic center, was created by a Type Iax SN.Individual SNe are thus injected at random positionsinside the cylindrical volume, one after another with afixed interval according to the assumed birth rate. EachSN has an ejecta mass of M ej = 10 M ⊙ and a kinetic en-ergy of E ej = 1 × erg (Poznanski 2013). This energyis deposited into a sphere with a radius of R SN = 4 pc,ignoring any intrinsic anisotropy. The analytic solutionwithin R SN is derived from Truelove & McKee (1999),in which the newly born SN is divided into two parts, theinner uniform density core region and the outer power-law density envelope region. The radius of the former is10 times that of the latter, and the power-law index isset as zero.2.4. Simulation Runs and Synthetic Emission Maps
In this work, we perform four runs of simulation, eachwith a unique combination of magnetic field strengthand SN explosion interval. Our fiducial simulation isrepresented by run
B80I1 , where B and I indicate themagnetic field and explosion interval, respectively. Thefiducial run has B = 80 µ G and I = 1 kyr. The otherthree 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 inferredfor the Galactic center (Section 2.1). Finally,
B80I2 hasan explosion interval of 2 kyr. The total elapsed timeis set to be 330 kyr for all four runs. In the fiducialsimulation, this is about the time when the top of thebubble approaches the edge of the simulation box. Thetime step is adaptive and ranges between 1 −
40 yr. Thesimulation parameters are summarized in Table 1. To facilitate comparison with the observations, wegenerate synthetic radio and X-ray maps for the finalsnapshot (i.e., t = 330 kyr) of the simulation. We includesynchrotron radiation and free-free emission in the ra-dio band (default at 1284 MHz, to be consistent with theMeerKAT observation), while for the X-ray band onlythermal emission from a collisionally-ionized, optically-thin plasma is considered.First we need to distinguish regions inside and outsidethe evolving bubble. This is realized by adding a tracerparameter, Q , evaluated at every pixel in the simulation,which obeys a simple conservation equation: ∂ ( ρQ ) ∂t + ∇ · ( ρQ v ) = 0 . (5) Q has a value of 1 for pure SN ejecta and 0 for theunpolluted ISM, and a value between 0–1 for pixels withmixed ejecta and ISM. We further calculate the Machnumber for every pixel. The synthetic maps only takeinto account pixels with a non-zero tracer parameter ora Mach number greater than 2. The latter condition isto ensure that pixels with a high Mach number but azero tracer parameter, such as those at or immediatelybehind the shock, are included.Synchrotron emissivity depends on the magnetic fieldstrength and the density of relativistic electrons. How-ever, the latter cannot be directly obtained from oursimulation and thus requires some working assump-tion. Here, we assume that the relativistic electrondensity at a given pixel of interest is proportional tothe local gas density (Orlando et al. 2007; Zhang et al.2017), normalized to have a mean energy density of0 . − across the bubble volume. This is com-patible with the estimated mean cosmic-ray energy den-sity of 10 eV cm − in the bubble (Heywood et al. 2019)and the empirical fact that relativistic electrons accountfor ∼
1% of the total cosmic-ray energy density in theGeV band (Blasi 2013). We calculate the synchrotron
M.F.Zhang et al.
100 uG−100 −50 0 50 100x (pc)−100−75−50−250255075100 y ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 ) 100 uG−100 −50 0 50 100x (pc)0255075100125150175200 z ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 ) 100 uG−100 −50 0 50 100y (pc)0255075100125150175200 z ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 ) Figure 1.
Initial distribution of gas density, plotted in logarithmic scale and in a unit of cm − . The white arrows indicate theinitial magnetic field distribution. The three panels are slices through the z = 0, y = 0 and x = 0 planes, respectively. emissivity and integrate along the light-of-sight (i.e., the y -axis) in each pixel to derive the synchrotron intensitymap. In this calculation the y -component of the mag-netic field is neglected due to the nature of synchrotronradiation.Radio free-free emission is calculated following thestandard formula of Longair (2011), which, at a givepixel, scales with density squared and is a functionof temperature. The X-ray emissivity of an optically-thin thermal plasma in collisional ionization equilibrium(Smith et al. 2001), also scaling with density squared, isextracted from ATOMDB , version 3.0.9, for which weadopt a solar abundance. The free-free and X-ray in-tensity maps are again derived by integrating along the y -axis. We find that self-absorption is negligible in boththe radio and X-ray bands, thanks to the relatively lowcolumn density involved. RESULTSIn this section we present the simulation results. Wefirst describe the formation and subsequent evolutionof the bubble in the fiducial run, showing that a goodagreement on the overall morphology of the bubble isachieved between the simulation and observation (Sec-tion 3.1). We then present the other three runs ofsimulations and examine the effect of varying magneticfield strength or SN birth rate on the bubble formation(Section 3.2). Lastly, we confront the synthetic emis-sion maps with the radio and X-ray observations (Sec-tion 3.3).3.1.
Bubble Formation and Evolution in the FiducialSimulation In Figures 2 and 3, we show the gas density and tem-perature maps of run
B80I1 . In each figure, the densityor temperature distribution is shown for a slice throughthe 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 (bottomrows) kyr.By design, 30 SNe have exploded by the time of 30kyr. The forward shock front of several youngest SNeare clearly revealed in the density map of the result, aswell as by the overlaid projected velocity vectors. Ahigh-density region forms and persists around the origin( x = y = z = 0), because of the steep gravitationalpotential even in the presumed absence of an SMBH.As the shocks propagate, they compress and heat theambient gas and also frequently collide with each other,eventually forming an expanding complex of post-shockgas with temperatures of 10 − K.By the time of 180 kyr, this hot gas complex has de-veloped into a bubble structure with a common denseshell, most clearly seen in the x − z and y − z planes.Inside the bubble, the density is low as a result of ex-pansion, while the temperature remains high due to re-peated shock heating. Numerous arc-like features areevident in the temperature map, especially in the x − z and y − z planes, which are the relic of individual SNshocks. At this stage, the bubble looks fat, with asimilar extent ( ∼
100 pc) along the three dimensions.However, the overall expansion starts to show a pref-erence along the vertical (positive z ) direction, withthe vertical expansion velocity of the shell now be-ing ∼
690 km s − , substantially larger than the aver-age expansion velocity of ∼
120 km s − in the x − y plane. This is primarily due to the collimation ef-fect by an ordered magnetic field (Insertis & Rees 1991; alactic Center Radio bubbles
10 km/s−100 −50 0 50 100x (pc)−100−75−50−250255075100 y ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 ) 10 km/s−100 −50 0 50 100x (pc)0255075100125150175200 z ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 ) 10 km/s−100 −50 0 50 100y (pc)0255075100125150175200 z ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 )50 km/s−100 −50 0 50 100x (pc)−100−75−50−250255075100 y ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 ) 500 km/s−100 −50 0 50 100x (pc)0255075100125150175200 z ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 ) 500 km/s−100 −50 0 50 100y (pc)0255075100125150175200 z ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 )500 km/s−100 −50 0 50 100x (pc)−100−75−50−250255075100 y ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 ) 500 km/s−100 −50 0 50 100x (pc)0255075100125150175200 z ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 ) 500 km/s−100 −50 0 50 100y (pc)0255075100125150175200 z ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 ) Figure 2.
Density-velocity distributions after 30 (top row), 180 (middle row) and 330 (bottom row) kyr, for simulation
B80I1 ,i.e., with initial magnetic field strength of 80 µ G and an explosion interval of 1 kyr. The gas density is plotted in logarithmicscale and in a unit of cm − . The white arrows indicate the velocity vector. The left, middle and right columns are slices throughthe z = 0, y = 0 and x = 0 planes, respectively. M.F.Zhang et al.
100 uG−100 −50 0 50 100x (pc)−100−75−50−250255075100 y ( p c ) z ( p c ) z ( p c ) y ( p c ) z ( p c ) z ( p c ) y ( p c ) z ( p c ) z ( p c ) Figure 3.
Temperature-magnetic field distributions after 30 (top row), 180 (middle row) and 330 (bottom row) kyr, forsimulation
B80I1 , The gas temperature is plotted in logarithmic scale and in a unit of Kelvin. The white arrows indicate themagnetic field vector. The left, middle and right columns are slices through the z = 0, y = 0 and x = 0 planes, respectively. alactic Center Radio bubbles −100 −50 0 50 100x (pc)−100−75−50−250255075100 y ( p c )
30 60 90 120B/μG −100 −50 0 50 100x (pc)0255075100125150175200 z ( p c ) z ( p c ) Figure 4.
Magnetic strength distributions after 330 kyr for simulation
B80I1 , The left, middle and right columns are slicesthrough the z = 0, y = 0 and x = 0 planes, respectively. Stone & Norman 1992; Rozyczka & Tenorio-Tagle 1995;Wu & Zhang 2019). Specifically, the SN shocks tendto push the semi-vertical magnetic field to the sides,greatly suppressing the magnetic field inside the bubbleand in the meantime amplifying the magnetic field nearthe bubble shell. In turn, the latter decelerates and evenhalts the horizontal expansion of the bubble. The verti-cal expansion, on the other hand, feels no such magneticconfinement, thus a high velocity along this direction re-mains. The relatively strong gravitational potential inthe x − y plane also contributes to retarding the hori-zontal expansion and facilitates the bubble collimationalong the z -axis.As a result, by the time of 330 kyr, the bubble be-comes much more elongated. The top of the bubblealmost reaches the edge of the simulation box ( z = 200pc), with a vertical expansion velocity still as high as ∼
600 km s − , whereas its horizontal extent has notgrown significantly since t = 180 yr. The width ofthe bubble at its base is about 120 pc, with a smallbut appreciable offset towards the positive x -axis, bothin agreement with the observed bubble. Arc-like fea-tures tracing the sequential SN shocks remain promi-nent throughout the bubble interior. Near some of thesearcs, locally enhanced magnetic fields are evident, whichis the result of shock compression, as illustrated in Fig-ure 4. The magnetic field strength takes a highest valueof 175 µ G across the bubble. Our simulation ends atthis point.3.2.
Comparison with Other Simulation Runs
Similarly, we show the snapshots of density and tem-perature maps of simulation runs
B80I2 , B50I1 and
B200I1 in Figures 5 and 6, all at t = 330 kyr. Thesethree simulations share some common features withthe fiducial simulation. In particular, a vertically- collimated, bubble-like structure is formed in all thesesimulations. The bubble is delineated by a denseouter shell with compressed magnetic field and has alow-density, high-temperature interior with vertically-oriented velocities and generally weak magnetic field.The bubble interior is not smooth, rather, it is filledwith chaotic small-scale structures, again due to thesequential SN shocks and mutual interactions betweenthem. Below we shall describe the more unique featuresin the individual simulations.Simulation B80I2 (top row in Figures 5 and 6) adoptsa lower explosion frequency than the fiducial case. Thisleads to a smaller energy injection rate, but is still suffi-cient 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 thanin the fiducial case (thus also narrower than the observedbubble), which remains the case even if we followed thebubble growth to a height of 200 pc. This occurs be-cause, given a weaker SN energy injection but the samemagnetic confinement, reduction in the horizontal ex-pansion is greater than in the vertical expansion. Wenote that at a further reduced explosion frequency, muchof the SN ejecta would not be able to escape from thestrong gravity near the mid-plane, and a bubble wouldnever form.In
B50I1 (middle row in Figures 5 and 6), which hasa weaker magnetic field compared to the fiducial run,the resultant vertical collimation is less effective andthus the bubble appears fatter. We note that a thinnerbubble could still be achieved, should a lower explosionfrequency be adopted in combination with the weakermagnetic field, for the reason explained above. How-ever, in this case it would take a much longer time forthe bubble to grow to the observed height of 190 pc.0
M.F.Zhang et al.
500 km/s−100 −50 0 50 100x (pc)−100−75−50−250255075100 y ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 ) 500 km/s−100 −50 0 50 100x (pc)0255075100125150175200 z ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 ) 500 km/s−100 −50 0 50 100y (pc)0255075100125150175200 z ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 )500 km/s−100 −50 0 50 100x (pc)−100−75−50−250255075100 y ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 ) 500 km/s−100 −50 0 50 100x (pc)0255075100125150175200 z ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 ) 500 km/s−100 −50 0 50 100y (pc)0255075100125150175200 z ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 )500 km/s−100 −50 0 50 100x (pc)−100−75−50−250255075100 y ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 ) 500 km/s−100 −50 0 50 100x (pc)0255075100125150175200 z ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 ) 500 km/s−100 −50 0 50 100y (pc)0255075100125150175200 z ( p c ) −1.8 −1.2 −0.6 0.0 0.6log(n/cm −3 ) Figure 5.
Simulated density-velocity images after 330 kyr. In the upper, middle and lower rows, we show the results of runs
B80I2 , B50I1 and
B200I1 , respectively. The x − z and y − z panels are slices through the center of the box along each axis,while the x − y panel shows the slice at z = 0. The background is the density distribution with a unit of log(cm − ), and thewhite arrows indicate the velocity. alactic Center Radio bubbles
100 uG−100 −50 0 50 100x (pc)−100−75−50−250255075100 y ( p c ) z ( p c ) z ( p c ) y ( p c ) z ( p c ) z ( p c ) y ( p c ) z ( p c ) z ( p c ) Figure 6.
Simulated temperature-magnetic field images after 330 kyr. In the upper, middle and lower rows, we show the resultsof runs
B80I2 , B50I1 and
B200I1 , respectively. The x − z and y − z panels are slices through the center of the box along eachaxis, while the x − y panel shows the slice at z = 0. The background is the temperature distribution with log scale, and thewhite arrows indicate the magnetic field. M.F.Zhang et al.
In contrast,
B200I1 results in a significantly thinnerstructure. The magnetic field in this run is so strongthat it can resist the compression of the SN shocks andconsequently there is little sweeping of magnetic field in-side the bubble. With the strong magnetic collimation,some SN ejecta are able to rapidly propagate along thefield lines, forming vertical protrusions (several of theseare captured in the x − z and y − z slides). The overallmorphology is obviously inconsistent with the observedbubble. 3.3. Comparison with Observations
Here we shall provide a more quantitative comparisonwith the radio and X-ray observations, with a focus onthe fiducial simulation, which has the best morphologi-cal agreement with the observed bubble.The synthetic synchrotron and free-free intensitymaps 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 thesynchrotron and free-free emission, which is partiallyowing to our assumption that the density of relativisticelectrons scales with the local gas density. However, thesynchrotron intensity is everywhere orders of magnitudehigher than the free-free counterpart in the simulation.This holds true even considering the uncertainties in theenergy density of the relativistic electrons and the mag-netic strength. Consequently, synchrotron dominatesthe total flux density at 1284 MHz, consistent withthe MeerKAT observation (Heywood et al. 2019). It isnoteworthy that both the hydrogen recombination line,H90 α , at 8309 MHz and the 8.4 GHz continuum arefound to trace the GCL (Nagoshi et al. 2019), which ex-hibits a loop-like structure spatially coincident with thenorthern radio bubble. This suggests that the thermalcomponent may have an increasingly larger contributiontoward higher frequencies, which can be due to a com-bined effect of substantial synchrotron cooling at higherfrequencies and the presence of ambient cooler gas nottaken into account in our simulation.The overall extent of the synthetic synchrotron emis-sion 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 190pc. Another interesting feature in the simulation isthe presence of numerous filaments both at the edge ofand inside the bubble, which closely resemble the NTFs(Yusef-Zadeh et al. 1984), although the ones in the sim-ulation appear thicker and fuzzier in general, which maybe partly owing to our moderate resolution. In the simu-lation, these filaments originate from the sequential SNshocks and their mutual interactions, and are associ- ated with locally amplified magnetic field (Figure 4).Their possible relation with the NTFs will be furtheraddressed in Section 4.The 1284 MHz flux density of the simulated bubbleis found to be 5375 Jy, which is to be contrasted withour rough estimate of the observed flux density in theMeerKAT image, 970 Jy, obtained by assuming a meanflux density of 3 mJy beam − across the projected areaof the bubble. We caution that the MeerKAT mosaicimage presented in Heywood et al. (2019) was not cor-rected for the primary beam attenuation and that theextended emission from the bubble suffers from poten-tial flux loss in the interferometric image (I. Heywood,private communication), thus our estimate should betreated as a lower limit of the true flux density. On theother hand, the simulated flux density depends heavilyon the assumed energy density of relativistic electrons.Therefore, the apparently large discrepancy between theobserved and simulated radio flux densities should betaken as a point for future improvement rather than afailure of the simulation.The synthetic 0.5–1.5 keV and 1.5–10 keV X-ray inten-sity maps are shown in Figure 8. Compared to its radiomorphology, the simulated bubble appears smoother inthe 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-rayobservations. This might be due to the fact that the shellis on average cooler than the bubble interior (Figure 3).Indeed, in the 0.5-1.5 keV map, which is more sensitiveto gas temperatures below ∼ cm − ; Ponti et al. 2019). The 1.5-10keV map also exhibits much fewer small-scale structuresin the bubble interior, except near the x − y plane wherethe gas density is high and the most recent SNe freshlydeposit a fraction of their kinetic energy. In particu-lar, remnants of two newly exploded SNe are evidentnear the center (marked in the left panel of Figure 8),although they are not clearly seen in the synthetic ra-dio map. An SNR evolving near Sgr A* will be heavilyshaped by the strong gravity, with a large part of theejecta pulled to the mid-plane, resulting in an appear-ance resembling the bipolar X-ray lobes detected in theinnermost 15 parsecs of the Galactic center (Ponti et al.2015, 2019).The thermal, kinetic and magnetic energy of the bub-ble is calculated by summing over all “bubble pixels”(Section 2.4), which is found to be 1.9, 1.2 and 0.1 × erg, respectively. The initial thermal and magnetic en- alactic Center Radio bubbles −100 −50 0 50 100x (pc)0255075100125150175200 z ( p c ) −4.8 −4.4 −4.0 −3.6 −3.2log S syn (Jy arcsec −2 ) −100 −50 0 50 100x (pc)0255075100125150175200 z ( p c ) −10.0 −9.5 −9.0 −8.5 −8.0log S ff (Jy arcsec −2 ) Figure 7.
B80I1 at 330 kyr. The red dotted line outlines the rim of thenorthern radio bubble.
Left:
Synchrotron emission;
Right:
Free-free emission. −100 −50 0 50 100x (pc)0255075100125150175200 z ( p c ) −8.8 −8.0 −7.2 −6.4log S (erg s −1 cm −2 arcsec −2 ) −100 −50 0 50 100x (pc)0255075100125150175200 z ( p c ) −8.8 −8.4 −8.0 −7.6 −7.2log S (erg s −1 cm −2 arcsec −2 ) Figure 8.
Synthetic 0 . − . . −
10 (right) keV X-ray intensity distribution in simulation
B80I1 at 330 kyr. Thered dotted line outlines the rim of the northern radio bubble, while the black circles highlight two young SNRs. To make iteasier to see some details, we excluded the values lower than 10 − erg s − cm − arcsec − . M.F.Zhang et al. ergy within the bubble volume are 0.7 and 1.1 × erg. A net decrease of the magnetic energy underscoresthe sweep-up of the magnetic field. Ponti et al. (2019)estimated a thermal energy of 4 × erg for the X-raychimneys (sum of the northern and southern halves),which is well matched by the simulated value of 1.9 × erg for the northern chimney.Ponti et al. (2019) also measured density and temper-ature profiles along selected Galactic longitude, l = 0 ◦ ,and Galactic latitude, b = 0 . ◦
7. For a direct comparison,we construct density and temperature profiles at l = 0 ◦ and b = 0 . ◦ l = 0 . ◦ b = − . ◦ x = 0 plane and z = 98 pc plane forcomparison. We calculate the mean density along theline-of-sight calculated as < n > = s R n dV R dV , (6)where n is the number density, dV = Adl , A = 625 pc is the projected area, and the line-of-sight integration( dl ) is from the farthest side to the nearest side of thebubble. We have adopted a constant projected area toapproximate the rather irregular spectral extraction re-gions in Ponti et al. (2019).Similarly, the mean temperature is calculated as < T > = R T n Λ( T, Z ) dV R n Λ( T, Z ) dV , (7)where Λ is the tabulated X-ray emissivity as a func-tion of temperature and metallicity extracted from ATOMDB . By examining the distribution of the SNejecta through the tracer parameter, we have verifiedthat the assumption of a uniform metallicity is a rea-sonable approximation.At l = 0 ◦ , the simulated density profile appears flatwith a value of 0 . − at z &
30 pc, which is in goodagreement with the observed density profile. Below z ∼
30 pc, the simulated density profile modestly increasestoward the mid-plane, while the observed profile risesmuch more steeply. This may be due partly to contam-ination from unresolved stellar objects and non-thermalextended features to the apparently diffuse X-ray emis-sion near the mid-plane (Zhu et al. 2018). The simu-lated temperature profile appears bumpy around a meanvalue of ∼ . z ∼
15 pc, which can be attributed to theheating by the two young SNe shown in Figure 8. Theobserved temperature profile, derived from a spectral fitusing a single-temperature model, appears flatter andhas a lower mean value of ∼ b = 0 . ◦
7, the simulated density profile showsroughly a “U”-shape, with a high value of 0.27(0.42) cm − at the eastern (western) edge, correspond-ing to the position of the bubble shell. The simulateddensity profile is in broad agreement with the obser-vation, except for the eastern edge where the observedX-ray emission shows no sign of limb-brightening. Onepossibility is that the soft X-ray emission from thedenser and cooler eastern shell has largely dropped outof the observation band (above 1.5 keV). The simulatedtemperature profile shows a roughly inverse “U”-shape,with values peaking at ∼ x = 15 pc andreaching a minimum of 0.25 keV at the western edge.The observed temperature profile, again derived froma spectral fit using a single-temperature model, ap-pears flat around a value of 0.8 keV. We note that theobserved density/temperature profiles at this latitudecover a wider range than the simulated profile, becausethe latter only takes into account the bubble volume.Ponti et al. (2019) did not provide an explicit total X-ray luminosity of the chimneys. A rough estimate of thisvalue can be made by adopting a cylinder of 150 pc inboth diameter and height, as assumed by (Ponti et al.2019), a mean density of 0.1 cm − and a mean tempera-ture of 1 × K (0.86 keV), which are representative ofthe X-ray chimneys. This leads to an estimated 1.5–10keV luminosity of ∼ × erg s − for the northernchimney, again well matched by the simulated value of2.0 × erg s − .For completeness, the synthetic radio and X-ray mapsof runs B80I2 , B50I1 and
B200I1 are shown in Fig-ure 10. While these maps exhibit some interesting fea-tures, it is immediately clear that none of them matchesthe observed bubble morphology (again approximatedby the red dotted line). DISCUSSION4.1.
The Origin and Fate of the Galactic Center RadioBubbles/X-ray Chimneys
The simulations presented in the previous sectionshow that an outflow driven by sequential SN explo-sions and collimated by a vertical magnetic field canprovide a reasonable explanation for the observed ra-dio bubbles/X-ray chimneys in the Galactic center. Inparticular, the simulations can well reproduce the over-all morphology, X-ray luminosity and thermal energy ofthe northern bubble.This scenario relies on two key ingredients: SN ex-plosions clustering in the nuclear disk to provide asemi-continuous energy input, and a vertical, moder-ately strong magnetic field to provide the collimation.Both ingredients are very likely available in the Galac- alactic Center Radio bubbles T e m p e r a t u r e ( k e V ) D e n s i t y ( c m − ) sim Tobs Tsim nobs n −100 −50 0 50 100x (pc)0.20.40.60.81.01.2 T e m p e r a t u r e ( k e V ) D e n s i t y ( c m − ) sim Tobs Tsim nobs n Figure 9.
Density and temperature profile of
B80I1 . The light blue dot and pink plus respectively indicate the temperatureand the density in the simulation. The blue dot and red plus respectively indicate the temperature and the density from theobservations. The observed values are manually estimated from Ponti et al. (2019).
Left:
The profile at Galactic longitude l = 0 ◦ . Right:
The profile at Galactic latitude b = 0 . ◦
7. Note that the observed density/temperature profiles cover a wider rangereaching beyond the bubble volume. −100 −50 0 50 100x (pc)0255075100125150175200 z ( p c ) −4 −3 −2 −1log S (Jy) −100 −50 0 50 100x (pc)0255075100125150175200 z ( p c ) −4.0 −3.2 −2.4 −1.6 −0.8log S (Jy) −100 −50 0 50 100x (pc)0255075100125150175200 z ( p c ) −2.4 −1.6 −0.8 0.0log S (Jy)−100 −50 0 50 100x (pc)0255075100125150175200 z ( p c ) z ( p c )
26 28 30 32log L (erg/s) −100 −50 0 50 100x (pc)0255075100125150175200 z ( p c ) Figure 10.
Upper:
Radio flux density at 1284 MHz by synchrotron mechanism.
Lower: ∼
10 keV X-ray luminosity ofthermal emission. The red dot line outlines the morphology of the northern radio bubble. The left, middle and right columnsshow the results of runs
B80I2 , B50I1 and
B200I1 , respectively. M.F.Zhang et al. tic center. Indeed, direct evidence for contemporarySN explosions in the Galactic center was provided byat least a few SNRs clearly visible in radio or X-rayimages (e.g., Ponti et al. 2015). Moreover, about twohundred emission-line objects have been detected in theGalactic center, most of which are likely evolved massivestars (Dong et al. 2012). These stars may belong to thesame population that gave rise to the SNe responsiblefor launching the bubbles. As for the magnetic field,it is widely thought that it is predominantly poloidalin the Galactic center, at least in regions outside thegiant molecular clouds (Ferri`ere 2009). In this regard,an SNe-driven, magnetically-collimated outflow shouldnaturally develop in the Galactic center, provided thecorrectness of our simulations.As mentioned in Section 1, a competing driver of alarge-scale outflow is the kinetic power from the centralSMBH, even though Sgr A* is by no means compara-ble with a classical AGN. While our simulations can-not automatically rule out an AGN-driven outflow, theyshare useful insight on the latter case. Compared tothe distributed SN explosions, energy input from theSMBH is highly concentrated. Thus an AGN-drivenoutflow on the hundred-parsec scale may either acquirea highly elongated shape in the case of a canonical jet-driven outflow (e.g., Zhang & Guo 2020), or inflate afat bubble in the case of a more isotropic wind sym-biotic with the hot accretion flow onto a weakly ac-creting SMBH (Yuan et al. 2015). Magnetic collimationmay also shape the wind-blown bubble, but one expectsthat the resultant structure is again a highly elongatedone. Thus matching the morphology of the radio bub-bles with an AGN wind-driven outflow may require somefine-tuning, which awaits a detailed investigation.We now turn to consider the fate of the radio bub-bles. In the framework of our simulations, the SNe-driven outflow is necessarily an evolving structure. Infact, at the end of our fiducial simulation, the top of thebubble still expands at a speed of ∼
600 km s − (Sec-tion 3.1). Provided a continuous energy injection fromfuture SNe, which is quite likely given the evolved mas-sive stars near the disk plane (Dong et al. 2012), thebubbles should continue to grow and gradually evolveinto a more “chimney”-like structure, as long as a mod-erately strong magnetic field persists to greater heights.Conversely, if SNe were temporarily shut off, one expectsthat the bubble/chimney would ultimately disperse andcollapse within a time not much greater than the sound-crossing time (a few hundred kyr).It is interesting to ask whether the radio bubbles/X-ray chimneys have a causal relation with the Fermi bub-bles (Su et al. 2010) and eROSITA bubbles (Predehl et al. 2020) found on much larger scales. We note that theage of the radio bubbles inferred from our simulationsis only a few hundred kyr, much shorter than the dy-namical timescale of a few Myr originally suggestedby Heywood et al. (2019). However, Heywood et al.(2019)’s estimate was based on the assumption of aconstant expansion velocity of the bubbles, which isimplausible, hence a shorter timescale is expected.The estimated age of the Fermi bubbles, on the otherhand, ranges from 1 Myr (Yang et al. 2013) to 1 Gyr(Crocker & Aharonian 2011). Thus, in the context ofour supernova-based model for the origin of the radiobubbles/chimneys, the radio bubbles would be a dy-namically younger and independent structure simplyevolving in the interior of the Fermi/eROSITA bubbles,which themselves were formed by older activities in theGalactic center.Alternatively, as suggested by Ponti et al. (2019), theX-ray chimney may be a channel that transports en-ergy from the Galactic center to the high-latitude re-gion currently occupied by the Fermi bubbles. In thiscase, the channel should have existed for tens of Myr, sothat star formation in the Galactic center can be suffi-cient to supply the total energy content of the Fermibubbles, ∼ erg (Carretti et al. 2013). However,such a picture contradicts with the capped morphol-ogy of the radio bubbles (though the southern bubbleis not obviously capped in X-rays; Ponti et al. 2021),which, according to our simulations, is naturally ex-plained as the expanding shell of a newly born outflow.This picture may be reconciled if star formation in theGalactic center has been episodic on a timescale of ∼ Origin of the Non-thermal Filaments
The origin of the NTFs has been extensively debatedsince their discovery nearly four decades ago. Proposedmodels for the NTFs include expanding magnetic loops(Heyvaerts et al. 1988), induced electric fields (Benford1988; Morris & Yusef-Zadeh 1989), thermal instabil- alactic Center Radio bubbles
Strength of the Galactic Center Magnetic Field
The magnetic field is a crucial component of theGalactic center environment. At present, the averagefield strength is still quite uncertain. The assumptionof energy equipartition between the magnetic field andrelativistic particles leads to estimates up to ∼ µ G in the more dif-fuse background. Crocker et al. (2010) derived a lowerlimit of ∼ µ G based on the diffuse γ -ray flux andsuggested a typical value of ∼ µ G in the central 400pc region.In our simulation
B50I1 , which adopts a field strengthof 50 µ G, an outflow can be developed, although the re-sultant bubble appears fatter due to the reduced mag-netic confinement compared to the fiducial simulation(Section 3.2). This lends some support to the abovelower limit.On the other hand, simulation
B200I1 , which assumesa field strength of 200 µ G, is obviously inconsistent withthe observation (Figure 10). This conclusion holds evenif the other parameter, the SN birth rate, were adjustedwithin a reasonable range. Qualitatively, at a lower SNbirth rate, the shock and ejecta of individual SNe wouldbe less resistant to the magnetic pressure, thus they areless likely to evolve into a mutual network. The resul-tant outflow hardly takes a bubble shape, rather it wouldconsist of many barrel-like structures, through which in-dividual SN ejecta propagate. Only a much higher SNbirth rate can counteract the magnetic pressure, but thiswould be inconsistent with the currently accepted starformation rate in the Galactic center ( ∼ . ⊙ yr − ).Therefore, our simulations provide a meaningful con-straint on the average magnetic field on 100 pc scales inthe Galactic center, 50 µ G . B . µ G.Our fiducial run
B80I1 demonstrates localized mag-netic field amplification across the bubble, reaching amaximum field strength of 175 µ G. It is expected thatthe global magnetic field would gradually restore to theinitial configuration after the termination of clusteringSN explosion and the dispersion/collapse of the outflow.4.4.
Caveats
Despite the satisfactory reproduction of the major ob-served properties of the radio bubbles/X-ray chimneys,some notable discrepancies exist between our simulationresults and the observations, which warrant the follow-ing remarks remarks.The observed edge-brightened radio bubbles have alow-surface-brightness interior, while in our simulationthe edge-interior contrast is less significant. A possiblecause is that we have ignored synchrotron cooling. Usinga magnetic field of 20 µ G, Heywood et al. (2019) deriveda synchrotron cooling time of 1–2 Myr by assuming that8
M.F.Zhang et al. the electron energy density distribution has a power-lawindex of 2. Based on the same method, we estimate acooling time of 250 kyr for 80 µ G, which is comparableto the evolution time of the bubble in our simulation.Hence the relativistic electrons produced at the earlystage and now filling the bubble interior should be sub-ject to radiative cooling, an effect that is not taken intoaccount but otherwise would enhance the edge-interiorcontrast.An alternative and more likely cause is the absence ofa cool gas shell in our simulation. The presence of coolgas (with a temperature of ∼ K) in the outer partof the GCL has been known for some time (Law 2010;Nakashima et al. 2019). This cool gas is not found inour simulations, owing to the very moderate radiativecooling even in the dense shell of post-shock gas. Thisis also the reason why the free-free emission predicted byour simulation is negligible compared to the synchrotron(Section 3.3). Hence the detected cool gas probably hasan external origin that is missing in the framework ofour simulation. Indeed it is almost certain that cool orcold gas exists in the nuclear disk/CMZ, and part of thisgas may be swept into the bubble shell. For example,Ponti et al. (2021) argued that a gas cloud associatedwith the bright 25 µ m source AFGL5376 has been ac-celerated and is now defining part of the wall of thebubble. An additional source of cool gas is the stellarwind of the massive stars distributed in the nuclear disk.In principle, the Galactic center outflow may also bedriven by stellar winds (Chevalier 1992). Stellar windsas an additional energy and momentum source have notbeen included in our simulation. We can give a roughestimate of the collective energy input from the massivestars in the Galactic center. The stellar winds should bedominated by the Wolf–Rayet stars, which have a typi-cal mass loss rate of 10 − M ⊙ yr − and a wind velocityof 2000 km s − . Thus the ∼