Fermi-LAT Research Articles

Fermi Blazars in the Zwicky Transient Facility Survey: Properties of Large Optical Variations

We analyze the optical light curve data, obtained with the Zwicky Transient Facility (ZTF) survey, for 47 γ-ray blazars monitored by the Large Area Telescope onboard the Fermi Gamma-ray Space Telescope (Fermi). These 47 sources are selected because they are among the Fermi blazars with the largest optical variations in the ZTF data. Two color–magnitude variation patterns are seen in them, with one being redder-to-stable-when-brighter (RSWB; in 31 sources) and the other being stable when brighter (in 16 sources). The patterns fit with the results recently reported in several similar studies with different data. Moreover, we find that the colors in the stable state of the sources share similar values, for which (after being corrected for the Galactic extinction) most sources are in a range of 0.4–0.55. This feature could be intrinsic and may be applied in, for example, study of the intragalactic medium. We also determine the turning points for the sources showing the RSWB pattern, after which the color changes saturate and become stable. We find a correlation between optical fluxes and γ-ray fluxes at the turning points. The physical implications of the correlation remain to be investigated, probably better with a sample of high-quality γ-ray flux measurements.

A Gamma-Ray Flare from TXS 1508+572: Characterizing the Jet of a z = 4.31 Blazar in the Early Universe

Blazars can be detected from very large distances due to their high luminosity. However, the detection of γ-ray emission of blazars beyond z = 3 has only been confirmed for a small number of sources. Such observations probe the growth of supermassive black holes close to the peak of star formation in the history of galaxy evolution. As a result from a continuous monitoring of a sample of 80 z > 3 blazars with the Fermi Large Area Telescope (Fermi-LAT), we present the first detection of a γ-ray flare from the z = 4.31 blazar TXS 1508+572. This source showed high γ-ray activity from 2022 February to August, reaching a peak luminosity comparable to the most luminous flares ever detected with Fermi-LAT. We conducted a multiwavelength observing campaign involving XMM-Newton, the Neil Gehrels Swift Observatory, the Effelsberg 100 m radio telescope, and the Very Long Baseline Array. In addition, we make use of the monitoring programs by the Zwicky Transient Facility and the Near-Earth Object Wide-field Infrared Survey Explorer at optical and infrared wavelengths, respectively. We find that the source is particularly variable in the infrared band on daily timescales. The spectral energy distribution collected during our campaign is well described by a one-zone leptonic model, with the γ-ray flare originating from an increase of external Compton emission as a result of a fresh injection of accelerated electrons.

Multiwavelength Investigation of γ-Ray Source MGRO J1908+06 Emission Using Fermi-LAT, VERITAS and HAWC

This paper investigates the origin of the γ-ray emission from MGRO J1908+06 in the GeV–TeV energy band. By analyzing the data collected by the Fermi Large Area Telescope, the Very Energetic Radiation Imaging Telescope Array System, and High Altitude Water Cherenkov, with the addition of spectral data previously reported by LHAASO, a multiwavelength study of the morphological and spectral features of MGRO J1908+06 provides insight into the origin of the γ-ray emission. The mechanism behind the bright TeV emission is studied by constraining the magnetic field strength, the source age, and the distance through detailed broadband modeling. Both spectral shape and energy-dependent morphology support the scenario that inverse Compton emission of an evolved pulsar wind nebula associated with PSR J1907+0602 is responsible for the MGRO J1908+06 γ-ray emission with a best-fit true age of T = 22 ± 9 kyr and a magnetic field of B = 5.4 ± 0.8 μG, assuming the distance to the pulsar d PSR = 3.2 kpc.

Diversity in Fermi/GBM Gamma-Ray Bursts: New Insights from Machine Learning

Classification of gamma-ray bursts (GRBs) has been a long-standing puzzle in high-energy astrophysics. Recent observations challenge the traditional short versus long viewpoint, where long GRBs are thought to originate from the collapse of massive stars and short GRBs from compact binary mergers. Machine learning (ML) algorithms have been instrumental in addressing this problem, revealing five distinct GRB groups within the Swift Burst Alert Telescope (BAT) light-curve data, two of which are associated with kilonovae (KNe). In this work, we extend our analysis to the Fermi Gamma-ray Burst Monitor catalog and identify five clusters using unsupervised ML techniques, consistent with the Swift/BAT results. These five clusters are well separated in the fluence-duration plane, hinting at a potential link between fluence, duration, and complexities (or structures) in the light curves of GRBs. Further, we confirm two distinct classes of KN-associated GRBs. The presence of GRB 170817A in one of the two KN-associated clusters lends evidence to the hypothesis that this class of GRBs could potentially be produced by binary neutron star mergers. The second KN-associated GRB cluster could potentially originate from neutron star–black hole mergers. Future multimessenger observations of compact binaries in gravitational waves and electromagnetic waves can be paramount in understanding these clusters better.

