Microquasar Cygnus X-3 as the PeVatron Powering the Cygnus Bubble

Recent detections of dozens of ultrahigh-energy (UHE) γ-ray sources (Z. Cao et al. 2024) have revealed a growing population of Galactic PeVatrons—cosmic accelerators capable of boosting particles to PeV energies. Most of these sources are spatially extended γ-ray structures associated with pulsar wind nebulae (PWNe), molecular clouds in the proximity of middle-aged supernova remnants (SNRs), and bubbles surrounding microquasars, hereafter referred to as “microquasar halos” (μQSHs). Generally, the γ-ray emission does not spatially coincide with the particle accelerator itself; rather, the γ-ray morphology reflects the complex interplay between the spatial distribution of relativistic particles and the ambient target medium. Consequently, identifying the specific acceleration sites remains challenging and requires deep phenomenological and theoretical studies based on comprehensive multiwavelength analysis.

In μQSHs, relativistic particles may originate within the compact binary system, at the parsec-scale termination of the outflow (jets or winds), or both. For most sources, current observations cannot distinguish between these sites. Cygnus X-3, however, stands out as a unique case; the reported UHE γ-ray variability (The LHAASO Collaboration et al. 2025), specifically its modulation at the 4.8 hr orbital period, indicates that acceleration and emission occur within a compact region, likely associated with a (sub)relativistic wind or jet inside the binary. The production of UHE γ-rays extended beyond 1 PeV inside the binary system points to the presence of a hadronic super-PeVatron, precluding any “leptonic” contribution; because of severe synchrotron and inverse Compton losses, even at the maximum possible theoretical acceleration rate, electrons cannot reach PeV energies (The LHAASO Collaboration et al. 2025). Notably, a hadronic PeVatron localized within the binary system does not preclude the further acceleration of PeV particles in a large-scale “jet-termination” scenario.

Within the binary environment, the accelerated protons lose only a small fraction of their energy and escape into the surrounding medium without significant deformation of their initial (acceleration) spectrum. Upon entering the interstellar medium (ISM), these particles propagate diffusively, interacting with the ambient gas to form a vast γ-ray “halo” extending over hundreds of parsecs. The detectability of such a structure depends primarily on the medium’s diffusion coefficient and the target gas density. Generally, in the vicinity of the source, the cosmic-ray (CR) diffusion coefficient must be substantially smaller than the “standard” Galactic value to prevent the rapid runaway of protons and maintain an adequate surface brightness distribution. Such an enhancement of diffusive confinement is physically expected in the surroundings of powerful cosmic-ray accelerators (M. A. Malkov et al. 2013; M. D’Angelo et al. 2016; L. Nava et al. 2016, 2019; B. Schroer et al. 2022). Finally, the detection of these extended features requires instruments highly sensitive to low-surface-brightness emission; in this regard, the LHAASO KM2A detector has uniquely demonstrated its potential to reveal diffuse UHE γ-ray sources throughout the Galactic plane.

Notably, it is likely that a halo around Cygnus X-3 has already been detected. Indeed, the Cygnus Bubble, reported by the LHAASO collaboration as an extended UHE γ-ray source with an angular radius of ≈6 ° toward the Cygnus region (LHAASO Collaboration 2024), matches this description. This statement, at first glance, may seem paradoxical as it opposes the widely accepted view that the Cygnus Bubble is associated with the Cygnus X star-forming region, specifically the Cygnus OB2 association at d₁ = 1.4 kpc (F. Aharonian et al. 2019; S. Menchiari et al. 2024; L. Härer et al. 2025). One might argue that assuming the larger distance of Cygnus X-3 (d₂ ≈ 9.7 kpc) increases the required γ-ray luminosity by ${({d}_{2}/{d}_{1})}^{2}$, or roughly 1.5 orders of magnitude.

However, this does not present a fundamental problem. Depending on the radial dependence of the product of the CR and gas densities, a distant source can be energetically competitive. The observed γ-ray flux, F_γ, resulting from the interaction of CRs with the ambient gas, is proportional to

$$\begin{eqnarray}&&{F}_{\gamma }\propto \frac{1}{{d}^{2}}{\int }_{0}^{R}{N}_{{\rm{CR}}}(r){n}_{{\rm{g}}}(r)4\pi {r}^{2}\,dr.\end{eqnarray} \tag{ 1 }$$

Assuming power-law spatial distributions for the CR and gas densities, ${N}_{{\rm{CR}}}(r)={N}_{0}{r}^{-{\alpha }_{1}}$ and ${n}_{{\rm{g}}}(r)={n}_{0}{r}^{-{\alpha }_{2}}$, and noting that for a fixed angular size ϑ, the linear size of the source scales as R = ϑd, the flux scales as⁷

$$\begin{eqnarray}&&{F}_{\gamma }\propto {N}_{0}{n}_{0}{d}^{-({\alpha }_{1}+{\alpha }_{2})+1}.\end{eqnarray} \tag{ 2 }$$

Equation (1) assumes that the volume filled by CRs extends at least to the radius R. This imposes a lower limit on the propagation speed—and correspondingly on the diffusion coefficient—to ensure that particles could reach the physical distance R = ϑd within the operational lifetime of the accelerator.

