© The Authors 2025

1 Introduction

Almost one and a half decades ago, the all-sky Fermi Large Area Telescope discovered the γ-ray Fermi bubbles (FBs) in the Milky Way halo in the GeV energy range (Su et al. 2010). More recently, the larger X-ray eROSITA bubbles (EBs) were discovered in the 0.6–1 keV regime (Predehl et al. 2020), also stretching into the halo (Liu et al. 2024). As both arise from the Galactic center, it is theorized that they have the same origin (see e.g. Sarkar 2024). The exact driving mechanism that generated a total thermal energy of about 2.6 × 10⁵⁶ erg (Predehl et al. 2020) is, however, still under debate.

The X-ray emission can be explained by bremsstrahlung from rapidly decelerated particles in the shock front of the superbubbles (Yang et al. 2022). The γ-ray emission of the FBs is non-thermal. Cosmic rays (CRs), accelerated in the halo, are necessary to produce the high-energy photons (Su et al. 2010). Different phenomena can provide enough energy for CR acceleration. An active galactic nucleus (AGN) jet, for example, can accelerate CR electrons in the Galactic center and transport them rapidly into the halo where they scatter with the interstellar radiation field by inverse Compton emission (Yang et al. 2022). Regular starbursts or AGN winds can also deliver enough energy for the acceleration of CRs in the Galactic center according to the hadronic wind model. Since CR leptons lose their energy faster than hadrons, only CR protons can reach the halo where they would get scattered by interstellar matter, producing neutral pions (π⁰) that decay into γ-rays, unless CR electron reacceleration in the Galactic halo occurs. The third possibility for CR acceleration in the Galactic halo is the in situ acceleration. Here, like in diffusive shock acceleration, CRs drifting through a magnetized plasma can resonantly generate magnetohydrodynamic waves via a so-called streaming instability, which leads to a strong scattering and a secular gain of CR energy after multiple reflections across the converging flow of a shock wave (Axford et al. 1977; Bell 1978; Blandford & Ostriker 1978). Another kind of in situ acceleration can be provided by the shock drift mechanism, resulting from gradient drift across the shock and to an E × B-electric field component, which is strongest for perpendicular shocks, to accelerate particles at a high rate (Jokipii 1982).

Hadronic models appear unlikely because both their spectrum of the π⁰ decayed γ-rays and their synchrotron radiation spectrum of the secondary leptons are too soft (Su et al. 2010; Ackermann et al. 2014). High altitude water Cherenkov (HAWC) observations of γ-rays in the northern FBs also disfavor hadronic models since their flux of high energy γ-rays is too low (Abeysekara et al. 2017).

Many numerical studies are able to reproduce the γ-ray FBs (e.g. Yang et al. 2012; Guo & Mathews 2012; Crocker et al. 2015; Sarkar 2019; Ko et al. 2020; Zhang & Guo 2020). The X-ray EBs, however, have a different emission energy, shape, size, and contain ten times more thermal energy than the FBs, making it necessary to significantly adjust all FB models. New approaches must also be able to explain the two different bubbles simultaneously. Due to these challenges, to our knowledge only Yang et al. (2022), Chang & Kiang (2024), and Tseng et al. (2024) were able to numerically reproduce both the EBs and FBs simultaneously.

All of these studies assume that the supermassive black hole (SMBH) Sagittarius A^∗ (Sgr A^∗) in the Galactic center underwent an active phase with near-Eddington luminosity in which it ejected relativistic jets into the Milky Way halo. The AGN jets within the models start a few million years ago and last from several hundred thousands to a few million years. After the active phase of the jets, they expand into an X-ray emitting spherical cocoon shape with γ-ray emitting CRs in the interior, thus resembling the EBs and FBs simultaneously.

Although AGN jet models are able to explain both the EBs and FBs at the same time, the present weakness of the radio source at Sgr A^∗ does not indicate an active phase of the central SMBH only a few million years ago. Furthermore, most AGN jet models assume that the jets are orthogonal to the Milky Way disk. Thus, the accretion disk of Sgr A^∗ would have needed to be co-planar with the Milky Way disk, since AGN jets are ejected perpendicular to the accretion disk of their SMBH. However, the Event Horizon Telescope Collaboration (2022) discovered that the disk of Sgr A^∗ has an inclination to the line of sight from the Earth and hence to the Milky Way disk. Therefore, relativistic jets of Sgr A^∗ should have been tilted to the Galactic disk. Numerical investigations by Sarkar et al. (2023) show that tilted jets can produce spherical bubbles only if the jets are strongly dissipated. Dutta et al. (2024) use numerical models of clumpy media to rule out dissipation in the clumpy Milky Way halo, meaning that the dissipation had to occur in the central molecular zone (CMZ) or by magnetically dominated jets (see also Sarkar 2024, and references therein). Tseng et al. (2024) can reproduce the EBs and FBs with tilted, CMZ dissipated jets, however, Dutta et al. (2024) argue that the imprint of recently dissipated jets would be visible in the CMZ today and that magnetic dissipation should have resulted in a stronger forward shock than observed in the EBs. Furthermore, Mou et al. (2023) find that tilted AGN jets cannot reproduce the asymmetric brightness distribution of the EBs. Thus, future studies need to explore in greater detail whether the EBs and FBs could still originate from dissipated jets.

Stellar feedback, especially from supernovae (SNe) could also reach luminosities high enough to lead to a thermal outflow, as seen in the M82 galaxy (Rieke et al. 1980). The CMZ in the Milky Way, however, is assumed to have a much lower nuclear SN rate of 0.02–0.15 per century (Crocker et al. 2011; Ponti et al. 2015). With an average SN energy of 10⁵¹ erg, a maximum luminosity of 4.8 × 10⁴⁰ erg s^–1 can be achieved. This luminosity could blow the FBs in a few tens of millions of years, which is why different models consider star formation as their origin (e.g. Crocker et al. 2015; Sarkar 2019). The total thermal energy of the EBs, however, would be reached after a minimum of about 170 Myr. This time scale is significantly longer than the predicted current age of the EBs that range from a few Myr to a maximum of 20 Myr (Predehl et al. 2020; Schulreich & Breitschwerdt 2022; Yang et al. 2022; Mou et al. 2023). Similar to Mou et al. (2023), who argue that SNe in the CMZ lose too much energy to drive a thermal outflow due to radiative cooling, we conclude that SNe alone are most likely not the origin of the EBs.