Joint Deconvolution of Astronomical Images in the Presence of Poisson Noise

We present a new framework for joint likelihood deconvolution (Jolideco) of a set of astronomical observations of the same sky region in the presence of Poisson noise. The observations may be obtained from different instruments with different resolution, and different point-spread functions (PSFs). Jolideco reconstructs a single flux image by optimizing the posterior distribution based on the joint Poisson likelihood of all observations under one of several prior distributions. Most notably, we employ a patch-based image prior that is parameterized via a Gaussian mixture model, which we train on high-signal-to-noise astronomical images, including data from the JWST and the GLEAM radio survey. This prior favors correlation structures among the reconstructed pixel intensities that are characteristic of those observed in the training images. It is, however, not informative for the mean or scale of the reconstruction. By applying the method to simulated data, we show that the combination of multiple observations and the patch-based prior leads to much improved reconstruction quality in many different source scenarios and signal-to-noise regimes. We demonstrate that with the patch prior Jolideco yields superior reconstruction quality relative to alternative standard methods such as the Richardson–Lucy method. We illustrate the results of Jolideco applied to example data from the Chandra X-ray Observatory and the Fermi Gamma-ray Space Telescope. By comparing the measured width of a counts-based and the corresponding Jolideco flux profile of an X-ray filament in SNR 1E 0102.2–7219, we find the deconvolved width of 0.″58 ± 0.″02 to be consistent with the theoretical expectation derived from the known width of the PSF.

New γ-Ray Sources Identified in All-sky Surveys Based on Fermipy’s Advanced Algorithm

We employ an efficient method for identifying γ-ray sources across the entire sky, leveraging advanced algorithms from Fermipy, and cleverly utilizing the Galactic diffuse background emission model to partition the entire sky into 72 regions, thereby greatly enhancing the efficiency of discovering new sources throughout the sky through multi-threaded parallel computing. After confirming the reliability of the new method, we applied it for the first time to analyze data from the Fermi Large Area Telescope (Fermi-LAT) encompassing approximately 15.41 yr of all-sky surveys. Through this analysis, we successfully identified 1379 new sources with significance levels exceeding 4σ, of which 497 sources exhibited higher significance levels exceeding 5σ. Subsequently, we performed a systematic analysis of the spatial extension, spectra, and light variation characteristics of these newly identified sources. We identified 21 extended sources and 23 sources exhibiting spectral curvature above 10 GeV. Additionally, we identified 44 variable sources above 1 GeV.

A Comprehensive Study on the Mid-Infrared Variability of Blazars

We present a comprehensive investigation of mid-infrared (MIR) flux variability at 3.4 μm (W1 band) for a large sample of 3816 blazars, using Wide-field Infrared Survey Explorer (WISE) data through December 2022. The sample consists of 1740 flat-spectrum radio quasars (FSRQs), 1281 BL Lac objects (BL Lacs), and 795 blazars of uncertain type (BCUs). Considering Fermi Large Area Telescope detection, we classify 2331 as Fermi blazars and 1485 as non-Fermi blazars. Additionally, based on synchrotron peak frequency, the sample includes 2264 low-synchrotron peaked (LSP), 512 intermediate-synchrotron peaked (ISP), and 655 high-synchrotron peaked (HSP) sources. We conduct a comparative analysis of short- and long-term intrinsic variability amplitude (σm), duty cycle (DC), and ensemble structure function (ESF) across blazar subclasses. The median short-term σm values were 0.181−0.106+0.153, 0.104−0.054+0.101, 0.135−0.076+0.154, 0.173−0.097+0.158, 0.177−0.100+0.156, 0.096−0.050+0.109, and 0.106−0.058+0.100 mag for FSRQs, BL Lacs, Fermi blazars, non-Fermi blazars, LSPs, ISPs, and HSPs, respectively. The median DC values were 71.03−22.48+14.17, 64.02−22.86+16.97, 68.96−25.52+15.66, 69.40−22.17+14.42, 71.24−21.36+14.25, 63.03−33.19+16.93, and 64.63−24.26+15.88 percent for the same subclasses. The median long-term σm values were 0.137−0.105+0.408, 0.171−0.132+0.206, 0.282−0.184+0.332, 0.071−0.062+0.143, 0.218−0.174+0.386, 0.173−0.132+0.208, and 0.101−0.077+0.161 mag for the same subclasses, respectively. Our results reveal significant differences in 3.4 μm flux variability among these subclasses. FSRQs (LSPs) exhibit larger σm and DC values compared to BL Lacs (ISPs and HSPs). Fermi blazars display higher long-term σm but lower short-term σm relative to non-Fermi blazars, while DC distributions between the two groups are similar. ESF analysis further confirms the greater variability of FSRQs, LSPs, and Fermi blazars across a wide range of time scales compared to BL Lacs, ISPs/HSPs, and non-Fermi blazars. These findings highlight a close correlation between MIR variability and blazar properties, providing valuable insights into the underlying physical mechanisms responsible for their emission.