Crucially, for a fixed accelerator power L_CR in a similar gas environment, the γ-ray flux does not follow the standard 1/d² dilution if the cosmic-ray distribution is maintained by continuous injection. From an observational perspective, given that a fixed angular size implies consistent instrumental performance, the association of the Cygnus Bubble with the distant Cygnus X-3 is physically consistent with the available LHAASO data (The LHAASO Collaboration et al. 2025). Below, we explore this possibility through rigorous calculations that strictly treat particle transport. Specifically, we model the diffusion of particles after their escape from the binary system into the ISM, utilizing a realistic gas distribution surrounding Cygnus X-3 on scales up to 1 kpc. This analysis accounts for the nonisotropic nature of the environment, considering propagation both along the Galactic disk and perpendicular to the Galactic plane.

The Letter is organized as follows: in Section 2, we describe the analytical and numerical methods used to calculate the CR distributions, the resulting γ-ray spectra, and the surface brightness profiles, based on the relevant assumptions regarding the gas distributions. In Section 3, we present the results for various injection parameters and discuss their physical implications. Finally, we summarize our findings, discuss the results, and provide concluding remarks in Section 4.

After being accelerated, particles (for simplicity, we consider only protons) escape the source and propagate diffusively through the surrounding ISM. Assuming spherical symmetry, with the accelerator located at the center and a homogeneous environment, the transport equation (F. A. Aharonian & A. M. Atoyan 1996; V. Bosch-Ramon et al. 2005) for CRs is

$$\begin{eqnarray}&&\frac{\partial N}{\partial t}-\frac{D}{{r}^{2}}\frac{\partial }{\partial r}\left({r}^{2}\frac{\partial N}{\partial r}\right)-\frac{\partial (bN)}{\partial E}=Q,\end{eqnarray} \tag{ 3 }$$

where N(E, r, t) is the differential number density of relativistic particles, D(E) is the spatial diffusion coefficient, b(E) = − dE/dt is the energy-loss rate, and Q(E, r, t) is the source function.

The dominant energy-loss process for CR protons is the inelastic collisions with ambient gas. The corresponding cooling time in a pure hydrogen medium is t_pp = ${(\kappa {\sigma }_{{\rm{inel}}}c{n}_{{\rm{H}}})}^{-1}\,\simeq 5.3\times 1{0}^{7}$ ${({n}_{{\rm{H}}}/1\,{{\rm{cm}}}^{-3})}^{-1}({\sigma }_{{\rm{inel}}}$/40 mb)⁻¹ yr, where c is the speed of light, n_H is the number density of ambient hydrogen atoms, κ ≃ 0.5 is the inelasticity coefficient, and σ_inel is the inelastic pp collision cross section.

The cross section σ_inel is a weak function of energy and, in the relativistic limit, can be presented, using the parameterization of E. Kafexhiu et al. (2014), as ${\sigma }_{{\rm{inel}}}(E)=30.7-0.96\mathrm{ln}x\,+0.18{(\mathrm{ln}x)}^{2}\,{\rm{mb}}$, with x = E/T_th, where T_th ≈ 0.2797 GeV, is the proton’s threshold kinetic energy for pion production.

The spatial diffusion coefficient is parameterized as

$$\begin{eqnarray}&&D(E)={D}_{0}{\left(\frac{E}{\,{\rm{1}}\,\,\,{\rm{PeV}}\,}\right)}^{\delta },\end{eqnarray} \tag{ 4 }$$

where D₀ and δ depend on the level and spectrum of magnetic turbulence in the ISM (A. W. Strong et al. 2007). Since the origin of interstellar turbulence is not fully understood, D₀ is treated as a free parameter. We consider both Kolmogorov (δ = 1/3) and Iroshnikov–Kraichnan (δ = 1/2) turbulence spectra. We do not take into account here that the diffusion coefficient transitions to scale as E² when the gyroradii of CRs are larger than the correlation length of interstellar magnetic turbulence (P. Subedi et al. 2017; G. Giacinti et al. 2018; O. Pezzi et al. 2022).

We assume continuous injection starting at t = 0, so that

$$\begin{eqnarray}&&Q(E,r,t)=q(E){\delta }^{3}({\boldsymbol{r}})\theta (t),\end{eqnarray} \tag{ 5 }$$

where δ³(r) is the Dirac delta function and θ(x) is the Heaviside step function. The injection spectrum is taken as

$$\begin{eqnarray}&&q(E)={Q}_{0}{E}^{-s}\exp \left(-\frac{E}{{E}_{0}}\right),\end{eqnarray} \tag{ 6 }$$

where s is the spectral index, and E₀ is the cutoff energy. The total kinetic power injected into CR protons is ${\dot{W}}_{{\rm{CR}}}={\int }_{{E}_{{\rm{\min }}}}^{\infty }q(E)EdE$; ${E}_{{\rm{\min }}}$ is the minimum proton energy.