On a different note, SMBHs such as Sgr A^∗ are assumed to tidally disrupt nearby stars on a regular basis due to their high differential gravitational forces (Lacy et al. 1982; Cheng et al. 2011). Only about half of the debris is accreted to the SMBH, so that tidal disruption events (TDEs) typically release large amounts of energy, scaling with the stellar and SMBH masses (Alexander 2005). In galaxies similar to the Milky Way, TDEs are expected to occur every 10–100 kyr (Magorrian & Tremaine 1999). Numerical simulations by Ko et al. (2020) show that the energy of regular TDEs is sufficient to blow superbubbles into the Milky Way halo, assuming a luminosity of 3 × 10⁴¹ erg s⁻¹. However, they only studied the FBs in simplified Milky Way halo models and did not consider the larger EBs.

In this paper, we perform numerical studies of the TDE model in which γ-rays are accelerated in situ, and X-rays result from a strong forward propagating shock, strengthening in a halo with decreasing density, thereby compressing and heating the gas. In Sect. 2, we describe our numerical setup. Section 3 shows the simulation results and compares the derived synthetic X-ray maps with observations. We summarize the main results in Sect. 4.

2 Methods

Hydrodynamical simulations of the TDE model are performed using the publicly available adaptive mesh refinement (AMR) code RAMSES¹ (Teyssier 2002). The code uses the second-order accurate monotonic upstream-centered scheme for conservation laws (MUSCL), which is an extension of Godunov’s method (Toro 1997) to solve the Euler equations with the aid of an approximate Harten-Lax-van Leer-contact (HLLC) Riemann solver (Toro 1992; Toro et al. 1994). To model the evolution of the superbubbles as realistically as possible, three-dimensional simulations are used. The simulation box is chosen large enough so that the bubbles never reach the boundaries during their evolution.

In a resolution study, we investigated the necessary refinement of the AMR grid to reduce numerical diffusion (Patankar 1980). The refinement level of each cell was determined by the pressure and density gradient to its neighboring cells to follow the shocks and contact discontinuities in great detail. We found that the shape of the superbubbles converges at a maximum refinement level of ∼9.8 pc. The minimum resolution in the AMR grid was set to 312.5 pc.

The energy and mass outflows from TDEs occur on a much smaller scale than the resolution of our simulations since stars are tidally disrupted within the Roche radius of an SMBH, which is smaller than a parsec for the SMBH at the Galactic center (Lacy et al. 1982). Thus, we modeled TDEs as simple energy and mass injections that are homogeneously distributed in the central eight grid cells. The edge length of the cubic injection region (in which the gas speed was set to be zero) was therefore ∼20 pc, whereas the entire cubic computational domain had an extension of 40 kpc. Given that TDEs in the Milky Way occur every 10–100 kyr on average, we injected energy and mass either every 10 or 100 kyr, to probe the two limiting cases. The energy of a single TDE, either 9.5 × 10⁵² or 9.5 × 10⁵³ erg in our model, was injected in a single time step, resulting in the same average luminosity of 3 × 10⁴¹ erg s⁻¹ in both cases. The chosen luminosity was based on energy estimates for the EBs (Predehl et al. 2020) and TDE events in the Milky Way (Ko et al. 2020). During our numerical investigations, we found that the injected mass hardly influences the shape of the bubbles but was necessary for preventing the interior density, and thus, due to our refinement criterion, also the simulation time step, from approaching progressively smaller values. We chose to set the injected mass to the respective initial density in the central cells of the selected Milky Way model.

The gas in our simulation obeys the ideal gas law. It is initially in an isothermal state. We do not consider self-gravity, but instead assume a hydrostatic equilibrium at the beginning of our simulations that is established by an external gravitational field counteracting the thermal pressure gradient. This assumption is necessary to keep the background halo at rest during the simulations. Gravitation as an external force field does not directly affect the dynamics of the bubbles, since its potential energy is much smaller than the energies of the TDEs and superbubbles, but it shapes the density distribution in the halo via the hydrostatic equilibrium, which indirectly feeds back on the bubble dynamics. The density distribution ρ(R, ɀ), dependent on the height ɀ above the Galactic disk and the distance R in radial disk direction, is then related to the gravitational potential Φ(R, ɀ) by $$\rho (R,ɀ)={\rho }_0\exp \left(\frac{{\text{Φ}}_0-\text{Φ}(R,ɀ)}{k_{\text{B}}T/\left(\mu m_{\text{H}}\right)}\right),$$(1)

with ρ₀, Φ₀, k_B, T, and µ = 0.62 being the the central density and gravitational potential, Boltzmann’s constant, the temperature, and the mean molecular weight for an ionized gas with solar abundances in units of the hydrogen atomic mass m_H, respectively.

The background density distribution was found to have a strong influence on the structure of the evolving superbubbles. Therefore, different axisymmetric Milky Way mass models were considered.

First, we used a simple exponential halo model (e.g. Cordes et al. 1991; Nakashima et al. 2018) with a distribution of $$\rho (R,ɀ)=n_0\mu m_{\text{H}}\exp \left(-\frac{R}{R_0}\right)\exp \left(-\frac{ɀ}{{ɀ}_0}\right).$$(2)

Its parameters include the central number density n₀ = 0.004cm⁻³, the scale radius R₀ = 7 kpc, and the scale height ɀ₀ = 2 kpc (Nakashima et al. 2018). The temperature chosen for this halo model was T = 10⁶ K. Equation (2) is plotted in the left panel of Fig. 1.