A Comprehensive Timing Analysis of Individual Pulses in X-Ray Bursts from SGR J0501+4516

The pulses in X-ray-burst (XRB) light curves from soft gamma-ray repeaters (SGRs) are generally thought to arise from magnetic crustal fractures or magnetic reconnection, reflecting the evolution of the energy release process in magnetars. In this study, we conduct a comprehensive timing analysis of 27 XRBs from SGR J0501+4156 detected by the Gamma-ray Burst Monitor on board Fermi. Utilizing a improved pulse-finding algorithm, we identify a total of 95 pulses and fit them using multiple FRED functions to obtain pulse-shape parameters based on the Markov Chain Monte Carlo method. We calculate the minimum variability timescales (MVTs) of the XRBs based on the shortest pulse; the distribution of MVTs follows a log-Gaussian function with a mean of 7.22−2.11+2.93 ms (1σ). The distributions of rise time, decay time, waiting time, width, skewness, and peakedness all follow the log-Gaussian function, and multiple power-law dependencies are observed between them; for example, a power-law positive correlation between decay time and rise time with 4.7σ, and a power-law negative correlation between pulse width and peakedness with 6.8σ. Besides, there is a positive correlation with 3.7σ between the number of pulses and burst duration. Our findings favor a magnetospheric origin, and some similarities with gamma-ray bursts imply that they have similar radiation mechanisms, e.g., magnetic reconnection processes.

On the Correlation between Young Massive Star Clusters and Gamma-Ray Unassociated Sources

Star clusters (SCs) are potential cosmic-ray accelerators and therefore are expected to emit high-energy radiation. However, a clear detection of gamma-ray emission from this source class has only been possible for a handful of cases. This could in principle result from two different reasons: either detectable SCs are limited to a small fraction of the total number of Galactic SCs, or gamma-ray-emitting SCs are not recognized as such and therefore are listed in the ensemble of unidentified sources. In this Letter we investigate this latter scenario by comparing available catalogs of SCs and H ii regions, obtained from Gaia and Wide-field Infrared Survey Explorer observations, to the gamma-ray GeV and TeV catalogs built from Fermi Large Area Telescope (LAT), H.E.S.S., and LHAASO data. The significance of the correlation between catalogs is evaluated by comparing the results with simulations of synthetic populations. A strong correlation emerges between Fermi-LAT-unidentified sources and H ii regions that trace massive SCs in the earliest (≲1–2 Myr) phase of their life, where no supernova explosions have happened yet, confirming that winds of massive stars can alone accelerate particles and produce gamma-ray emission at least up to GeV energies. The association with TeV energy sources is less evident. Similarly, no significant association is found between Gaia SCs and GeV nor TeV sources. We ascribe this fact to the larger extension of these objects but also to an intrinsic bias in the Gaia selection toward SCs surrounded by a lower target gas density, which would otherwise hinder the detection in the optical wave band.

Evidence for hybrid gamma-ray emission from the supernova remnant G150.3+4.5