We do not explicitly address the underlying acceleration mechanisms of CRs, but note that super-Eddington accreting X-ray binaries may accelerate protons up to several tens of PeV (e.g., E. Peretti et al. 2025; J. Wang et al. 2025). Throughout this Letter, we fixed the cutoff energy at E₀ = 10 PeV.

The solution of Equation (3) for the continuous injection is (S. I. Syrovatskii 1959)

$$\begin{eqnarray}\begin{array}{rcl}N(E,r,t) & = & {\displaystyle \int }_{0}^{t}\frac{b({E}^{{\prime} })}{b(E)}q({E}^{{\prime} })\frac{1}{{[4\pi \lambda (E,{E}^{{\prime} })]}^{3/2}}\\ & & \times \exp \left[-\frac{{r}^{2}}{4\lambda (E,{E}^{{\prime} })}\right]d{t}^{{\prime} },\end{array}\end{eqnarray} \tag{ 7 }$$

where ${E}^{{\prime} }$ is defined through

$$\begin{eqnarray}&&t-{t}^{{\prime} }={\int }_{E}^{{E}^{{\prime} }}\frac{d{E}_{1}}{b({E}_{1})},\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&\lambda (E,{E}^{{\prime} })={\int }_{E}^{{E}^{{\prime} }}\frac{D({E}_{1})}{b({E}_{1})}d{E}_{1}.\end{eqnarray} \tag{ 9 }$$

For ultrarelativistic protons, because of the weak energy dependence of σ_inel, we have ${E}^{{\prime} }=E\exp [(t-{t}^{{\prime} })/{t}_{{\rm{pp}}}]$, thus (F. A. Aharonian & A. M. Atoyan 1996)

$$\begin{eqnarray}&&\lambda (E,{E}^{{\prime} })=D(E)(t-{t}^{{\prime} })\frac{\exp [\delta (t-{t}^{{\prime} })/{t}_{{\rm{pp}}}]-1}{\delta (t-{t}^{{\prime} })/{t}_{{\rm{pp}}}}.\end{eqnarray} \tag{ 10 }$$

In the limit t ≪ t_pp, Equation (7) is reduced to

$$\begin{eqnarray}&&N(E,r,t)=\frac{q(E)}{D(E)}\frac{{\rm{erfc}}(r/{r}_{{\rm{d}}})}{4\pi r}.\end{eqnarray} \tag{ 11 }$$

Here, ${\rm{erfc}}(x)$ is the complementary error function and ${r}_{{\rm{d}}}=\sqrt{4D(E)t}$ is the diffusion radius.

2.1. Applying to Cygnus X-3

Due to the large distance of Cygnus X-3 and the strong Galactic extinction along its line of sight, a detailed, high-resolution study of the local ISM surrounding the source is currently lacking. To establish a baseline for the γ-ray production calculations, we first assume a spatially uniform ambient hydrogen density, n_H, and then proceed to a more realistic treatment that accounts for the Galactic disk’s global vertical structure.

Following empirical models of the atomic hydrogen (H i) distribution, we adopt a stratified density profile that accounts for the increase of the gas scale height with Galactocentric radius in the outer Galaxy (P. M. W. Kalberla & L. Dedes 2008; P. M. W. Kalberla & J. Kerp 2009). At the Galactocentric radius of Cygnus X-3 (R ≈ 10 kpc), the weakening of the Galactic gravitational potential leads to a vertically extended gas disk with a characteristic scale height of z_c = 300 pc (P. M. W. Kalberla & L. Dedes 2008).

We model this vertical dependence using a simple exponential profile

$$\begin{eqnarray}&&{n}_{{\rm{H}}}(z)={n}_{0}\exp \left(-\frac{| z-{z}_{{\rm{mid}}}| }{{z}_{c}}\right),\end{eqnarray} \tag{ 12 }$$

where n₀ is the midplane density and z_mid = 0.

This vertically extended gas distribution provides an extended hadronic target for CRs injected by the Cygnus X-3 jets. In this way, the model bridges the lack of local high-resolution ISM data with the well-established large-scale structure of the Milky Way, enabling a more realistic estimate of the diffuse UHE γ-ray luminosity produced via pp interactions in the Cygnus region.

Once the distribution of CRs around the accelerator is obtained, the emissivity per hydrogen atom (in units of photons s⁻¹ sr⁻¹ GeV⁻¹) of γ-rays generated from the inelastic pp collisions is computed as following (E. Kafexhiu et al. 2014)

$$\begin{eqnarray}&&{q}_{\gamma }({E}_{\gamma },r,t)=\frac{c}{4\pi }{\int }_{{E}_{\gamma }}^{\infty }\frac{d{\sigma }_{\gamma }}{d{E}_{\gamma }}({E}_{\gamma },E)N(E,r,t)dE,\end{eqnarray} \tag{ 13 }$$

where dσ_γ/dE_γ(E_γ, E) is the differential production cross section of γ-rays, which is provided by the package aafragpy (S. Koldobskiy et al. 2021). Given that the emissivity and gas distribution are known, the γ-ray intensity is

$$\begin{eqnarray}&&{I}_{\gamma }({E}_{\gamma },\vartheta ,\varphi ,t)={\varepsilon }_{{\rm{M}}}{\int }_{{\rm{L.O.S.}}}{n}_{{\rm{H}}}(z){q}_{\gamma }({E}_{\gamma },r,t)\,dl,\end{eqnarray} \tag{ 14 }$$