Secondly, the β-model halo proposed by Miller & Bregman (2013) was considered. It is based on O VII and O VIII absorption line observations and assumes a modified King profile, which shows good agreement with the density distribution of other galaxies (O’Sullivan et al. 2003). It is given by $$\rho (R,ɀ)=n_0\mu m_{\text{H}}{\left[1+{\left(\frac{R}{R_0}\right)}^2+{\left(\frac{ɀ}{{ɀ}_0}\right)}^2\right]}^{-3\beta /2},$$(3)

with n₀ = 0.46cm⁻³, R₀ = 420 pc, ɀ₀ = 260 pc, and a galaxydependent factor, constrained to β = 0.71 for the Milky Way. The corresponding plot is shown in the central panel of Fig. 1. As for the exponential model, the temperature was set to T = 10⁶ K.

Finally, the multicomponent model of McMillan (2017) was investigated, according to which the total gravitational potential of the Milky Way can be assembled by summing up six individual components. Only the density distribution obtained via the Poisson equation is shown here since the corresponding gravitational potential cannot be described analytically. The density distribution of the hydrostatic equilibrium needed to be calculated from the numerically given gravitational potential via Eq. (1) with ρ₀ = 0.03 cm⁻³ µmH (Zhang & Guo 2020) and T = 2.32 × 10⁶ K (corresponding to a thermal energy of 0.2 keV). The first two components are the thin and thick stellar disks that are described by $${\rho }_{\text{d}}(R,ɀ)=\frac{{\text{Σ}}_0}{2{ɀ}_{\text{d}}}\exp \left(-\frac{\left|ɀ\right|}{{ɀ}_{\text{d}}}-\frac{R}{R_{\text{d}}}\right),$$(4)

where Σ₀ = 896 M_⊙ pc⁻² is the central surface density, ɀ_d = 300 pc is the scale height, and R_d = 2.5 kpc is the scale length for the thin disk, and Σ₀ = 183 M_⊙ pc⁻², ɀ_d = 900 pc, and R_d = 3.02 kpc for the thick disk. The next components are the H I and H₂ gas disks which follow $${\rho }_{\text{g}}(R,ɀ)=\frac{{\text{Σ}}_0}{4{ɀ}_{\text{g}}}\exp \left(\frac{R_{\text{m}}}{R}-\frac{R}{R_{\text{g}}}\right){\mathrm{sech}}^2\left(\frac{ɀ}{2{ɀ}_{\text{g}}}\right),$$(5)

with Σ₀ = 53.1 M_⊙ pc⁻², ɀ_g = 85 pc, R_g = 7 kpc, and R_m = 4kpc for the H I disk, and Σ₀ = 2180 M_⊙ pc⁻², ɀ_g = 45 pc, R_g = 1.5 kpc, and R_m = 12 kpc for the H₂ disk, respectively. The bulge of the Milky Way was modeled as $${\rho }_{\text{b}}(R,ɀ)=\frac{{\rho }_{0, \text{b}}}{{\left(1+\sqrt{R^2+4{ɀ}^2}/r_0\right)}^{\alpha }}\exp \left[-{\left(\frac{\sqrt{R^2+4{ɀ}^2}}{r_{\text{cut}}}\right)}^2\right],$$(6)

with ρ_0,b = 9.93 × 10¹⁰ M_⊙ kpc⁻³, r₀ = 0.075 kpc, α = 1.8, and r_cut = 2.1 kpc. The last component is the dark matter halo which is described by $${\rho }_{\text{h}}(R,ɀ)=\frac{{\rho }_{0, \text{h}}}{\left(r/r_{\text{h}}\right){\left(1+r/r_{\text{h}}\right)}^2},$$(7)

where $r=\sqrt{R^2+{ɀ}^2},{\rho }_{0, \text{h}}=0.00854$ M_⊙ pc⁻³, and r_h = 19.6 kpc. The density distribution for this model is shown in the right panel of Fig. 1.

All three Milky Way density models were studied in separate numerical simulations under the same assumptions. To keep the computational effort low, radiative cooling is not included in our model. It would lead to a thinner and denser shell and would have an impact on local structures but would not strongly influence the global shape of the superbubbles. Furthermore, magnetic fields are neglected. An initially disk-parallel magnetic field would wrap around the bubbles so that it would be tangential to the contact discontinuity, thereby strongly inhibiting Rayleigh-Taylor (RT) instabilities (Breitschwerdt et al. 2000). Furthermore, magnetic fields would decelerate the expansion due to magnetic tension forces. Our simulations should therefore be taken as a lower limit with respect to the dynamical time scale. Significant changes to the shape of the bubbles are not expected, since the energy input by TDEs vastly exceeds the magnetic energy stored in the disk field. As bubble expansion progresses, and the ambient density decreases (see above), the frozen-in magnetic field diminishes accordingly.

From the most promising simulation results, synthetic eROSITA maps were calculated, showing the expected X-ray emission of the simulated superbubbles. A collisional ionization equilibrium source model for X-ray emission of the intergalactic medium (IGM) is considered, using the pyXSIM² implementation (ZuHone et al. 2014) of the Cloudy algorithm (Chatzikos et al. 2023). For the cooler, less relevant regimes outside of the main shock, a different IGM model including photoionization is used (Khabibullin & Churazov 2019). Absorption was simulated using the Tuebingen-Boulder ISM absorption model (Wilms et al. 2000), a variable foreground column density NH, based on data from the HI4PI Collaboration (2016), and the abundance table by Feldman (1992) using a constant solar metallicity. Different metal- licities do not change the general X-ray bubble structure but shift the intensity distribution slightly. With a lower metallicity, which is expected in the Milky Way halo, especially the background emission becomes more significant. Of the initially generated X-ray photons that reach the internal observer positioned in the Solar System, only a random sample was projected into an all-sky map using pyXSIM. Then, the instrumental effects were replicated with the Simulated Observations of X-ray Sources (SOXS)³ (ZuHone et al. 2023) software package, using the response matrix files, ancillary response files, point-spread functions, and particle background file of eROSITA, distributed with the Simulation of X-ray Telescopes (SIXTE)⁴ software (Dauser et al. 2019). We also assumed that the filter of every eROSITA camera was active, as this is the standard observing mode (Predehl et al. 2021). For the exposure time within the all-sky survey, we used a constant value, large enough to see the significant features in the final maps, because the reproduction of the varying exposure time of eROSITA (Predehl et al. 2021) is beyond the scope of this paper. Lastly, only photons with energies between 0.6 keV and 1 keV were selected for the synthetic X-ray maps using the Python implementation healpy (Zonca et al. 2019) of HEALPix⁵ (Górski et al. 2005).