The supernova remnant (SNR) G150.3+4.5 was first identified in radio, exhibiting a hard GeV spectrum and a ~1.5º radius. Radio observations revealed a bright arc with an index of ~−0.40, which stands in contrast to the index of ~−0.69 for the rest. This arc is coincident with the point-like Fermi source 4FGL J0426.5+5434 and KM2A source 1LHAASO J0428+5531. The rest of the SNR has a hard GeV spectrum and a soft TeV spectrum, implying a spectral cut-off or break near 1 TeV. Since there is no X-ray counterpart and no pulse signal detected, the gamma-ray (γ-ray) emission mechanism from the SNR and the point-like source appear puzzling. In this work, we reanalyse the γ-ray emission using 14 yr data recorded by Fermi Large Area Telescope and find that the spectrum of the northern half-sphere is compatible with a broken power law with a break at 146 ± 11 eV and photon indices of ΓNorthlobe = 1.54 ± 0.04stat ± 0.07syst (2.28 ± 0.08stat ± 0.12syst) below (above) the break. In addition, the southern half-sphere can be described well with a single power law with ΓSouthlobe =1.95 ± 0.07stat ± 0.09syst. Since the southern half-sphere is well correlated with CO emission, we propose that the γ-ray emission of the northern half-sphere could be dominated by relativistic electrons via inverse-Compton processes, while the southern half-sphere is dominated by cosmic rays via hadronic processes. 4FGL J0426.5+5434 may result from the illumination of a cloud by escaping cosmic rays or recent shock-cloud interaction. Observations from LHAASO-KM2A thus favour the possibility of a cosmic-ray PeVatron candidate, however, leptonic scenarios cannot be ruled out. Further multi-wavelength observations are warranted to confirm the hadronic nature of 1LHAASO J4028+5531.

Revisit of GeV Gamma-Ray Emission from Orion B with the Fermi Large Area Telescope

We revisit the γ-ray emission above 300 MeV towards the massive star-forming region of Orion B by adopting 14 yr observations with the Fermi Large Area Telescope and utilizing the updated software tools. The extended γ-ray emission region around Orion B is resolved into two components (region I and region II). The γ-ray spectrum of region I agrees with the predicted γ-ray spectrum assuming the cosmic ray (CR) density is the same as that of Alpha Magnetic Spectrometer (AMS-02) measured locally. The γ-ray emissivity of region II appears to be deficit at low energy band (E < 3 GeV). Through modeling we find that CR densities exhibit a significant deficit below 20 GeV, which may be caused by a slow diffusion inside the dense region. This is probably caused by an increased magnetic field whose strength increases with the gas density.

GRB 231129C: Another Thermal Emission Dominated Gamma-Ray Burst

This study presents detailed time-integrated and time-resolved spectral analysis of the Fermi Gamma-ray Burst Monitor observations of the bright GRB 231129C. The results reveal its distinct spectral characteristics, featuring a hard low-energy spectral index (α) and soft high-energy spectral index (β), similar to GRB 090902B, suggesting a possible dominance of thermal emission. Further analysis indicates that 92% of the spectral indices exceed the synchrotron “line of death,” with the hardest index at α ∼ +0.44. Simultaneously, 53% of the spectra can be well fitted by the nondissipative photosphere model, supporting a potential origin from a nondissipative photosphere. Additionally, we observe strong correlations between the spectral index α and peak energy E p with flux. For the α−F relationship, we employ F = F 0 e (3.00±0.10)α to describe it, whereas the E p−F relationship requires a smoothly bending power-law function. Based on the framework proposed by Hascoët et al. and Gao & Zhang, the jet characteristics of this burst were studied, revealing that both methods support the suitability of a pure fireball model for this GRB at small initial jet radii.

Likely Detection of GeV γ-Ray Emission from Pulsar Wind Nebula G32.64+0.53 with Fermi-LAT

Abstract In this study, we report the likely GeV γ-ray emissions originating from the pulsar PSR J1849-0001's pulsar wind nebula (PWN) G32.64+0.53. Our analysis covers approximately 14.7 yr of data from the Fermi Large Area Telescope Pass 8. The position of the source and its spectrum matches those in X-ray and TeV energy bands, so we propose that the GeV γ-ray source is indicative of PWN G32.64+0.53. We interpret the broadband spectral energy distribution (SED) using a time-dependent one-zone model, which assumes that the multiband nonthermal emission of the target source can be generated by synchrotron radiation and inverse Compton scattering (ICS) of the electrons/positrons. Our findings demonstrate that the model substantially elucidates the observed SED. These results lend support to the hypothesis that the γ-ray source originates from the PWN G32.64+0.53 powered by PSR J1849-0001. Furthermore, the γ-rays in TeV bands are likely generated by electrons/positrons within the nebula through ICS.