where the integration is along the line of sight (L.O.S.), whose direction is defined by the angular separation from the accelerator ϑ and the azimuthal angle φ, and ε_M = 2.0 is the nuclear enhancement factor accounting for the contribution to γ-ray production of nuclei heavier than hydrogen in both CRs and ambient gas (M. Mori 2009; M. Kachelriess et al. 2014). The nuclear enhancement is calculated assuming that the CR and ISM compositions are the same as the local ones. When the gas distribution is homogeneous, the γ-ray intensity is independent of φ. On the other hand, when the gas distribution is nonhomogeneous perpendicular to the Galactic plane, we can define the azimuthally averaged γ-ray intensity as ${(2\pi )}^{-1}{\int }_{0}^{2\pi }I({E}_{\gamma },\vartheta ,\varphi ,t)\,d\varphi$. The γ-ray flux is

$$\begin{eqnarray}&&\begin{array}{l}{F}_{\gamma }({E}_{\gamma },t)=\,{\displaystyle \int }_{0}^{\arcsin (R/d)}d\vartheta \,\sin \vartheta \cos \vartheta {\displaystyle \int }_{0}^{2\pi }d\varphi \,{I}_{\gamma }({E}_{\gamma },\vartheta ,\varphi ,t),\end{array}\end{eqnarray} \tag{ 15 }$$

where d is the distance to the accelerator, and R is the dimension of the emission region. We assume that R = 1 kpc, which corresponds to an angular radius of ∼6° for the emission region, given that the distance to Cygnus X-3 is d = 9.67 kpc (M. J. Reid & J. C. A. Miller-Jones 2023).

The CRs near their sources can excavate a cavity with a radius of ∼10–50 pc, resulting in partial evacuation of gas from the cavity (B. Schroer et al. 2022). We demonstrate here that this evacuation effect has no significant impact on our estimates of CR acceleration efficiency in Section 3. The diffusion radius of CRs is ${r}_{{\rm{d}}}=\sqrt{4D{t}_{{\rm{age}}}}\simeq 600\,{\rm{pc}}$ at 1 PeV with D = 3 ×10²⁹ cm² s¹ and t_age = 100 kyr (refer to Section 3 on the related discussion of model parameters). For the stationary injection, the CR radial distribution is proportional to 1/r when r ≪ r_d according to Equation (11). Assuming a homogeneous gas distribution, the γ-ray emissivity is also proportional to 1/r; then the γ-ray flux originating from within radius r is ${F}_{\gamma }\,(\lt r)\propto {\int }_{0}^{r}{{r}^{{\prime} }}^{-1}4\pi {{r}^{{\prime} }}^{2}d{r}^{{\prime} }=2\pi {r}^{2}$, and thus, ${F}_{\gamma }(\lt {r}_{1})/{F}_{\gamma }(\lt {r}_{2})\,={({r}_{1}/{r}_{2})}^{2}$. Now, taking r₁ = 50 pc and r₂ = 200 pc, we can see that F_γ( < 50 pc)/F_γ( < 200 pc) = 1/16, which means that the γ-ray flux originating from within 50 pc is much smaller than that from within 200 pc. Taking the evacuation into account, the contribution to γ-ray emission from within 50 pc is much smaller, and the γ-ray flux mainly originates from between 50 pc and r_d. As a consequence, our estimates on the efficiency are not affected significantly.

Before discussing the implications of our results, we briefly discuss the selection of model parameters. First, we assume that the injection spectral index s = 2.0 is fixed, which is consistent with the prediction of diffusive shock acceleration theory in strong shocks (e.g., see M. A. Malkov & L. O. Drury 2001). The strong termination shocks driven by fast outflows (winds and jets), which can reach mildly relativistic velocities, in ultraluminous X-ray sources such as the microquasar Cygnus X-3 (A. Veledina et al. 2024) could accelerate CRs effectively to several tens of PeV (e.g., see E. Peretti et al. 2025; J. Wang et al. 2025), while other underlying acceleration mechanisms can not be excluded. Second, since we do not know exactly how long the injection process sustains, we assume that the elapsed time denoted by t_age since the injection starts from t = 0 is a free parameter, and we choose four different values, i.e., t_age = 100, 200, 300, and 400 kyr, given that the companion of Cygnus X-3 is a Wolf–Rayet star (M. H. van Kerkwijk et al. 1992, 1996). Then, we leave only the diffusion coefficient normalization D₀, the ambient gas density n_H (n₀) for the homogeneous (nonhomogeneous) gas distribution, and the injection kinetic power ${\dot{W}}_{{\rm{CR}}}$ the other three free parameters. The γ-ray flux and intensity are proportional to the product ${n}_{{\rm{H}}}{\dot{W}}_{{\rm{CR}}}$, which can be determined by the observed data once we know D₀. As a fiducial value, we assume that n_H = 1.0 cm⁻³ (n₀ = 1.0 cm⁻³) for the homogeneous (nonhomogeneous) gas distribution. Once choosing D₀ = 3 × 10²⁹ cm² s⁻¹, we find our model can reasonably explain the observed data for both Iroshnikov–Kraichnan (δ = 1/2) and Kolmogorov (δ = 1/3) turbulence, by tuning ${\dot{W}}_{{\rm{CR}}}$.