Fig. 1 Density distributions of the different Milky Way models dependent on the height ɀ and the radius R away from the Galactic center. Two-dimensional slices through the origins of the three-dimensional simulation boxes are shown. The exponential halo model (left), the β-model halo (middle), and the multicomponent model (right) are described by Eqs. (2), (3), and (4–7), respectively. The models differ by several orders of magnitude in their density distribution.

3 Results

To start our analysis, we first performed simulations in a simple homogeneous ambient medium and compared it to the corresponding analytical solution for a wind-blown bubble (Castor et al. 1975; Weaver et al. 1977). Our numerical results were in very good agreement with the theory. Next, we included a simplified exponential and a simplified β-model halo in our simulation, which can be analytically solved using the Kompaneets formalism (Kompaneets 1960; Baumgartner & Breitschwerdt 2013; Ko et al. 2020; Schulreich & Breitschwerdt 2022). Although the ellipsoidal structure of the exponential halo and the more spherical structure of the β-model halo bubbles were successfully reproduced in these simulations, the numerically resulting bubbles were smaller than the theoretical predictions at any given time step. This deviation occurred because we included fewer assumptions in our model than Kompaneets (1960). For example, we accounted for the surrounding pressure, which led to a noticeable deceleration and resulted in smaller bubble sizes in our model. Therefore, we can conclude that our numerical model reproduced the analytically given features. Hence, more complex and realistic Milky Way models could be implemented into the model.

3.1 Exponential halo model

In the exponential halo model, mass and energy injections were initially set to happen at each time step instead of every 10 or 100 kyr to study the general shape of the bubble, which was not significantly altered by changing the TDE rate with similar luminosity. The time step varied so that energy and mass were injected every 100–1000 yr. The simulation was stopped at 12.5 Myr because the approximate present-day size of the EBs was reached. In total, an energy of 1.2 × 10⁵⁶ erg was injected.

The result can be seen in the top-left panel of Fig. 2. Although the forward shock agrees well with the spherical shape of the EBs, the shell and especially the contact discontinuity of the superbubbles are strongly eroded by fast growing RT instabilities, which arise because the density gradient of the exponential halo is opposite to the shock acceleration (Rayleigh 1882; Taylor 1950).

Following the in situ model for CR acceleration, the contact discontinuity should resemble the FBs in our model. After the CR protons and electrons gain energy through first-order Fermi acceleration at the strong termination shock (ℳ > 20, with ℳ = v/c, where v is the local flow speed and c the local sound speed) visible in the inner 3–5 pc, and at individual shocks propagating into the halo due to consecutive TDEs, they can penetrate the low-density interior of the bubbles. The γ-ray emission could then come from inverse Compton scattering of high energy electrons on the cosmic microwave background (or other low-energy) photons, synchrotron emission from electrons and positrons in strong magnetic fields, CR proton interactions with low-energy protons leading to π⁰-decay, and bremsstrahlung.

Through the strong erosion of the contact discontinuity in the exponential halo model bubble, this γ-ray emission would certainly indicate instabilities in observations even if small perturbations were smoothed out. However, Fermi suggests that the FBs do not incorporate large instabilities. Similarly, eROSITA finds fairly smooth edges of the EBs, not matching the strong deformations seen in this bubble model. This model therefore deviated from the observational EBs and FBs enough to not be considered for further investigations with more realistic TDE rates. However, we note that above a cut-off wavenumber, the strong RT instabilities could be inhibited by magnetic fields (Breitschwerdt et al. 2000), which are not considered in our model.

Fig. 2 Resulting density distributions of the numerical simulations for different Milky Way models and TDE rates. The snapshots are taken at the end of the simulations when the forward shocks reach the size of the EBs and show the detailed bubble structure, including the elapsed time since the first TDE. The different subfigures correspond to the results of the exponential halo model with energy and mass injections every 100–1000 yr (top left), the β-model halo with TDEs every 10 kyr (top right), the β-model halo with TDEs every 100 kyr (bottom left), and the multicomponent model with TDEs every 100 kyr (bottom right).

3.2 β-model halo

For the β-model halo, more realistic TDE rates of 10⁻⁴ yr⁻¹ and 10⁻⁵ yr⁻¹ were considered. The numerical result of the higher TDE rate is shown in the top-right panel of Fig. 2. In a comparison between the forward shock of the simulated bubbles and the X-ray EBs of Predehl et al. (2020), one can see that the spherical and symmetric shape in the northern and southern Milky Way halo is well matched. The shell has fewer and weaker RT instabilities than the exponential halo model bubble because its density gradient is not as strong. It therefore resembles the observations with unperturbed edges better.

Furthermore, the bubbles grow for 16 Myr until they reach the size of the EBs in a β-model halo, leading to a total injected energy of 1.5 × 10⁵⁶ erg. A comparison to the actual evolution time of the EBs is difficult since it is highly uncertain. Predictions range from a few Myr for an AGN jet model (Yang et al. 2022) to ∼20 Myr for either a different AGN jet model with a circumgalactic medium wind (Mou et al. 2023), or theoretical calculations with today’s shock velocity (Predehl et al. 2020) and the Kompaneets formalism (Schulreich & Breitschwerdt 2022). The evolution time of our model agrees well with the upper limit of this range. We note that TDEs can also happen before the derived evolution time in our model as the lifetime of Sgr A^∗ is much longer. However, as stars typically gather in clusters, TDEs often occur concurrently. Thus, a quiescent phase with no Galactic outflow and thus no bubbles, followed by an episode of successive TDEs could explain why the EBs and FBs only started to be driven a few tens of Myr ago, and were not present before. Since the real times of past TDEs are unknown, we chose to model the simplest possibility, meaning that a TDE quiet period ended 16 Myr ago and was followed by a phase of continuous TDEs.