Two Models for the Orbital Modulation of Gamma Rays in Cyg X-3

We model the currently available γ-ray data on Cyg X-3 from the Fermi Large Area Telescope. Thanks to Cyg X-3’s very strong γ-ray activity during 2018–2021, the data quality has significantly improved. We study the strong orbital modulation of the γ-rays observed at high γ-ray fluxes. The modulation, as found earlier, is well modeled by anisotropic Compton scattering of the donor blackbody emission by relativistic electrons in a jet strongly misaligned with respect to the orbital axis. We confirm that this model fits well both the average γ-ray modulation light curve and the spectrum. However, we find that if the jet were aligned with the spin axis of a rotating black hole, it would undergo geodetic precession with a period of ∼50 yr. However, its presence is ruled out by both the γ-ray and radio data. Therefore, we consider an alternative model in which the average jet direction is aligned, but it is bent outside the orbit owing to the thrust of the donor stellar wind, and thus precesses at the orbital period. The γ-ray modulation then appears as a result of the variable Doppler boosting of synchrotron self-Compton jet emission. This model also fits the data well. However, the fitted bending angle is much larger than the theoretical one based on the binary and wind parameters as currently known. Thus, both models disagree with important aspects of our current theoretical understanding of the system. We discuss possible ways to find the correct model.

A New Puzzling Periodic Signal in GeV Energies of the γ-Ray Binary LS I+61°303

LS I+61°303 is a high-mass X-ray binary system comprising a massive Be star and a rapidly rotating neutron star. Its spectral energy distribution across multiwavelengths categorizes it as a γ-ray binary system. In our analysis of LS I+61°303 using Fermi Large Area Telescope observations, we not only confirmed the three previously discussed periodicities of orbital, superorbital, and orbital–superorbital beat periods observed in multiwavelength observations, but also identified an additional periodic signal. This newly discovered signal exhibits a period of ∼26.3 days at a ∼7σ confidence level. Moreover, the power spectrum peak of the new signal gradually decreases as the energy increases across the energy ranges of 0.1–0.3, 0.3–1.0, and 1.0–500.0 GeV. Interestingly, a potential signal with a similar period was found in data obtained from the Owens Valley Radio Observatory 40 m telescope. We suggest that the newly discovered periodic signal may originate from a coupling between the orbital period and the retrograde stellar precession period.

Energy-dependent intrinsic time delay of gamma-ray bursts on testing Lorentz invariance violation

High-energy photons of gamma-ray bursts (GRBs) might be emitted at different intrinsic times with energy dependence at the source. In this letter, we expand the model from previous works on testing the Lorentz Invariance Violation (LV) with the observed GRB data from the Fermi Gamma-ray Space Telescope. We reanalyze the previous data with the full Bayesian parameter estimation method and get consistent results by assuming that the time delays are due to an LV term and a constant intrinsic time delay term. Subsequently, we neglect the LV effect and only consider the intrinsic time delay effect. We assume a common intrinsic time delay term along with a source energy correlated time delay of high-energy photons. We find that the energy-dependent emission times can also explain the observed GRB data of high-energy photon events. Finally, we integrate these two physical mechanisms into a unified model to distinguish and evaluate their respective contributions using the observed GRB data.

CTA and SWGO can discover Higgsino dark matter annihilation

Thermal Higgsino dark matter (DM), with a mass near 1.1 TeV, is one of the most well-motivated and untested DM candidates. Leveraging recent hydrodynamic cosmological simulations that give DM density profiles in Milky Way analog galaxies we show that the linelike gamma-ray signal predicted from Higgsino annihilation in the Galactic Center could be detected at high significance with the upcoming Cherenkov Telescope Array (CTA) and Southern Wide-field Gamma-ray Observatory (SWGO) for all but the most pessimistic DM profiles. We perform the most sensitive search to-date for the linelike signal using 15 years of data from the Fermi Large Area Telescope, coming within an order one factor of the necessary sensitivity to detect the Higgsino for some Milky Way analog DM density profiles. We show that H.E.S.S. has subleading sensitivity relative to Fermi for the Higgsino at present. In contrast, we analyze H.E.S.S. inner Galaxy data for the thermal wino model with a mass near 2.8 TeV; we find no evidence for a DM signal and exclude the wino by over a factor of two in cross section for all DM profiles considered. In the process, we identify and attempt to correct what appears to be an inconsistency in previous H.E.S.S. inner Galaxy analyses for DM annihilation related to the analysis effective area, which may weaken the DM cross-section sensitivity claimed in those works by around an order of magnitude. Published by the American Physical Society 2024