The radial profile of γ-ray intensity can give crucial information on the spatial distribution of CRs and their injection history. To obtain this radial flux distribution, we extracted the photon distribution within the 6 ° radius region from Figure 1 of LHAASO Collaboration (2024), which contained a total of 66 photon-like events with energies exceeding 400 TeV, with an estimated CR background of 9.5. Using the energy spectrum from Figure 3 of LHAASO Collaboration (2024), we calculated the integrated energy above 400 TeV within the same 6° radius region. After subtracting the CR background, we converted the photon distribution into a radial flux distribution for energies above 400 TeV.

Meanwhile, we estimated the diffuse γ-ray emission utilizing the proton and helium spectra measured by DAMPE (Q. An et al. 2019; F. Alemanno et al. 2021) and LHAASO (Z. Cao et al. 2025, 2026). For the gas distribution, we employed the HI4PI survey data (HI4PI Collaboration et al. 2016) and the CfA ¹²CO data (T. M. Dame et al. 1987). The column density of neutral hydrogen was calculated by integrating over the entire velocity range using (T. L. Wilson et al. 2013)

$$\begin{eqnarray}&&{N}_{{\rm{HI}}}=1.8\times 1{0}^{18}\int {T}_{{\rm{mb}},{\rm{HI}}}\,dV.\end{eqnarray} \tag{ 16 }$$

The molecular hydrogen column density was derived by

$$\begin{eqnarray}&&{N}_{{{\rm{H}}}_{2}}=X\,{W}_{{}^{12}{\rm{CO}}}=X\int {T}_{{\rm{mb}}{,}^{12}{\rm{CO}}}\,dV.\end{eqnarray} \tag{ 17 }$$

Here, we adopted a mean CO-to-H₂ conversion factor of X = 2.0 × 10²⁰ cm⁻² K⁻¹ km s (A. D. Bolatto et al. 2013). The total hydrogen column density is ${N}_{{\rm{H}}}={N}_{{\rm{HI}}}+2{N{}_{{\rm{H}}}}_{2}$, and we assumed a helium abundance of N_He = 0.1N_H in the ISM. We then used aafragpy (S. Koldobskiy et al. 2021) to calculate the diffuse γ-ray emission expected from the local CR proton and helium spectra.

Figure 1 shows the fits of our model to the flux of γ-rays with energies larger than 400 TeV from the Cygnus Bubble within a radius of 6° as observed by LHAASO (LHAASO Collaboration 2024) and the radial profile of integrated γ-ray intensity above 400 TeV (i.e., ${\int }_{\,{\rm{400}}\,\,\,{\rm{TeV}}\,}^{\infty }{I}_{\gamma }(\vartheta ,{E}_{\gamma },t)d{E}_{\gamma }$) in the left and right panels, respectively, for a homogeneous gas distribution. In the left panel, the blue dashed–dotted–dotted line shows the diffuse Galactic γ-ray flux within a 6 ° radius region around Cygnus X-3, while the black and gray lines show the summation of our model and diffuse γ-ray fluxes for Iroshnikov–Kraichnan and Kolmogorov turbulence, respectively, when t_age = 100 (solid lines), 200 (dashed lines), 300 (dotted–dashed lines), and 400 (dotted lines) kyr. Moreover, in the right panel, the corresponding lines show the corresponding radial profiles of integrated γ-ray intensity above 400 TeV. As shown in Figure 1, our model can reasonably explain the observed γ-ray flux and intensity radial profile for energies above 400 TeV. In fact, given the parameters we chose as discussed above, the radial profile of integrated γ-ray intensity can be fitted reasonably without fine-tuning the parameters, once we fit our model to the observed γ-ray flux by tuning ${\dot{W}}_{{\rm{CR}}}$, whose best-fit values are listed in Table 1. Moreover, we find that the required injection kinetic power ${\dot{W}}_{{\rm{CR}}}$ decreases when t_age increases, as expected. However, when t_age > 400 kyr, we find ${\dot{W}}_{{\rm{CR}}}$ no longer decreases apparently, when increasing t_age further. Due to the interplay of continuous injection into and escaping of CRs from the emission region, their distribution tends to N(E, r, t) = q(E)/[4πrD(E)] when r ≪ r_d, which is satisfied as t_age is large enough. Given that the kinetic luminosity of Cygnus X-3 is estimated to be L_kin = 5 × 10³⁹erg s⁻¹ (A. Veledina et al. 2024; J. Wang et al. 2025), the required acceleration efficiency ($={\dot{W}}_{{\rm{CR}}}/{L}_{{\rm{kin}}}$) is 0.7%–1.6% according to Table 1. Even though we adopt a lower gas density, for instance, n_H = 0.1 cm⁻³, the required acceleration efficiency is 7%–16%, thus Cygnus X-3 has enough energy budget for accelerating CRs to about 10 PeV.