The higher TDE rate β-model halo bubbles still had potential for improvement because they exhibited quite strong RT instabilities at the contact discontinuity that almost penetrated the shell although being weaker than in the exponential halo model. This would result in γ-ray emission at very high latitudes, different to the FBs which only reach up to ±50^◦ to ± 60^◦. Thus, the model could still be improved.

The bottom-left panel of Fig. 2 shows the results of the lower realistic TDE rate (10⁻⁵ yr⁻¹ ) simulation. No strong deformation of the contact discontinuity by RT instabilities is visible anymore, since there are fewer but stronger events. The forward shock shape hardly changes because the total injected energy remains the same.

To better compare the two-dimensional slices to the eROSITA and Fermi all-sky surveys, we calculated synthetic X-ray maps from the three-dimensional simulation results most similar to the EBs and FBs and plotted them in a Hammer-Aitoff projection using Galactic coordinates, analogous to Predehl et al. (2020). Because our model did not include background X-ray radiation, we added it as a constant value, so that the minimum brightness arrived at 4.66 photons s⁻¹ deg⁻². This value represents the mean of the observed brightness profiles (Predehl et al. 2020, Fig. 2b) above α_G = 90^◦ and below α_G = −90^◦. The map of the β-model halo bubble with a TDE rate of 10⁻⁵ yr⁻¹ is shown in the top panel of Fig. 3. Since our model does not consider CRs, we are not able to obtain synthetic γ-ray observations. We do, however, argue that the γ-ray emission should occur at the contact discontinuity of the simulated bubbles. This would appear in Fig. 3 at the inner edge of the brighter (pink-white) region. The EB and FB observations of Predehl et al. (2020) and Su et al. (2010) are indicated with red and yellow data points, respectively. One can see that the overall structure is well matched but the simulated forward shock and especially the contact discontinuity are a bit too wide.

There are several reasons for this. First of all, the actual times of past TDEs in the Milky Way are unknown. There could have occurred multiple TDEs concurrently or subsequently, or no TDEs for a longer period of time, altering the evolution and possibly the shape of the contact discontinuity. Second, the energy and mass of the past TDEs are unknown. If the TDEs had an higher overall luminosity than assumed, the bubbles would need less time to spread into the halo, possibly leading to a narrower forward shock and contact discontinuity, better resembling both superbubbles. Third, the Milky Way temperature and density distribution are highly uncertain. We try to overcome this by considering different Milky Way models, however, all models still may be too simplistic to match the actual Milky Way density distribution, lacking clumps or other inhomogeneities that can alter the forward shock evolution. Fourth, although neglecting cooling and the effects of magnetic fields may be justified for calculating the overall structure of the bubbles, their inclusion will certainly have an influence locally and on the detailed structure. Taking all of these limitations into account, the simulation results match the observations fairly well.

This is supported by the brightness profiles of our synthetic X-ray maps, shown in the top panel of Fig. 4. At latitudes ±40^◦, ±50^◦, and ±60^◦, these can be compared to the eROSITA data in Predehl et al. (2020, Fig. 2b), which are shown in Fig. 4 in red and blue for the northern and southern EB, respectively. As our bubbles are almost perfectly symmetric with respect to the Galactic plane, we refrain from showing the southern latitudes separately. The background radiation was again added as a constant value, similarly to Fig. 3. The observed northern EB data was furthermore shifted by 16.5^◦ to the east to account for the observed EB shift. In a comparison of the brightness profiles, our simulations and the observations agree fairly well. At all shown latitudes, the brightness profiles of our synthetic β-model halo bubbles and especially the northern EB have peaks at ±40^◦ longitude. Since our model does not consider any asymmetries, the eastern peak of the observations is brighter and wider than in our maps because it is influenced by the very bright North Polar Spur, which lies on the outer edge of the EBs. This is also why the bright region in the east appears to be more on the outside in the observations. The western peak, without North Polar Spur, is slightly dimmer but similarly wide compared to our model. A decrease in brightness between those peaks can be seen in the northern EB and in our data. As the observed southern EB is not very bright, a detailed comparison to it is more difficult. The brightness increase of the southern EB, however, still seems to be at similar longitudes as the northern EB and our synthetic observations. Compared to the brightness profiles of synthetic X-ray maps of other studies that use an AGN jet model, such as Yang et al. (2022), Mou et al. (2023), Chang & Kiang (2024) and Tseng et al. (2024), our brightness profiles seem to match fairly well the eROSITA data, especially in terms of peak visibility and position. Of course, we note that all numerical models still had many initial assumptions that have a significant impact on the brightness.

Since there is no termination shock visible in the resulting bubbles (bottom-left panel of Fig. 2) due to the low TDE rate, we furthermore wanted to study whether the TDEs lead to a shock at all. Otherwise, the γ-ray emission of the FBs may be left unexplained as CRs could not be sufficiently accelerated. To study this in detail, we applied a shock tracer algorithm. Following a similar approach as Schaal & Springel (2015) and Pfrommer et al. (2017), we found shocked cells by checking if

• (i)

  ∇ ·v < 0,

• (ii)

  ∇T · ∇ρ > 0, and

• (iii)

  ℳ > 1.0.

These conditions were checked for a central simulation box slice, also considering the adjacent slices to allow for a correct evaluation of the finite differences in all spatial directions. Figure A.1 (provided in the appendix) shows the detected shocked cells for the TDE shock evolution marked in gray. The snapshots were taken from 18–18.1 Myr, slightly after the observed size of the EBs had been reached.