Searching for Rapid Pulsations in Solar Flare X-Ray Data

Most studies of quasiperiodic pulsations (QPPs) in solar flares have identified characteristic periods in the 5–300 s range. Due to observational limitations, there have been few attempts to probe the <5 s period regime and understand the prevalence of such short-period QPPs. However, the Fermi Gamma-ray Burst Monitor (GBM) has observed approximately 1500 solar flares to date in high-cadence 16 Hz burst mode, providing us with an opportunity to study short-period QPPs at X-ray energies. We systematically analyze every solar flare observed by Fermi/GBM in burst mode, estimating the prevalence of QPPs in multiple X-ray energy bands. To better understand these results, we complement this with an analysis of synthetic solar flare lightcurves, both with and without oscillatory signals present. Using these synthetic lightcurves, we can understand the likely false-alarm and true-positive rates in the real solar GBM data. We do not find strong evidence for widespread short-period QPPs, indicating either a low base occurrence rate of such signatures or that their typical signal-to-noise ratios must be low—less than 1—in Fermi/GBM data. Finally, we present a selection of the most interesting potential QPP events that were identified in the GBM solar X-ray data.

Detection of Increasing γ-Ray Emissions from 4FGL J1718.5+4237 with Fermi-LAT

We report a gradual brightening γ-ray source, 4FGL J1718.5+4237, in 0.1–500.0 GeV, which may be associated with a blazar NVSS J171822+423948 with a redshift of ∼2.7. We analyzed 15.25 yr of γ-ray data recorded by the Large Area Telescope on board the Fermi Gamma-ray Space Telescope and detected significant γ-ray emissions in the direction of the blazar with a test statistic (TS) of ∼135. Based on timing analysis using a 1 yr time bin, we have observed a gradual brightening in γ-ray emissions from the target. In our analysis, we categorize them into two states: Quiet (TS ∼ 0) and Loud (TS ∼ 226) states, with the distinction occurring in 2016 August (MJD 57602.69). From the Quiet state to the brightest period (the last data point), the γ-ray flux in 0.1–500.0 GeV increased by more than 12-fold from <0.2 × 10−8 photons cm−1 s−1 to 2.6 × 10−8 photons cm−1 s−1. Additionally, we studied the spectral properties in detail for the Loud state and the overall data. While no significant variation was detected, both exhibited a spectral index Γ of ∼2.6. The origin of the brightening γ-ray emissions from the target is not yet clear. Future long-term multi-wavelength observations and studies will provide insight into the astrophysical mechanisms of the target.

Can the Symmetric Fermi and eROSITA Bubbles Be Produced by Tilted Jets?

The Fermi Gamma-Ray Space Telescope reveals two large bubbles in the Galaxy, extending nearly symmetrically ∼50° above and below the Galactic center (GC). Previous simulations of bubble formation invoking active galactic nucleus (AGN) jets have assumed that the jets are vertical to the Galactic disk; however, in general, the jet orientation does not necessarily correlate with the rotational axis of the Galactic disk. Using three-dimensional special relativistic hydrodynamic simulations including cosmic rays (CRs) and thermal gas, we show that the dense clumpy gas within the Galactic disk disrupts jet collimation (“failed jets” hereafter), which causes the failed jets to form hot bubbles. Subsequent buoyancy in the stratified atmosphere renders them vertical to form the symmetric Fermi and eROSITA bubbles (collectively, Galactic bubbles). We find that (1) despite the relativistic jets emanating from the GC at various angles ≤45° with respect to the rotational axis of the Galaxy, the Galactic bubbles nonetheless appear aligned with the axis; (2) the edge of the eROSITA bubbles corresponds to a forward shock driven by the hot bubbles; (3) followed by the forward shock is a tangling contact discontinuity corresponding to the edge of the Fermi bubbles; (4) assuming a leptonic model we find that the observed gamma-ray bubbles and microwave haze can be reproduced with a best-fit CR power-law spectral index of 2.4; The agreements between the simulated and the observed multiwavelength features suggest that forming the Galactic bubbles by oblique AGN failed jets is a plausible scenario.