Zoom In Zoom Out Reset image size

Table 1. The CR Injection Kinetic Power for a Homogeneous Gas Distribution

tage	${\dot{W}}_{{\rm{CR}}}$ [1037 erg s−1]
(kyr)	Kraichnana	Kolmogorova
100	8.0	7.2
200	5.6	4.8
300	4.8	4.0
400	4.4	3.6

Note. ^aThe CR injection spectrum is an exponentially cutoff power-law function with the spectral index s = 2.0 and the cutoff energy E₀ = 10 PeV, while the minimum proton kinetic energy is ${E}_{{\rm{\min }}}=\,{\rm{1}}\,\,\,{\rm{GeV}}\,$. We have assumed that the gas distribution is homogeneous and its density is n_H = 1 cm⁻³. The diffusion coefficient normalization D₀ = 3 × 10²⁹ cm² s⁻¹ for both Iroshnikov–Kraichnan and Kolmogorov turbulence.

Download table as:  ASCIITypeset image

As we have discussed at the end of Section 2, the cavity excavated by the CR pressure (B. Schroer et al. 2022) has no significant impact on the γ-ray flux for an emission region with a size of 1000 pc. However, such a low-density cavity embedded in the emission region can leave an imprint on the γ-ray intensity profile. Here, we briefly discuss the aftermath due to the evacuation effect. We assume the cavity has a radius R_cavity = 100 pc, and the evacuated gas accumulates in a thin shell with a thickness ΔR = 5 pc located at the outer surface of the cavity. Therefore, the gas density in the shell is ${n}_{{\rm{shell}}}={n}_{{\rm{H}}}/[1-{(1-{\rm{\Delta }}R/{R}_{{\rm{cavity}}})}^{3}]=7.0{n}_{{\rm{H}}}$, where n_H = 1 cm⁻³ is the gas density outside of the cavity, which we assume is not affected by the evacuation effect. Figure 2 shows the integrated γ-ray intensity radial profile above 400 TeV taking the evacuation effect into account. All parameters are the same as those discussed in the previous paragraph, except the gas distribution. As shown in Figure 2, there is a characteristic peak around ϑ = 0$\mathop{.}\limits^{{\unicode{x000b0}}}$6 (=R_cavity/d) in the intensity profile, while we verify that the γ-ray flux is almost not affected, at least for energies above 1 TeV for each set of parameters (not shown). However, such a feature can not be resolved by the present LHAASO data, and the future high angular resolution observations may reveal if such a feature exists. Hereafter, we will not consider the evacuation effect.

Zoom In Zoom Out Reset image size

While a homogeneous gas distribution is not realistic and should only be regarded as an average over the emission region, we also consider a more realistic nonhomogeneous gas distribution, which has a finite scale height vertical to the Galactic plane as prescribed by Equation (12). We assume that the midplane gas density n₀ = 1.0 cm⁻³ and the scale height z_c = 300 pc. For such a gas distribution, the average gas density within the emission region is about 0.38 cm⁻³. In Figure 3, the left and right panels show the fits of our model to observed γ-ray flux and intensity radial profile, respectively, while the best-fit values of ${\dot{W}}_{{\rm{CR}}}$ are listed in Table 2, assuming D₀ = 3 × 10²⁹ cm² s⁻¹ for both Iroshnikov–Kraichnan and Kolmogorov turbulence. Similar to the case for a homogeneous gas distribution, our model for the nonhomogeneous gas distribution can also reasonably explain the observed γ-ray flux and intensity radial profile for energies above 400 TeV. According to Table 2, the required acceleration efficiency is 1.6%–3.2% for the nominal kinetic luminosity L_kin = 5 ×10³⁹ erg s⁻¹ of Cygnus X-3. Therefore, our results suggest that the γ-rays with energies above 400 TeV from Cygnus Bubble observed by LHAASO (LHAASO Collaboration 2024) may originate from the μQSH forming around the super-PeVatron microquasar Cygnus X-3, similar to the other five Galactic microquasars reported recently by LHAASO (LHAASO Collaboration et al. 2025).

Zoom In Zoom Out Reset image size

Table 2. The CR Injection Kinetic Power for a Nonhomogeneous Gas Distribution

tage	${\dot{W}}_{{\rm{CR}}}$ [1037 erg s−1]
(kyr)	Kraichnana	Kolmogorova
100	16.0	14.4
200	12.0	10.4
300	10.4	8.8
400	9.6	8.0

Note. ^aThe CR injection spectrum is an exponentially cutoff power-law function with the spectral index s = 2.0 and the cutoff energy E₀ = 10 PeV, while the minimum proton kinetic energy is ${E}_{{\rm{\min }}}=\,{\rm{1}}\,\,\,{\rm{GeV}}\,$. We have assumed that the gas distribution is nonhomogeneous and is given by Equation (12) with n₀ = 1 cm⁻³ and z_c = 300 pc. The diffusion coefficient normalization D₀ = 3 × 10²⁹ cm² s⁻¹ for both Iroshnikov–Kraichnan and Kolmogorov turbulence.

Download table as:  ASCIITypeset image