One can see in Fig. A.1a that the newly injected TDE creates a shock front in the very center of the slice. In addition, ellipsoidal shock fronts from previous TDEs are still visible, extending to at least 4 kpc in the slice shown. The new TDE shock then propagates rapidly through the interior due to its high energy and therefore catches up and collides with a slower pre-existing shock in Fig. A.1b. As a result of the collision, the TDE shock splits into a transmitted shock, which continues to propagate toward the shell, and a reverse shock, which moves back to the center. In Figs. A.1c–e, the forward shock continues on its path only slower, whereas the reverse shock reaches the center in Fig. A.1e and reverses a second time. In Fig. A.1f, both shocks still move outwards, while a new TDE shock is already traced in the center of the box because 100 kyr have passed since the last TDE.

Although we showed that there are shocks traveling through the interior of the bubble, their Mach number of 𝓜 = 2–3 (Fig. 5, top panel) might be too low to accelerate CRs to high enough energies. However, our model assumed one of the worst cases for building high-Mach TDE shocks, since the TDEs all had the same energy, and happened only every 100 kyr. Therefore, stronger shocks due to occasional higher energy TDEs or consecutive TDEs from stellar clusters were never encountered in our model. Furthermore, once the bubbles have expanded, it is to be expected that shock waves due to ordinary SN explosions will easily travel out into the halo, catching up with previous shocks, reinforcing them, and strengthening as they travel down a density gradient, as has been proposed in the context of CR-driven galactic winds (Dorfi & Breitschwerdt 2012).

Figure 5 also shows the velocity of the shocks. The TDE shocks reach a vertical speed of ~10⁴ km s⁻¹, while the velocity is an order of magnitude lower for the forward shock. The recently injected TDE has effects on all fluid quantities shown and can therefore be easily identified in all panels at ɀ = 0 kpc. Other properties of our results show that the interior has a low density, a high temperature, and is almost isobaric, whereas the shell has a high density, a high Mach number, and a lower temperature than the interior, perfectly aligning with typical properties of superbubbles.

For a further comparison to the EBs, the shell thickness of the simulated superbubble, which is the separation distance between the forward shock (dashed orange line) and the contact discontinuity (dashed green line), can be determined from the density, Mach number, and temperature profile to be ≳1 kpc. Predehl et al. (2020) estimate the EB shell to be ~2 kpc thick by comparing geometric models to surface brightness maps at different Galactic latitudes. The synthetic observations of Fig. 3 also support the conclusion that the shell is too thin in our model and thus, the FBs are too large. However, as stated earlier, we have not considered magnetic fields in our simulations, which would certainly increase the shell thickness because magnetic pressure would reduce the overall compression. The real magnetic field in the Milky Way and especially the Galactic center is not known with certainty as it is rather complex (see review of Beck & Wielebinski 2013). However, poloidal components (e.g. Ferrière 2009) could be swept up by the TDE shocks in our model, leading to a thicker shell in both directions, parallel and perpendicular to the Galactic plane, which would resemble the FBs better. Other changes such as a slightly different density distribution or a non-constant TDE rate could also impact the shell thickness significantly. We leave the detailed investigation of the impact of these different parameters for a future TDE study.

To summarize, we find that our β-model halo bubbles reproduce the main features of the EBs and FBs fairly well, given the numerous unknown and uncertain initial parameters. Furthermore, we traced TDE shocks in the interior that could explain the CR γ-ray emission. However, the FB shape does not perfectly match this model. An improvement can be made by including magnetic fields that would increase the shell thickness, or occasional high-energy or clustered TDEs that would affect the inner structure significantly. Furthermore, more realistic simulations could also include short-lived collimated jets that are observed for some TDEs (e.g. Swift 1644+57/GRB 110328, see Bloom et al. 2011; Shao et al. 2011). These jets would induce a direction to the rather omnidirect unbound stellar debris of the TDEs and thus lead to a more ellipsoidal FB shape.

Fig. 3 Synthetic observations of the full three-dimensional simulation results using a Hammer–Aitoff projection in Galactic coordinates. The color-coded intensity is the surface brightness B in logarithmic units for the photons detected per time and solid angle and corresponds to detected X-ray photons with an energy between 0.6 keV and 1 keV. Since the cosmic X-ray background was not included in our model, we added it as a constant value to the maps so that we arrive at the same minimum brightness of 4.66 photons s−1 deg−2 as Predehl et al. (2020). The colorbar limits (obtained from Tseng et al. 2024) are the same as in the observations of Predehl et al. (2020). The main underlying model is the thermal collisional ionization equilibrium model of Cloudy. Both, the geometry for an observation from the Solar System and instrumental effects of the eROSITA telescope were considered for these figures. The FBs would appear at the contact discontinuity, which lies at the inner side of the brighter region (pink-white color) according to our model. The approximate shell of the observed EBs (red) and FBs (yellow) are indicated using data points obtained from Predehl et al. (2020, Fig. 3). We note that these points have a reading error of about ±5° and ±1 h because the edges of the observed EBs and FBs are not clearly defined everywhere. The top panel shows the β-model halo bubble (cf. Fig. 2 bottom left) and the bottom panel shows the multicomponent model halo bubble (cf. Fig. 2 bottom right).

Fig. 4 Brightness profiles of our synthetic observations at different latitudes. The top panel shows the surface brightness B for the β-model halo bubble and the bottom panel for the multicomponent model bubble for the latitudes ±40°, ±50°, and ±60° along the longitudes αG. Only the northern brightness profiles of our synthetic observations (black) are shown and compared to the data of the northern (red) and southern (blue) EB from Predehl et al. (2020). The northern eROSITA data was shifted by 16.5° to the east. Furthermore, a constant background X-ray radiation was added to our data similarly to Fig. 3 to better compare it to the observations.

Fig. 5 Simulated vertical profiles of the Mach number 𝓜, the density ρ, the pressure p, the absolute value of the flow speed in ɀ-direction υɀ, and the temperature T measured along a ray through the origin (injection region) of TDEs occurring at a rate of 10−5 yr−1 in the β-model halo at l6 Myr. The green and orange dotted lines indicate the position of the contact discontinuity and the forward shock, respectively.

3.3 Multicomponent Milky Way model

Since the lower TDE rate of 10⁻⁵ yr⁻¹ did not change the bubble structure significantly in the β-model but resembled the FBs slightly better than the higher TDE rate model due to fewer RT instabilities, the multicomponent model was studied with the lower TDE rate only. The simulation results are shown in the bottom-right panel of Fig. 2. The shape of the bubbles differs significantly from the previously presented models as the Milky Way density gradient does not decrease as steeply in the ɀ-direction in this model. The resulting bubbles are smaller but wider when they reach the size of the EBs. They also show no constriction in the mid-plane, thus not matching well with the EBs that are more spherical. The TDEs take 14.5 Myr to form the superbubbles in this model, totaling an energy injection of 1.4 × 10⁵⁶ erg, similar to the β-model.

Once again we produced a synthetic X-ray map, which is shown in the bottom panel of Fig. 3. It emphasizes that the simulated superbubbles in a multicomponent Milky Way model are too wide at the mid-plane, spanning over 4 h in longitude, while being as large as the EBs (red) in latitude. This can be seen once again in the brightness profiles in Fig. 4, especially at +50° and +60° latitude, where the brightness peaks are not that far separated or not separated at all, except at low latitudes as seen at +40°, where it is similar to the β-model. If the multicomponent model superbubbles grew larger, brightness peaks would also appear in the higher latitude data, potentially aligning better with the observations. Despite that, the peaks at +40° would separate further and thus might no longer match observations. Therefore, the previous β-model halo bubbles seem to reproduce the EBs better. Furthermore, the contact discontinuity in this model, which lies at the inner edge of the brighter region in Fig. 3, is too wide, but shows a similar size in latitude to the FBs (yellow). However, it does not reproduce the FBs much better than the β-model.

The shock propagation is very similar to the TDE shocks in the β-model halo, which is why we refrain from presenting the shock evolution once again. The strength of the shocks, however, is different from the previous models since the density in the interior of the bubbles is different. That is why we still show one-dimensional plots through the x = y = 0 line in Fig. 6 for this model, including the Mach number. Compared to the previous β-model plots, the quantities do not have strong gradients in the shell, as no RT instability formed in the contact discontinuity (Fig. 2). The bubbles also have a higher temperature and a lower density. The TDE shocks start with a very high Mach number of almost 𝓜 = 10 but also lose their energy more quickly – shocks are only found in the very center of the superbubbles. Therefore, the γ-rays would need to be accelerated more centrally by the initial high Mach TDE shock. The vertical speed of the forward shock is with ~500 km s⁻¹ also lower in this model. The interior of the bubbles again has a low density, a high temperature, and is almost isobaric, whereas the shell has a high density and a low temperature. The shell thickness was estimated from the density, the Mach number, and the flow speed to be 2 kpc, better matching the estimates of Predehl et al. (2020).

Concluding, the bubbles blown in the multicomponent model with our setup do not match any better with the EBs and FBs than the bubbles in the β-model do. Even though the shell is thicker and the FBs therefore smaller, the simulated FBs and EBs are too wide and the forward shock is not spherical.

On a different note, the asymmetric shape of the EBs remains without any explanation in our investigations since our set-up is axisymmetric and has a symmetry with respect to the Galactic plane. An asymmetric density distribution or a postulated circumgalactic wind (Mou et al. 2023) may be able to explain this property. Although Mou et al. (2023) use an AGN outflow as energy source, their considered outflow has similarities to the thermal energy outflow produced by the TDEs in our model since it has a constant 19 Myr-long low luminosity. The evolution of a TDE blown bubble with asymmetric features needs to be investigated in detail in future studies.

Fig. 6 Similar to Fig. 5 but for the multicomponent model.

4 Conclusions

Three-dimensional hydrodynamic simulations using a TDE model were carried out to reproduce the EBs and the FBs. Our findings can be summarized as follows:

1. Regular energy injections with a rate of 10⁻⁵ yr⁻¹ and an average luminosity of 3 × 10⁴¹ erg s⁻¹ at the Galactic center by for example TDEs can blow superbubbles into the Galactic halo, having a comparable spherical and smooth shape, size, and evolution time to the EBs. Simulations with a β-model halo reproduce the typical EB properties better than simulations with an exponential halo or multicomponent Milky Way model;

2. The general cocoon-like structure of the EBs around the FBs can be reproduced by our model. We predict that the γ-ray emission occurs at the contact discontinuity of the simulated bubbles;

3. Tidal disruption events create shock fronts inside the bubbles that may be able to continuously accelerate CRs to energies high enough to explain the FB γ-rays, be they leptonic or hadronic in origin. The shocks found in our study are not certain to be strong enough to do so but could be strengthened by more consecutive or stronger TDEs, and regular SN explosions in the disk, however, this needs to be investigated further;

4. Synthetic observations of our simulated superbubbles show that the X-ray emission occurs mostly in the dense shell. Together with the FBs at their contact discontinuity, they confirm the points 1 and 2 once again;

5. The TDE model can reproduce the brightness profiles observed with eROSITA fairly well. The brightness peaks at about ±40° longitude with a decrease in brightness between them are visible in both the observations and the synthetic maps. Furthermore, the absolute brightness of our generated observations is quite close to that of the actual observations;

6. Finally, we would like to stress that the agreement between the observations of the EBs and FBs, and our hydrodynamic models is fairly good, considering the tremendous impact of the detailed local gas density and pressure distributions in the Galactic halo on the shape of the bubbles, as is also well known from the modeling of large superbubbles in the ISM.

Nevertheless, some discrepancies between model and observations need further discussion. We note that especially the FBs are not matched perfectly by our model. We argue that this offset can be explained by the uncertainty of several input parameters that have a strong impact on the forward shock and contact discontinuity shape. Besides the Milky Way temperature and density distribution, for which we investigate several models but which are still quite unconstrained, the actual times, masses, and energies of past TDEs are unknown. Furthermore, CRs, radiative cooling, self-gravity, and magnetic fields may have an impact on the detailed and local structures including the RT instabilities, whereas the global structure should remain unchanged. Due to these uncertainties, more comprehensive simulations may still deviate from the detailed EB and FB shape. Only a thorough study of a large parameter space, including more Milky Way models, a non-constant TDE rate, a variable TDE energy, and short-lived TDE jets that would induce a direction to the outflow may give better insights particularly into the evolution of the FBs with a TDE model.