The spatial diffusion coefficient plays a vital role in determining the spatial distribution of CRs, and hence, the resulting γ-ray morphology. Therefore, we discuss it further. The spatial diffusion coefficient in the ISM is D(E =1 GeV) ∼ 3 × 10²⁸ cm² s⁻¹, based on the investigations on the propagation of CRs in the Galaxy and on the diffuse Galactic γ-ray emission (A. W. Strong et al. 2007). When extrapolating the empirical diffusion coefficient to E = 1 PeV, D(E = 1 PeV) ∼ 3 × 10³¹ cm² s⁻¹(3 × 10³⁰ cm² s⁻¹) for Iroshnikov–Kraichnan (Kolmogorov) turbulence, which is 2 (1) orders of magnitude larger than D₀ = 3 × 10²⁹ cm² s⁻¹ we obtained. Although self-generated turbulence from CR streaming instability can suppress diffusivity around their sources, this suppression is confined to a region of a few tens of parsecs (B. Schroer et al. 2022). On the other hand, strong extrinsic turbulence could be injected by the relativistic jets of Cygnus X-3, producing an extended region of suppressed diffusivity. The spatial diffusion coefficient of CRs according to quasi-linear theory (R. Schlickeiser 1989) for r_g < l_c is

$$\begin{eqnarray}&&\begin{array}{l}D(p)\simeq \frac{\beta c{r}_{{\rm{g}}}}{3}{\left(\frac{{l}_{{\rm{c}}}}{{r}_{{\rm{g}}}}\right)}^{1-\delta }\frac{1}{\eta }\\ \,=\,\left\{\begin{array}{ll}1.85\times 1{0}^{28}\,{{\rm{cm}}}^{2}\,{{\rm{s}}}^{-1}\beta {\left(\frac{pc}{{\rm{PeV}}}\right)}^{1/2}{\left(\frac{B}{\,{\rm{3}}\,\,\,{\rm{\unicode{x000B5}G}}\,}\right)}^{-1/2}{\left(\frac{{l}_{{\rm{c}}}}{{\rm{pc}}}\right)}^{1/2}{\eta }^{-1},\quad & \delta =1/2,\\ 2.19\times 1{0}^{28}\,{{\rm{cm}}}^{2}\,{{\rm{s}}}^{-1}\beta {\left(\frac{pc}{{\rm{PeV}}}\right)}^{1/3}{\left(\frac{B}{\,{\rm{3}}\,\,\,{\rm{\unicode{x000B5}G}}\,}\right)}^{-1/3}{\left(\frac{{l}_{{\rm{c}}}}{{\rm{pc}}}\right)}^{2/3}{\eta }^{-1},\quad & \delta =1/3,\end{array}\right.\end{array}\end{eqnarray} \tag{ 18 }$$

where c is the speed of light, β the ratio of CR velocity to c, r_g = pc/eB = 0.36 pc (pc/PeV)(B/3 µG)⁻¹ the gyroradius of CRs, l_c the correlation length of magnetic turbulence, and η = (δB/B)² the magnetic turbulence level. For r_g > l_c, D(p) ∝ βp² (P. Subedi et al. 2017; G. Giacinti et al. 2018; O. Pezzi et al. 2022). The typical correlation length of magnetic turbulence in the interarm regions is 10 pc ≲ l_c ≲ 100 pc (M. Haverkorn et al. 2008; J. F. Hollins et al. 2017). In the spiral arms, stellar sources dominate the energy injection for the turbulence cascade, and the typical correlation length is l_c ∼ 1 pc (M. Haverkorn et al. 2008). According to M. J. Reid & J. C. A. Miller-Jones (2023), Cygnus X-3 is located in the outer spiral arm. If Cygnus X-3 dominates the turbulence energy injection for the emission region with size R = 1000 pc, and if we assume l_c ≃ 1 pc, then according to Equation (18), the turbulence level should be η ∼ 1/10 for D₀ = 3 × 10²⁹ cm² s⁻¹, assuming the magnetic field strength B = 3 µG. In such a situation, however, Equation (18) is not applicable for p > 1 PeV/c. If we assume a larger correlation length l_c ≃ 10 pc, Equation (18) is applicable up to p = 10 PeV/c, and the required turbulence level is η ∼ 1/3. Though there are still many uncertainties in our understanding of the interaction between PeV CRs and magnetic turbulence (Y. Hu 2026), the transition of the scattering regime should be incorporated into the modeling of PeV CR propagation, which is not taken into account in the present work.

Finally, we briefly discuss the possibility of whether Cygnus X-3 can account for the γ-ray flux of the 6° radius Cygnus Bubble as observed by LHAASO, not limited to energies above 400 TeV. In order to account for the observed wideband γ-ray flux from 1 TeV to 2 PeV, a soft CR injection spectrum above 1 TeV is needed. Here, we assume an Iroshnikov–Kraichnan turbulence spectrum, i.e., δ = 1/2, and the diffusion coefficient normalization D₀ = 3 × 10²⁹ cm² s⁻¹ as obtained previously. Furthermore, we assume the CR injection duration t_age = 1 Myr. In order to fit to the observed flux, as shown in Figure 4, the CR injection spectral index we obtained is s = 2.45, assuming an exponential cutoff energy E₀ = 10 PeV. The required CR injection kinetic power above 1 TeV is ${\dot{W}}_{{\rm{CR}}}(\gt \,{\rm{1}}\,\,\,{\rm{TeV}}\,)=5.6\times 1{0}^{38}\,{\rm{erg}}\,{{\rm{s}}}^{-1}$, which is about 11% of the kinetic luminosity of Cygnus X-3. If we extrapolate the CR injection spectrum to 1 GeV, then the required CR injection kinetic power above 1 GeV is higher than the kinetic luminosity of Cygnus X-3. However, a harder CR injection spectrum below 1 TeV can not be excluded; thus, the LHAASO observation can not constrain the low-energy spectrum.

Zoom In Zoom Out Reset image size

In this work, we proposed a simple propagation model, which is based on the diffusion Equation (3), to explain the origin of γ-rays with energies above 400 TeV coming from the direction of Cygnus Bubble reported recently by LHAASO (LHAASO Collaboration 2024) within a region with a radius of 6°. With only a few free parameters, our model can reasonably explain the observed γ-ray flux and integrated intensity radial profile above 400 TeV within the 6° radius region, assuming that CRs are injected continuously by the microquasar Cygnus X-3 for a duration of 100 kyr with an energy spectral index s = 2.0 and an exponential cutoff energy E₀ = 10 PeV. The CR spatial diffusion coefficient at E = 1 PeV we obtained is D₀ = 3 × 10²⁹ cm² s⁻¹ for both Iroshnikov–Kraichnan and Kolmogorov turbulence, which is a plausible value in accordance with the quasilinear theory, though the transition of the scattering regime in the magnetic turbulence is not taken into account. Given that the kinetic luminosity of Cygnus X-3 is about L_kin = 5 × 10³⁹ erg s⁻¹ (A. Veledina et al. 2024; J. Wang et al. 2025), our results imply that the CR acceleration efficiency is 0.7%–1.6% (1.6%–3.2%) for the homogeneous (nonhomogeneous) gas distribution assuming a density n_H = 1.0 cm⁻³ (n₀ = 1.0 cm⁻³ with a scale height of 300 pc). Thus, our results suggest that Cygnus X-3, which is capable of accelerating CRs to 10 PeV, can explain UHE photons above 400 PeV coming from the direction of the Cygnus Bubble within a 6° radius region, despite its large distance (d = 9.67 kpc) from the Earth. This scenario presents a unique case where we simultaneously detect both the primary accelerator, Cygnus X-3, and the surrounding “CR halo” formed by the historical accumulation of particles injected by this source into the ISM.

While the UHE emission above 400 TeV is the main focus of this work, the lower-energy (GeV–TeV) emission of the Cygnus Bubble can be treated as foreground radiation from an extended γ-ray source linked to the Cygnus OB2 association. In other microquasars, such as SS 433, V4641 Sgr, and GRS 1915+105, the γ-ray observations have shown that acceleration happens further out in the jet termination shocks (A. U. Abeysekara et al. 2018; H. E. S. S. Collaboration et al. 2024; R. Alfaro et al. 2024; LHAASO Collaboration et al. 2025; A. Acharyya et al. 2026). Our results demonstrate that Cygnus X-3 is among these objects. In this case, Cygnus X-3 would be a “dual” source: the orbitally modulated PeV photons come from the very compact inner region, while the more extended jet termination regions can be responsible for multi-TeV CRs. Instead of just a stellar-wind cavity, the Cygnus Bubble at energies above 400 TeV can be seen as a massive microquasar nebula.

The future deployment of next-generation Imaging Atmospheric Cherenkov Telescopes (IACTs), such as the Cherenkov Telescope Array (CTA; Cherenkov Telescope Array Consortium et al. 2019), the Astrofisica con Specchi a Tecnologia Replicante Italiana (ASTRI; S. Vercellone et al. 2022), and the proposed Large Array of Cherenkov Telescopes (LACT; S. Zhang 2026), will be instrumental in validating the Cygnus X-3 injector hypothesis by providing superior angular resolution (<0$\mathop{.}\limits^{{\unicode{x000b0}}}$05) in the UHE band. While LHAASO has effectively identified the Cygnus region as a super-PeVatron, its current resolution remains insufficient to fully disentangle the complex line-of-sight superposition between the foreground Cygnus Cocoon and the background Cygnus X-3 environment. IACTs will allow for the spatial resolution of a compact, point-like core at the microquasar’s coordinates and the mapping of energy-dependent morphology, where a shrinking emission size at higher energies would serve as a classic signature of a discrete injector.

R.z.Y. is supported by the National Natural Science Foundation of China under grants 12393854 and 12588101, and by the natural science funding of Sichuan Province under grant 2025ZNSFSC0065. R.z.Y. gratefully acknowledges the support of Cyrus Chun Ying Tang Foundations and of the studio of Academician Zhao Zhengguo, Deep Space Exploration Laboratory. F.A. acknowledges the support from Science and Technology Department of Sichuan Province.

All authors contributed equally.