Concluding all findings, we showed that TDEs, which are expected to happen regularly in the Milky Way, or other regular energy injection events at the Galactic center, are a possible origin of the EBs and FBs. The fact that we demonstrated that it is possible to match the observations by our TDE model does not rule out that other mechanisms discussed in the literature contribute to the expansion and the emission. It is quite likely that jets were blown out during episodes of strong accretion onto the SMBH, although not in a continuous fashion. In addition, the star formation rate in the Galactic center is known to be higher than in the solar neighborhood, so that SNe, as already been mentioned, would also add their share.

Acknowledgements

The authors thank Chung-Ming Ko and Romain Teyssier for providing further information on their past paper and the efficiency of the RAMSES code, respectively. We further want to thank Bert Vander Meulen for the helpful discussions regarding our synthetic X-ray maps, Mattia Pacicco for test simulations of the impact of magnetic fields and useful discussions about the X-ray maps, as well as Nina Sartorio for their valuable suggestions and assistance in refining the manuscript. TS acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme DustOrigin (ERC-2019-StG-851622). Some of the results in this paper have been derived using the healpy and HEALPix package (Górski et al. 2005; Zonca et al. 2019). This work also made use of the Python packages NumPy (Harris et al. 2020), Matplotlib (Hunter 2007), MPI for Python (Dalcin et al. 2011), and the analysis and visualization software yt (Turk et al. 2011).

Appendix A Shocks inside the superbubble

To check whether the injected TDEs lead to shocks that may be able to accelerate CRs, we applied a shock tracer algorithm (see Sect. 3.2 for details) to the β-model halo bubble, as shown in Fig. A.1, and studied how these TDE shocks evolve over time.

Fig. A.1 Shock evolution between two TDEs in the β-model halo bubble. The shocked cells were found using a shock tracer algorithm and were marked by a gray dot. The right column shows a zoom-in of the inner 5 kpc of the left column figures. The elapsed time since the last TDE at 18 Myr is indicated. The next TDE happens shortly before snapshot (f).

References

All Figures

Fig. 1 Density distributions of the different Milky Way models dependent on the height ɀ and the radius R away from the Galactic center. Two-dimensional slices through the origins of the three-dimensional simulation boxes are shown. The exponential halo model (left), the β-model halo (middle), and the multicomponent model (right) are described by Eqs. (2), (3), and (4–7), respectively. The models differ by several orders of magnitude in their density distribution.

In the text

Fig. 2 Resulting density distributions of the numerical simulations for different Milky Way models and TDE rates. The snapshots are taken at the end of the simulations when the forward shocks reach the size of the EBs and show the detailed bubble structure, including the elapsed time since the first TDE. The different subfigures correspond to the results of the exponential halo model with energy and mass injections every 100–1000 yr (top left), the β-model halo with TDEs every 10 kyr (top right), the β-model halo with TDEs every 100 kyr (bottom left), and the multicomponent model with TDEs every 100 kyr (bottom right).

In the text

Fig. 3 Synthetic observations of the full three-dimensional simulation results using a Hammer–Aitoff projection in Galactic coordinates. The color-coded intensity is the surface brightness B in logarithmic units for the photons detected per time and solid angle and corresponds to detected X-ray photons with an energy between 0.6 keV and 1 keV. Since the cosmic X-ray background was not included in our model, we added it as a constant value to the maps so that we arrive at the same minimum brightness of 4.66 photons s−1 deg−2 as Predehl et al. (2020). The colorbar limits (obtained from Tseng et al. 2024) are the same as in the observations of Predehl et al. (2020). The main underlying model is the thermal collisional ionization equilibrium model of Cloudy. Both, the geometry for an observation from the Solar System and instrumental effects of the eROSITA telescope were considered for these figures. The FBs would appear at the contact discontinuity, which lies at the inner side of the brighter region (pink-white color) according to our model. The approximate shell of the observed EBs (red) and FBs (yellow) are indicated using data points obtained from Predehl et al. (2020, Fig. 3). We note that these points have a reading error of about ±5° and ±1 h because the edges of the observed EBs and FBs are not clearly defined everywhere. The top panel shows the β-model halo bubble (cf. Fig. 2 bottom left) and the bottom panel shows the multicomponent model halo bubble (cf. Fig. 2 bottom right).

In the text

Fig. 4 Brightness profiles of our synthetic observations at different latitudes. The top panel shows the surface brightness B for the β-model halo bubble and the bottom panel for the multicomponent model bubble for the latitudes ±40°, ±50°, and ±60° along the longitudes αG. Only the northern brightness profiles of our synthetic observations (black) are shown and compared to the data of the northern (red) and southern (blue) EB from Predehl et al. (2020). The northern eROSITA data was shifted by 16.5° to the east. Furthermore, a constant background X-ray radiation was added to our data similarly to Fig. 3 to better compare it to the observations.

In the text

Fig. 5 Simulated vertical profiles of the Mach number 𝓜, the density ρ, the pressure p, the absolute value of the flow speed in ɀ-direction υɀ, and the temperature T measured along a ray through the origin (injection region) of TDEs occurring at a rate of 10−5 yr−1 in the β-model halo at l6 Myr. The green and orange dotted lines indicate the position of the contact discontinuity and the forward shock, respectively.

In the text

	Fig. 6 Similar to Fig. 5 but for the multicomponent model.
In the text

	Fig. A.1 Shock evolution between two TDEs in the β-model halo bubble. The shocked cells were found using a shock tracer algorithm and were marked by a gray dot. The right column shows a zoom-in of the inner 5 kpc of the left column figures. The elapsed time since the last TDE at 18 Myr is indicated. The next TDE happens shortly before snapshot (f).
In the text