© The Authors 2025

1. Introduction

Ultra-luminous X-ray sources (ULXs; see review King et al. 2023) are loosely defined as point-like sources located away from the galactic nucleus and characterized by fluxes, which, assuming isotropic emission, translate to an X-ray luminosity L_X≳10³⁹ erg s⁻¹. Over the past decade our understanding of ULXs has vastly expanded, with almost 2000 ULXs having been identified thus far (Walton et al. 2022). Moreover, after the discovery of pulsating ULXs (e.g., Bachetti et al. 2014; Israel et al. 2017b; Amato et al. 2023) and various spectral studies (e.g., Pintore et al. 2014; Koliopanos et al. 2017), it is now widely accepted that ULX engines are fueled by super-Eddington accretion onto stellar-mass compact objects (Israel et al. 2017a), instead of sub-Eddington accretion onto intermediate-mass black holes as early studies suggested (Colbert & Mushotzky 1999).

Extreme accretion goes hand in hand with strong winds and outflows (Fabrika et al. 2015; Kaaret et al. 2017; Weng & Feng 2018). Outflows have been seen in Galactic X-ray binaries hosting both neutron stars and black holes accreting below the Eddington limit (e.g., Ponti et al. 2012; Vincentelli et al. 2023). However, perhaps one of the early arguments for super-Eddington accretion came from the studies of a low-mass X-ray binary, Cyg X-2. King & Ritter (1999) argued that the system had survived a phase of extreme accretion (10⁻⁶−10⁻⁵ M_⊙ yr⁻¹) during which about only 0.1 M_⊙ was accreted and 2−3 M_⊙ were expelled, suggesting that Cyg X-2 had once gone through a ULX phase. For ULXs, observing campaigns with the Reflection Grating Spectrometer of XMM-Newton revealed the presence of blueshifted absorption features consistent with outflow velocities of 0.2−0.3 c (e.g., Pinto et al. 2017, 2020). It should also be noted that the presence of these winds and their properties are also supported by simulations (Kitaki et al. 2021). These radiative-driven outflows are optically thick and can reprocess the accretion luminosity, while due to their inhomogeneity they imprint variability onto the X-ray emission and spectral hardness (Middleton et al. 2015; Gúrpide et al. 2021). The basic idea is that these optically thick outflows block escaping radiation along certain lines of sight, allowing a distant observer to only view the central engine along a narrow funnel. This geometric structure would cause collimation of the radiation – commonly referred to as beaming – along the funnel and block radiation for an edge-on observer. Since beaming enhancement would be proportional to mass-accretion rate (King 2009), the majority of ULXs could perhaps be hidden from us due to their unfavorable orientation.

One of those hidden ULXs is the Galactic source SS 433, where studies have argued that the central engine accretes at a rate of 10⁻⁵−10⁻⁴ M_⊙ yr⁻¹ (Shkovskii 1981). The properties of the outflows in SS 433 were recently studied by Middleton et al. (2021). The authors modeled broadband X-ray spectra from NuSTAR and demonstrated that emission is consistent with reflection from the outflow funnel walls with an outflow velocity of 0.15−0.3 c and a wind-cone half-opening angle of 10 degrees.

Cygnus X-3 was recently characterized as another hidden ULX by IXPE observations (Veledina et al. 2024a). Specifically, the authors found that the energy-independent linear polarization is orthogonal to the direction of the radio ejections, indicating the presence of a collimating outflow with a half-opening angle less than 15 degrees. Yet, subsequent studies have also suggested that a face-on observer would observe Cygnus X-3 with an apparent luminosity of 5·10³⁹erg s⁻¹, making the system a bona fide ULX (Veledina et al. 2024b).

Quasi-spherical, shocked ionized structures are typically observed in association with ULXs (see, e.g., Pakull & Mirioni 2002; Soria et al. 2021; Gúrpide et al. 2022). These structures can extend for more than 10² pc in diameter, implying a million year lifetime of the associated engine. These observations suggest that ULX-driven winds not only result in relatively collimated jets, but typically they expand with a wide opening angle forming a bubble structure (Belfiore et al. 2020; Zhou et al. 2022, 2023). Wind bubbles are extremely interesting objects in the context of high-energy astrophysics as they can provide an ideal environment for diffusive shock acceleration (DSA). In fact, the evolution of wind bubbles results in the formation of strong shocks (Weaver et al. 1977; Koo & McKee 1992), where cosmic rays (CRs) can be accelerated. In the context of DSA, the high velocity of ULX winds is of particular interest as the acceleration timescale of CRs is expected to scale as the inverse of the second power of the wind speed (Drury 1983). This would make ULX winds extremely efficient accelerators up to the PeV range, which is believed to be the range of maximum energy for Galactic CR sources (Blasi 2013; Gabici et al. 2019). The nature of the Galactic CR accelerators at PeV energies, often referred to as PeVatrons (Cristofari 2021), is one of the main open questions in the CR community that has not been answered yet. However, recent gamma-ray observations performed by LHAASO (Cao et al. 2024) identified several Galactic sources whose non-thermal radiation clearly implies the presence of PeV particles. Some of these sources coincide with Galactic super-Eddington accreting objects (LHAASO Collaboration 2024), which could be “hidden” ULXs. This supports the hypothesis that super-Eddington sources are a new class of PeVatrons.

Fast winds from stellar-mass compact objects and microquasars were known as sites of non-thermal emission (see, e.g., Romero & Orellana 2005; Escobar et al. 2021, 2022; Abaroa et al. 2023). In addition, Abaroa et al. (2024), studying the timescales of relevant high-energy processes, already identified microquasar jets powered by super-Eddington accretion as a promising site for proton acceleration up to the PeV range (see also Sridhar et al. 2024 for similar implications in super-Eddington hypernebulae).

In this work, we developed a model of DSA taking place in the wind bubbles excavated by ULXs. We did so by solving the energy-dependent and space-dependent transport equation and analyzing the possible multi-messenger implications in terms of CRs, gamma rays, radio, and high-energy neutrinos. Our investigation shows that ULX wind bubbles are extremely promising cosmic accelerators that in several conditions can energize CR protons up to the PeV range. Therefore, we propose Galactic ULXs as a novel class of PeVatron candidates through DSA.

The manuscript is organized as follows. In Sect. 2, we outline the ULX wind-bubble structure and the acceleration and transport model for non-thermal particles. In Sect. 3, we present our results in terms of accelerated particle spectra and maximum energies achieved at ULX wind shocks. In Sect. 4, we describe how we applied our model to SS 433, a hidden Galactic ULX, and discuss our predictions in light of recent gamma-ray observations at very high energies. We also analyzed the possible neutrino and radio counterparts of our ULX bubble model. In Sect. 5, we discuss the contribution of Galactic ULXs to the observed cosmic-ray flux on Earth at energies close to the knee of the spectrum. Finally, we conclude in Sect. 6.

2. Model

When a compact object accretes mass from a nearby companion, the accretion flow is typically characterized by a certain angular momentum. This results in the formation of an accretion disk. Due to magnetic and viscous forces between disk annuli, there is an outwards transfer of angular momentum and the disk material approaches the compact object (Shakura & Sunyaev 1973, 1976). If the accretion rate exceeds the Eddington limit, ${\dot{M}}_{\mathrm{Edd}}=L_{\mathrm{Edd}}/{\eta }_{\mathrm{rad}}{c}^2$, where L_Edd is the Eddington luminosity and η_rad(∼0.1) is the energy conversion efficiency, the radiation pressure in the innermost part of the disk exceeds the thermal pressure and the disk inflates. In particular, the disk thickness grows linearly with the radius, thereby developing a funnel-like structure. The accretion rate in the innermost region of the disk is bounded at the Eddington rate, while the excess mass flux is ejected in the form of a fast outflow (see e.g. Abaroa et al. 2023, Abaroa & Romero 2024, and references therein). The formation of winds from super-Eddington (or supercritical) accretion was hypothesized by Shakura & Sunyaev (1973) and thoroughly investigated in pioneering work in the 1970s (Meier 1979, 1982a, b, c). Recently, progress in numerical simulations allowed a more detailed comprehension of the wind dynamics such as the typical launching radius, the angular properties of the outflow, and the connection of the kinetic power to the X-ray luminosity (Kitaki et al. 2021).

Ultra-luminous X-ray source activity may result in many sites of interest for particle acceleration and non-thermal processes, such as shocks developed nearby the wind launching radius or other shocks developed at relatively small scales due to the collision between the wind of the donor star and the outflow from the accretion disk (Chatzis et al. 2022; Abaroa et al. 2023). In addition, the accretion onto the compact object might also result in the launching of a highly collimated jet reaching at least mildly relativistic velocities (Blandford & Znajek 1977; Blandford & Payne 1982). As discussed in Globus & Levinson (2016), the existence of wide-opening-angle disk winds and jets, which are typically associated with micro-quasar activity, is not mutually exclusive. Besides, the jet itself is also a relevant site for particle acceleration and emission of non-thermal radiation (H. E. S. S. Collaboration 2024; Abeysekara et al. 2018; Safi-Harb et al. 2022).

In this work, we focused on the wind-blown bubbles inflated by the disk winds, and particularly on the DSA scenario at their wind termination shocks that usually take place at tens of parsecs from the central engine. The simultaneous presence of a disk wind and a collimated jet, even if possible, as attested in the case of SS433 (Watson et al. 1983; Middleton et al. 2021), can be expected to deform and elongate the structure of the bubble along the jet axis, as observed in the case of W50 (Churazov et al. 2024). In our work, we neglected the possible interactions between co-existing disk winds and jets, which may lead to additional sites of particle acceleration, and we focused on the wind bubble on its own in order to understand its multi-messenger signatures.

Figure 1 illustrates a not-to-scale representation of the system where we highlight the key ingredients of a ULX and the associated wind bubble where particle acceleration and transport take place. A compact object such as a stellar-mass black hole (BH) or a neutron star (NS) in super-Eddington accretion develops a puffed part in the innermost region of the accretion disk¹, while a baryon-loaded fast wind with wide opening angle is launched at high speed. Such a wind, as described in Sect. 2.1, is a complex, stratified structure characterized by a strong shock surrounding a fast-wind region. The latter is then surrounded by a hot and extended shell of shocked wind material and by an additional outer thin-shell of swept-up shocked ambient medium. Particles can be accelerated by DSA at such a strong shock and being advected and diffusing throughout the whole bubble. In Sect. 2.2, we highlight the details of the model.

Fig. 1. ULX sketch showing various components (not to scale) that are relevant to our model. Two scenarios for the density distribution in the shocked wind that are investigated are displayed.

2.1. Wind bubble structure

We worked under the assumption that the ULX wind blows approximately steadily in time, with a wide opening angle in a spatially uniform interstellar medium (ISM), thereby developing a bubble structure. A wind bubble forms as follows. At a given time t₀, the ULX starts blowing a fast wind with wide opening angle. The wind is supersonic; consequently, it drives a forward shock (FS) expanding into the external medium. During the expansion the wind sweeps up the external medium, which gets shocked and accumulated in a thin layer behind the FS. Such a layer is separated from the wind material by the contact discontinuity (CD), the physical boundary of the wind. As the system expands, the collision of the wind with the external medium results in the formation of a second shock inside the wind structure, known as wind termination shock (TS), which has its upstream region facing the central engine. We note that, in the reference frame of the central engine, the TS expands radially as the FS does, whereas it travels backward (as a reverse shock) only as seen from the reference frame of the fast wind. After an initial phase of free expansion with velocity V_w, as soon as the swept-up external matter becomes roughly comparable with the wind mass, the system enters the deceleration phase. This happens roughly on a timescale of $t_{\mathrm{dec}}\approx [3\dot{M}/(4\pi m_p{n}_{\mathrm{ext}}{V}_w^3){]}^{1/2}\approx 3\hspace{0.5em}\mathrm{yr}{\dot{\mathrm{M}}}_{-6}^{1/2}{\mathrm{n}}_{\mathrm{ext},0}^{-1/2}{\mathrm{V}}_{\mathrm{w},9}^{-3/2}$, where $\dot{M}$ (${\dot{M}}_{-6}$) is the wind mass-loss rate (in units of 10⁻⁶M_⊙yr⁻¹), n_ext (n_ext,0) is the density of the external medium (in units of cm⁻³), V_w,9 is the terminal wind speed in units of 10⁹ cm s⁻¹, and m_p is the proton mass. In this phase, the FS and the TS expand self-similarly with scaling R_fs∼t^3/5 and R_ts∼t^2/5 (Weaver et al. 1977). In particular, the temporal evolution of the shocks can be written as follows (Koo & McKee 1992):

$$R_{\mathrm{fs}}(t)\approx 120\mathrm{pc}{t}_6^{3/5}({\dot{E}}_{39}{n}_{\mathrm{ext},0}^{-1}{)}^{2/10},$$(1)

$$R_{\mathrm{ts}}(t)\approx 14.5\mathrm{pc}{t}_6^{2/5}({\dot{E}}_{39}{n}_{\mathrm{ext},0}^{-1}{)}^{3/10}{V}_{w,9,}^{-1/2}$$(2)

where t₆ represents the time in units of Myr and ${\dot{E}}_{39}$ is the kinetic power of the outflow ($\dot{M}{V}_w^2/2$) expressed in units of 10³⁹ erg s⁻¹.

As one can see from Eq. (1), the FS rapidly slows down in the ISM. Such a shock is also expected to become radiative as the shocked ambient medium (SAM) typically has a short cooling time. Eq. (2) implies that the TS slows down even faster than the FS, possibly stalling in the reference frame of the central engine. However, as the innermost fast wind launched by the ULX and cooled adiabatically impacts on it at a very high velocity, such a shock is expected to be a strong shock.

The late stage of the bubble, t≫t_dec, is extremely interesting for particle acceleration and transport as the two shocks are evolving very slowly in time, while strong shock conditions are present at the TS. Therefore, DSA can take place efficiently at such a shock, while, approximately, the bubble can be assumed as a steady system. On the other hand, the FS, being radiative and slowing down, is not expected to be an efficient accelerator (see, e.g., del Valle et al. 2022, and references therein). The proper motion of the source might result in the distortion of the bubble shape in the late deceleration (Bosch-Ramon & Barkov 2011). While this can impact the source morphology, it is unlikely to produce an observable impact on the gamma-ray spectrum. Therefore, we leave it for a follow-up investigation.

To sum up, the bubble is a spherically symmetric onion-like structure characterized by the three typical discontinuities: the innermost TS (R_ts) separating fast cool wind from the hot shocked wind, the CD (R_cd) separating the hot shocked wind from the SAM, and the outermost FS (R_fs).

The wind bubble is assumed to be energy-conserving, namely in conditions such that the shocked wind does not radiate. Thus, the wind velocity is u(r)∝r⁻² for r>R_ts. In the innermost region, for r<R_ts the wind velocity is with good approximation constant, with u(r)=V_w. The number density profile is then computed from the mass continuity equation, given the wind velocity profile, as $n_1(r)=\dot{M}/[4\pi m_p{r}^2{V}_w]$ for r<R_ts, where for simplicity we have assumed a composition dominated by ionized hydrogen. For r>R_ts we explored two alternative scenarios:

• –

  S1: The standard thin-shell approximation in which the shocked wind is perfectly separated from the SAM, as represented by the small clumps in Fig. 1 (left hand side), by the CD. In this context, the density of the shocked wind would be constant and equal to n₂=4 n₁(R_ts).

• –

  S2: The target material is uniformly distributed throughout the shocked wind region (thick shell), as represented in the same figure (right hand side) by the uniformly distributed clumps in the shell extending from R_ts up to R_fs. This scenario of clumpy shocked wind (CSW), in which n₂ is a free parameter, can be the result of several dense molecular ISM clumps swallowed by the system, partial evaporation of the SAM (Weaver et al. 1977), or the consequence of instabilities taking place at the contact discontinuity and allowing a sizable fraction of the SAM material to be flooded over by the shocked wind (Bosch-Ramon et al. 2010).

2.2. Particle acceleration and transport in the bubble

We focused on the late-time stages of the wind-bubble evolution, where the shock velocities in the reference frame of the central engine are low compared to V_w. In this context, the particle transport can be approximated as stationary, and the only relevant site for DSA is the TS. On the other hand, at the FS, particles are assumed to freely escape without experiencing any re-acceleration. As mentioned previously, this assumption is justified since the FS is expected to be rapidly slowing down and radiative (Drury 1983).

We assume spherical symmetry. We point out that such an assumption is unjustified in the vicinity of the compact object due to angular dependence of the outflow velocity (Kitaki et al. 2021) and the vicinity of the donor companion star (Abaroa et al. 2023). However, on parsec scales from the launching region it is reasonable to expect these small-scale irregularities to play a negligible role in the overall geometry of the bubble (Menegazzi et al. 2024). Nonetheless, these spatial scales might be potentially interesting for the injection of mechanical turbulence in the fast cold wind. However, a deep investigation of the turbulence properties goes beyond the scope of this work. We modeled particle acceleration and transport of high-energy protons in the wind bubble with the following stationary transport equation (Morlino et al. 2021; Peretti et al. 2022, 2023; Mukhopadhyay et al. 2023):

$$\nabla ·[D\nabla f]-u·\nabla f+\frac{1}{3}p\frac{\partial f}{\partial p}\nabla ·u-\frac{f}{{\tau }_{\mathrm{loss}}}+Q=\frac{\partial f}{\partial t}=0,$$(3)

where f=f(r,p)=dN/d³pdV is the particle distribution function depending on radius r and momentum p. D(r,p) is the diffusion coefficient, u(r) is the velocity profile, and τ_loss(r,p) is the energy-loss timescale for pp interactions. The latter process can be approximated by a catastrophic energy-loss term (i.e., particles of momentum p are removed from the system on a timescale τ_loss instead of moving to lower energies), since the proton in one collision typically loses a significant fraction (κ_pp≈0.5) of its energy (e.g., Kelner et al. 2006). We stress that the stationary approximation is well justified, since in the late deceleration non-thermal particles evolve on timescales shorter than the dynamical time of the system. The injection term Q has the following expression:

$$Q=\eta \frac{n_1(R_{\mathrm{ts}})V_w}{4\pi p^2}\delta [p-p_{\mathrm{inj}}]\delta [r-R_{\mathrm{ts}}],$$(4)

where η is the fraction of the flux of all plasma protons arriving at the TS, n₁(R_ts)V_w, that is injected in the DSA. The two delta functions state that only particles at the wind termination shock position and with a certain momentum p_inj are accelerated. We note that p_inj is generally located in the high-energy tail of the thermal Maxwellian of the plasma, and it does not play any role on the CR pressure provided that it is not assumed to be much larger than m_pc. We assume p_inj=m_pc. Finally, we notice that Takeuchi et al. (2013) and Kobayashi et al. (2018) pointed out that ULX winds might be characterized by a substantial clumpiness, especially close to the launching region. As they are not expected to have a large filling factor over the whole multi-parsec-sized volume, their impact on the transport of particles in the wind bubble is not expected to have an observable impact.

Our model focuses on the hadronic acceleration and radiation, while we limit our discussion of primary electrons and their radio emission to Sect. 4. We also leave the inclusion of the secondary electron-positron populations and their multiwavelength emission to future investigations, as the secondary-to-primary electrons ratio is expected to be smaller than unity for ULXs in the energy range of interest for radio emission. It is worth mentioning that, as the inferred magnetic field in these bubbles is on the order of a few 10 μG, the synchrotron cooling of electrons is at least an order of magnitude more efficient than the inverse Compton cooling on the cosmic microwave background. Consequently, one can expect the leptonic emission to be relevant in radio and X-rays while being subdominant in the gamma-ray domain.

We estimated the magnetic-field pressure (B²/8π) in the fast cool wind as a fraction, ϵ_B, of the wind kinetic energy density ($\rho V_w^2/2$). We computed the magnetic-field amplitude in the downstream region, assuming that only the components perpendicular to the shock normal are compressed, such that the downstream magnetic field is amplified by a factor $\sqrt{11}$ with respect to the upstream one (Marcowith & Casse 2010). The downstream magnetic-field amplitude reads as follows:

$$B=8.5\mu \mathrm{G}{ϵ}_{B,-1}^{1/2}{\dot{M}}_{-6}^{1/2}{V}_{w,9}^{1/2}{R}_{\mathrm{ts},1}^{-1},$$(5)

where R_ts,1 is the TS radius normalized to 10 pc. Such a magnetic field is assumed spatially constant in the downstream. CRs are scattered off magnetic irregularities, such as Alfvén waves. The associated diffusion coefficient is computed in the context of the quasi-linear theory as (Blasi 2013)

$$D(r,p)\approx \frac{v(p)}{3}\{\begin{array}{cc}{r}_L^{2-\delta }(r,p)l_c^{\delta -1},\hfill & r_L\le l_c\hfill \\ r_L^2(r,p)l_c^{-1},\hfill & r_L\ge l_c\hfill \end{array},$$(6)

where v(p) is the particle velocity, r_L is the Larmor radius, l_c is the coherence length of the magnetic field, and δ is the spectral index of the turbulence cascade. We highlight that the diffusion coefficient is assumed to scale as $r_L^2$ when r_L≳l_c (Subedi et al. 2017). In this work, we mostly focused on δ=3/2 as prescribed for magneto-hydrodynamical turbulence (Kraichnan 1965). However, in Sect. 3.2.1, we develop our analytic discussion in the context of Bohm diffusion, δ=1, in order to provide some general upper limits to the maximum energy of accelerated particles.

We solve Eq. (3) following the semi-analytic approach developed in Morlino et al. (2021) and Peretti et al. (2022, 2023). This is based on solving Eq. (3) in the upstream and downstream regions separately and finally joining the two spatial solutions at the wind-shock location. We refer the interested reader to the above-mentioned works for additional information on the solving algorithm.

As it is useful for qualitative understanding, we report the form of the solution at the wind TS:

$$f_{\mathrm{ts}}(p)=C{p}^{-s}{e}^{-{\mathrm{\Gamma }}_1(p)}{e}^{-{\mathrm{\Gamma }}_2(p)},$$(7)

where s≈4 is the power-law index, while the Γ_i (i=1, 2) are two monotonic functions of momentum that set the maximum energy due to the finite size of the upstream and downstream regions. The normalization factor C is computed assuming that the pressure of the accelerated particles is a fraction ξ_CR≲0.1 of the wind ram pressure at the TS so that the shock structure would not be modified by the CR pressure (test particle approximation). In general, for fixed ξ_CR faster winds (with the same $\dot{M}$) favor higher cosmic-ray pressures at the TS. Clearly, ξ_CR and the parameter η in Eq. (4) are different parametrizations of the same physical quantity. Hereafter, we limit our discussion to ξ_CR as its meaning is more intuitive.

We computed the production spectra of gamma rays following the approach described in Kelner et al. (2006). We report the details of this calculation in Appendix A.

Finally, as it can be useful for general considerations of the contribution to the observed CRs, we report the expression of the escaping flux of particles from the FS of a ULX in the following (see also Peretti et al. 2022, 2023):

$$J_{\mathrm{esc}}(p)=\frac{\mathrm{d}N}{{\mathrm{d}}^{3}p\mathrm{d}t}=\frac{\pi R_{\mathrm{ts}}^2{V}_w[1-\mathcal{L}(p)]}{1-exp\{-R_{\mathrm{ts}}{V}_w(1-\frac{R_{\mathrm{ts}}}{R_{\mathrm{fs}}})/[4D_2(p)]\}}{f}_{\mathrm{ts}}(p),$$(8)

where the function ℒ accounts for the reduction to the escaping flux due to the pp catastrophic energy losses. As pp losses have a negligible observational impact on the bulk of accelerated particles, one generally expects ℒ≪1.

3. Results

In what follows, we discuss the results of our model. We remind the reader that, as we are mostly interested in the gamma-ray domain, we regard wind bubbles as hadronic-only sources. In Sect. 3.1, we discuss the spectra of non-thermal accelerated particles characterizing ULXs, while in Sect. 3.2, we focus on the maximum energy achievable by non-thermal particles in these objects.

3.1. Timescales and spectra of high-energy particles

We studied the non-thermal hadronic signatures of ULX bubbles by selecting a set of representative parameters for a prototype source (see Table 1).

Even though the typical timescales of the main physical mechanisms do not have a direct mathematical impact on the solution of our stationary model, they are extremely useful in guiding a qualitative understanding of the dominant processes in different energy ranges. Figure 2 illustrates the comparison of the timescales of the main physical processes taking place at the wind TS. The age of the system (solid black line) is much longer than the acceleration timescale (${\tau }_{\mathrm{acc}}\approx s{D}_1/V_w^2$), which is represented in the plot by the dashed blue line, and always subdominant with respect to the shortest escape time. At low energies, the advection of particles from the TS – which happens on a timescale of τ_adv=(R_fs−R_ts)/〈u₂〉 – is the dominant escape mechanism, whereas at energies of ≳3·10⁷ GeV, the diffusion, whose characteristic timescale is τ_diff=(R_fs−R_ts)²/D₂, takes over and sets an upper limit to the maximum energy. Here, the subscript “2” states that the quantity is computed in the downstream region of the TS. We do not show the adiabatic timescale as, for a downstream wind profile of ∼r⁻², adiabatic energy losses do not play a role in the downstream region. We finally note that the pp loss timescale, τ_pp≈(nκ_ppσ_ppc)⁻¹, is not included in the figure, as it is many orders of magnitude longer than all other timescales. For example, in the context of the thin-shell approximation, the density of the shocked wind can be written as $n_2\approx {10}^{-5}{\dot{M}}_{-6}^{-1}{V}_{w,9}{R}_{\mathrm{ts},1}^2{\mathrm{cm}}^{-3}$, resulting in a loss timescale longer than the age of the Universe. Assuming a higher density, as prescribed for S2, would lead to timescales of the order of 5-50 Myr. Only the density of the SAM layer in S1 might provide a reasonably short timescale. However, the residency time of particles in such a small layer would be reduced by some orders of magnitude due to its reduced radial extension. This means that the ULX bubble is far from being a hadronic calorimeter.

Fig. 2. Typical timescales regulating behavior of non-thermal particles in ULX wind bubble. The acceleration timescale (blue dashed) is compared with the timescales of the main escape processes: advection (yellow dot-dashed) and diffusion (green dotted). The age of the system is also shown as a solid black line.

Figure 3 illustrates the spectral shape of the proton distribution as computed at different radii throughout the downstream region. In particular, the solid black line represents the solution at the shock, while the other line styles and colors represent the proton spectra at progressively larger radii starting from the TS and moving outward to the FS. One can immediately notice that, as suggested from the timescales, the maximum energy² at the shock surpasses the PeV range, reaching the value of E_max≈8 PeV. The radial dependence of the spectrum shows the impact of diffusion, which becomes more and more relevant at larger radii. Finally, as the energy losses in the system are subdominant, the escape spectrum is not substantially different than the spectrum of accelerated particles at the TS.

Fig. 3. Proton spectra at different radii. The spectrum at the TS (solid black) is compared with the spectra at progressively larger radii approaching the FS (from the red dotted to the blue dashed).

3.2. Maximum energy in ULX wind bubbles

The maximum energy that particles can reach at a cosmic accelerator is a fundamental property in the context of the origin of the highest energy both Galactic and extra-Galactic CRs. In what follows, we first provide an analytic argument aimed at understanding the maximum energy accessible to particles accelerated at wind termination shocks of different wind-blown bubbles depending on their main parameters. The qualitative discussion is then followed by a quantitative analysis of the model outcomes and a parameter-space scan.

3.2.1. Analytic argument

In order to support our claim of ULX wind bubbles being new Galactic PeVatron candidates, we provide a simple analytic argument resulting in a general upper limit. Differently from the rest of this manuscript, in this subsection we assume the simplest and most optimistic micro-physical conditions provided by the Bohm diffusion, δ=1.

We derive a physically motivated upper limit on the maximum energy at the TS equating the upstream diffusion length D₁/V_w to a sizable fraction χ of the TS radius R_ts. This results in the following condition:

$$E_{\max }=\frac{3q}{c}\chi {ϵ}_B^{1/2}{\dot{M}}^{1/2}{V}_w^{3/2}\simeq 2.4\hspace{0.5em}\mathrm{PeV}\chi {ϵ}_{B,-1}^{1/2}{\dot{M}}_{-6}^{1/2}{V}_{w,9}^{3/2}.$$(9)

We highlight that while the maximum energy has a strong dependence on the wind speed, it is only mildly dependent on the mass-loss rate. Eq. (9) depends on a few parameters of the system allowing a direct qualitative comparison with other Galactic sources characterized by wind bubbles, such as those inflated by young massive stellar clusters (YMSCs). Figure 4 illustrates the parametric dependence of the maximum energy, where we assume χ=1 and ϵ_B=0.1. Interestingly, differently from the standard Hillas plot where the maximum energy is set by the source size and magnetic field, here the parameter space can be reduced to $\dot{M}$ and V_w, which are the main macro-physical parameters of these objects. The range of parameters shown in the figure are taken from King et al. (2023) and Stevens & Hartwell (2003) for ULXs and YMSCs, respectively. The dashed lines identify three different maximum energies given by Eq. (9).

Fig. 4. $\dot{M}-V_w$ parameter space of typical Galactic wind bubbles and achievable maximum energy computed according to Eq. (9).

We remind the reader that these contours were computed assuming the most optimistic conditions for the bubble, set by the Bohm diffusion. As a result, Eq. (9) should be interpreted as an upper limit to the true maximum energy. Our results suggest that only the most extreme YMSCs could possibly access the PeV range; whereas, because of their faster outflows, ULXs stand as an extremely promising class of Galactic PeVatrons.

3.2.2. Parameter-space scan

In our wind-bubble model, the power of the system is set by $\dot{M}$ and V_w, while ϵ_B and l_c are related to the properties of the magnetic field and affect the diffusion of particles. While the mass-loss rate and terminal wind speed can be inferred from observations, the magnetic field is typically unknown, thereby making its associated parameters highly unconstrained. In particular, the coherence length could extend from a small up to a sizable fraction of the size of the system. Similarly, the magnetic-to-ram pressure ratio could range from ∼10% down to a few percent (and possibly lower values). For each of the four mentioned parameters, we set up a variability range and explored the resulting maximum energy while keeping the age of the system to 1 Myr constant.

Table 2 illustrates the impact of the four parameters on E_max (expressed in PeV). In particular, we explored three possible values of l_c (0.01 pc, 0.1 pc, 1 pc) and ϵ_B (0.01, 0.03, 0.1) and two values of $\dot{M}$ (10⁻⁶–10⁻⁵ M_⊙ yr⁻¹) and V_w (10⁴–10⁵ km s⁻¹). It is possible to observe that the qualitative trend expressed by Eq. (9) is approximately confirmed for the three parameters characterizing such relation, even though the diffusion scenario adopted here, δ=3/2, is different from the one of the simple analytic argument leading to Eq. (9). In particular, in order to give a quantitative idea, the normalization at 10 GeV of the upstream diffusion coefficient at the shock can be written as $D_1(10\mathrm{GeV})\approx 2·{10}^{25}{\mathrm{cm}}^2{\mathrm{s}}^{-1}{R}_{\mathrm{ts},1}^{1/2}{l}_{\mathrm{c},-1}^{1/2}{ϵ}_{\mathrm{B},-1}^{-1/4}{\dot{M}}_{-6}^{-1/4}{V}_{w,9}^{-1/4}$, where the lengths are expressed in decades of parsecs.

Conversely to the other parameters, it is possible to notice a nontrivial dependence of the E_max on l_c, which seems to have a local maximum around l_c≈10⁻¹ pc at low velocities. This behavior can be ascribed to a transition in the diffusion regime. In fact, in the majority of realizations, the maximum energy is already reached in the small pitch angle scattering regime, where $E_{\max }\sim l_c^{1/2}$. While, for l_c≈1 pc and V_w=10⁴ km s⁻¹, the maximum energy occurs in the resonant scattering branch of the diffusion coefficient where $E_{\max }\sim l_c^{(1-\delta )/(2-\delta )}$.

The parameter space scan combined with the previous analytic argument suggests that typical ULXs can reach the PeV range for a wide range of mass-loss rates and terminal wind speeds. In particular, the maximum energy is an increasing function of these macroscopic parameters. The magnetic-field amplitude, parametrized by ϵ_B, also influences the maximum energy monotonically, whereas the dependence of the latter on the coherence length is nontrivial as the maximum energy might occur in different regimes of the diffusion coefficient. Finally, we noticed that particularly powerful ULXs could access several tens of PeV.

4. Multi-messenger emission from ULXs: Application to SS 433

The Galactic source SS 433 has been identified as a plausible misaligned (hidden) Galactic ULX (see, e.g., Begelman et al. 2006 and King et al. 2023) that is powered by a super-Eddington accretion onto a compact object (most likely a black hole, as discussed in, e.g., Poutanen et al. 2007 and Middleton et al. 2021). SS 433 has been observed in very high-energy gamma-rays (E_γ≳1 TeV) by H.E.S.S. (H. E. S. S. Collaboration 2024), HAWC (Abeysekara et al. 2018), and more recently LHAASO (LHAASO Collaboration 2024). In this section, we present our model predictions for the case of SS 433 using the representative set of parameters listed in Table 1 as inferred to be appropriate for such an object (King et al. 2023).

As the energy losses are not expected to affect the proton distribution (see discussion about Fig. 2), differences in the gamma-ray spectrum can only arise from different distributions of target material in the bubble for a fixed set of parameters. Fig. 5 illustrates the difference in the gamma-ray spectra obtained in the two scenarios, S1 (green curves) and S2 (blue curves), describing the downstream region of the TS (see sketch in Fig. 1). In particular, for S2 we assumed a downstream density of n_CSW=0.1 cm⁻³. For each scenario the shaded area represents the model predictions obtained by varying the terminal wind speed from 0.1 c to 0.2 c while keeping the cosmic-ray efficiency at the shock, ξ_CR=0.1, constant. The corresponding maximum energy is found to range from 3 PeV up to 8 PeV.

Fig. 5. Spectral energy distribution of gamma-ray emission from SS 433. Colored markers indicate flux measurements obtained with H.E.S.S. and HAWC. More recent data obtained by LHAASO are shown with black markers. The upper limits imposed by Fermi-LAT on the GeV emission from the eastern (E) and western (W) jets are overplotted with colored downward-pointing arrows. The predicted (hadronic) gamma-ray spectrum from the bubble is overplotted for the two scenarios presented in Fig. 1. The shaded region shows the expected range of gamma-ray emission when varying the terminal wind speed by a factor of two.

The model predictions for the gamma-ray emission from the whole bubble were visually compared with HAWC (Abeysekara et al. 2018) and LHAASO (LHAASO Collaboration 2024) data, as well as the Fermi-LAT upper limits and H.E.S.S. data (H. E. S. S. Collaboration 2024) collected from the eastern (E) and western (W) jets launched from the compact object powering the whole system. We note that as the model prediction is produced from the whole bubble embedding both jets, its contamination to the E and W lobes observed by H.E.S.S. should be scaled down by a factor accounting for the different angular extensions of the two regions. Consequently, the best impression of the actual contribution of the bubble to the overall emission can be deduced by comparing the model prediction with the LHAASO data.

Interestingly, unless the most optimistic scenario is assumed, the wind-bubble model prediction can provide a dominant contribution – of hadronic origin – to the observed emission only above 100 TeV. It is intriguing that, in such a range, leptonic models are expected to fall more rapidly than the observed spectrum (see, e.g., Fig. 1 in LHAASO Collaboration 2024). Moreover, the morphology of the gamma-ray source changes as a function of photon energy, with lower energy emission (1−100 TeV) following the jet geometry; whereas, the highest energy emission (>100 TeV) appears to be more spherical (see Fig. 1 in LHAASO Collaboration 2024). The angular extent of the >100 TeV emission translates to a minimum radius of r₃₉≈30 pc (taking the distance to the source to be 5.5 kpc in agreement with Blundell & Bowler 2004 and Lockman et al. 2007). Interestingly, using Eqs. (1) and (2), we find that R_ts≃3 pc and R_fs≃35−47 pc for the set of parameters assumed. These results support the bubble interpretation and warrant further investigation.

Finally, it is intriguing to notice that Fermi-LAT might soon be able to constrain the bubble emission in the GeV range from the hadronic scenarios discussed in this section. In particular, the detection of a relatively soft emission in the GeV band (softer than the rising part of an inverse Compton component) combined with a morphology similar to the LHAASO emission >100 TeV would be an additional strong indication supporting the bubble interpretation.

As a candidate Galactic PeVatron, SS 433 is expected to produce a high-energy neutrino flux comparable with the gamma-ray one. In particular, in the multi-pion regime one expects the production rates of neutrinos and gamma rays to scale as $2E_{\nu }^2{Q}_{\nu }(E_{\nu })\approx 3E_{\gamma }^2{Q}_{\gamma }(E_{\gamma })$, where E_γ≈2E_ν. Consequently, in the TeV range one can expect the neutrino flux from SS 433 to be on the order of $E_{\nu }^2{\mathrm{\Phi }}_{\nu }^{\mathit{SS}433}(E_{\nu })\approx {10}^{-12}-{10}^{-13}\mathrm{TeV}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}$. It is interesting to explore how upcoming neutrino observatories located in strategic positions to observe Galactic objects, such as KM3NeT (KM3NeT Collaboration 2024), will be able to constrain the possible hadronic emission from the wind bubble inflated by SS 433. Following KM3NeT Collaboration (2024), we considered the 12-year discovery potential of KM3NeT-ARCA for a point-like source with E⁻² spectrum: ${\mathrm{\Phi }}_{\mathrm{DP}}^{\mathit{PL}}\approx 5·{10}^{-10}{\mathrm{GeV}}^{-1}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}$. Moreover, we assumed σ_PSF=0.2deg as a representative uncertainty for track-like events (Ambrosone et al. 2024). This allowed us to estimate the discovery potential of the extended SS 433 wind bubble as

$$E_{\nu }^2{\mathrm{\Phi }}_{\mathrm{DP}}^{\mathrm{SS}\hspace{0.5em}433}(E_{\nu })\approx E_{\nu }^2{\mathrm{\Phi }}_{\mathrm{DP}}^{\mathit{PL}}(E_{\nu })\sqrt{\frac{{\sigma }_{\mathrm{PSF}}^2+{\sigma }_{\mathrm{src}}^2}{{\sigma }_{\mathrm{PSF}}^2}}\approx {10}^{-12}\mathrm{TeV}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1},$$(10)

where we assumed σ_src=r₃₉≈0.32deg. This suggests that during its first decade KM3NeT might already be able to detect the bubble of SS 433 if the actual flux from the source is at the level of the most optimistic realizations shown in Fig. 5.

The radio domain can be an interesting alternative window for testing particle acceleration taking place at the TS of ULX bubbles. In the following, we present an analytic argument providing an order-of-magnitude estimate of the low-energy leptonic emission from primary electrons accelerated in these systems.

We began our estimation assuming a primary electron population at the shock ${\mathit{dn}}_e/d{\gamma }_e\sim {\gamma }_e^{2-s}\exp [-{\gamma }_e{m}_{\mathit{ec}}^2/E_{\max ,\mathrm{e}}]$. We assumed s=4 as prescribed by test-particle DSA, and we estimated the maximum energy by comparing the acceleration with the radiative-loss timescale of electrons. The latter is governed by synchrotron losses, which dominate over the inverse Compton cooling on the CMB, since $U_B\simeq 4·{10}^{-12}{B}_1^2{\mathrm{erg}\mathrm{cm}}^{-3}\simeq 10U_{\mathrm{CMB}}$, where B₁ is the magnetic field expressed in units of 10 μG (typical value for the shocked wind bubble; see Eq. 5). The synchrotron-loss timescale then reads

$${\tau }_{\mathrm{syn},\mathrm{e}}\approx 0.2\hspace{0.5em}\mathrm{Gyr}{B}_1^{-2}{\gamma }_{e,3}^{-1},$$(11)

and the maximum electron energy achieved is

$$E_{\max ,\mathrm{e}}\approx 0.5V_{w,9}^{2/(3-\delta )}{l}_{c,-1}^{(1-\delta )/(3-\delta )}{B}_1^{-\delta /(3-\delta )}\mathrm{PeV}.$$(12)

Unless different background photon fields or stronger magnetic fields are present, ULXs can also be electron PeVatrons. Therefore, for our qualitative argument we assume E_max,e=1 PeV.

As we are interested in estimating the flux at radio wavelengths, we focused on relatively low electron energies. Using ${\nu }^{*}=1\mathrm{GHz}{\nu }_9^{*}$ as the observing frequency, we estimated the typical Lorentz factor of radiating electrons. Following Rybicki & Lightman (1986), we obtain ${\gamma }_e^{*}\approx 6·{10}^3{\nu }_9^{*1/2}{B}_1^{-1/2}$ (or $E_e^{*}\approx 3{\nu }_9^{*1/2}{B}_1^{-1/2}\mathrm{GeV}$). Since the cooling timescale of these electrons (see Eq. 11) is much longer than their advection timescale, f_e,ts is, with good approximation, spatially constant throughout the whole wind bubble. This also suggests a radial morphology of the radio emission.

Using a delta-function approximation for the single-particle synchrotron emissivity (Ghisellini 2013), j_S(ν)=j₀δ(ν−ν_S), where ${\nu }_S=\mathit{eB}{\gamma }_e^2/(2\pi m_{e}c)$ and $j_0=(4/3){\sigma }_{T}c{U}_B{\gamma }_e^2$, we may write the synchrotron flux at ν as

$$F_S(\nu )\approx j_S(\nu )\frac{\mathrm{d}{n}_e}{\mathrm{d}{\gamma }_e}\frac{\mathrm{d}{\gamma }_e}{\mathrm{d}\nu }\frac{\frac{4}{3}\pi R_{\mathrm{fs}}^3}{4\pi D_L^2},$$(13)

where the emitting volume is taken to be equal to the volume of the shocked wind with $R_{\mathrm{fs}}^3\gg R_{\mathrm{ts}}^3$, and D_L is the luminosity distance of the source. Assuming, as in the case of protons, that non-thermal electrons take a fraction, ξ_e≈10⁻²ξ_CR, of the ram pressure at the shock, the synchrotron flux can be rewritten as

$$\begin{array}{cc}\hfill {\nu }^{*}{F}_S({\nu }^{*})& \approx \frac{{\sigma }_T}{6\pi m_e^2{c}^3}\frac{E_e^{*}{U}_B}{D_L^2}\frac{{\xi }_e\dot{M}{V}_w}{{\mathrm{\Lambda }}_e}\frac{R_{\mathrm{fs}}^3}{R_{\mathrm{ts}}^2}\hfill \\ \hfill & \approx {10}^{-14}\frac{E_e^{*}{B}_1^2}{D_{L,3.74}}\frac{{\xi }_{e,-3}{\dot{M}}_{-6}{V}_{w,9}}{{\mathrm{\Lambda }}_{e,1.32}}\frac{R_{\mathrm{fs},2}^3}{R_{\mathrm{ts},1}^2}\mathrm{erg}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1},\hfill \end{array}$$(14)

where Λ_e=ln[E_max,e/(m_ec²)]≈21 and D_L is normalized to 10^3.74 pc≈5.5 kpc, the distance of SS 433. Interestingly, our order-of-magnitude prediction is in good agreement with the NVSS observation of the bubble at around 1.4 GHz (Condon et al. 1998). This could be a hint for a substantial contribution of the wind bubble to the spherical component of the nebula W50 surrounding SS 433 (see Churazov et al. 2024 for a detailed presentation of a model including both an isotropic and collimated wind).

Although the association of the spherical part of W50 and at least part of its radio emission with a supernova remnant cannot be discarded, the qualitative argument presented in this subsection is a plausible alternative and illustrates that a leptonic population in the wind bubble can already be probed through radio observations. A more detailed comparison with existing radio data would require the full solution of Eq. (3) and account for the secondary electron population injected through the decay of charged pions resulting from pp interactions. We leave this investigation for a follow-up study.

5. ULXs as Galactic PeVatrons

The maximum energies achievable at the TS of ULX bubbles makes them extremely interesting PeV accelerators. For this reason, we also explore the role of ULXs as sources of Galactic CRs in the knee region.

The possible contribution of ULXs to the CR flux at Earth could be relevant if at least one of the following two scenarios is true: (1) there are enough ULXs in our Galaxy to allow their escaping CR flux to diffuse and fill the whole Galactic volume; (2) the Earth is close enough to a ULX to allow CRs to reach us without escaping the Galaxy.

We begin by considering the first scenario. According to Kovlakas et al. (2020), the average number of observable ULXs in a galaxy can be inferred from astronomical observations as follows:

$${\mathcal{N}}_{\mathrm{ULX}}^{(\mathrm{obs})}=0.{45}_{-0.09}^{+0.06}\times \frac{\mathrm{SFR}}{{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}}+3.3_{-3.2}^{+3.8}\times \frac{M_{★}}{{10}^{12}{\mathrm{M}}_{\odot }},$$(15)

where SFR is the star-formation rate and M_★ is the stellar mass of the Galaxy. Using the typical SFR and stellar mass of the Milky Way (Licquia & Newman 2015), the expected number of observable ULXs in our Galaxy is on the order of unity. We highlight that the total number of ULXs can be larger than the observed one as these sources are geometrically beamed. In particular, assuming a standard opening angle of θ_f=20deg for the funnel, the fraction of observable ULXs will be rescaled with the factor ℱ=(1−cos θ_f)≈0.06 (see also Lasota & King 2023 and Kayanikhoo et al. 2025 for comparable estimates).

Bearing in mind that the expected number of observable ULXs in the Galaxy should be extremely limited, while the actual number could be larger, we evaluate the total number of ULXs necessary to contribute to the CR flux at PeV energies (i.e., at the knee of the CR spectrum) in the following. In order to do so, we assumed isotropic diffusive transport in the Galaxy with the Galactic diffusion coefficient $D_{\mathrm{Gal}}(E)=D_0{E}_{\mathrm{GeV}}^{{\delta }_G}=3·{10}^{28}{E}_{\mathrm{GeV}}^{{\delta }_G}{\mathrm{cm}}^2{\mathrm{s}}^{-1}$, where δ_G is the energy dependence of the Galactic diffusion coefficient (see, e.g., Aguilar 2016 and Génolini et al. 2017). The energy flux of CRs in stationary conditions in the Galaxy can be written as

$$E^2{\mathrm{\Phi }}_{\oplus }(E)=\frac{H_{\mathrm{Gal}}^2}{D_{\mathrm{Gal}}(E)}\frac{c}{4\pi }{\mathcal{N}}_{\mathrm{ULX}}\frac{E^2{J}_{\mathrm{esc}}(E)}{V_{\mathrm{Gal}}},$$(16)

where ${\mathcal{N}}_{\mathrm{ULX}}$ is the total number of ULXs in the Galaxy, $V_{\mathrm{Gal}}=2\pi H_{\mathrm{Gal}}{R}_{\mathrm{Gal}}^2$ is the Galactic volume, R_Gal the Galactic radius, and H_Gal≈4 kpc is the Galactic halo size (see, e.g., Evoli et al. 2018). Under the assumption of s=4 and neglecting the specific shape of the exponential high-energy cut-off, the escaping CR power from a ULX, E²J_esc(E), can be written as

$$E^2{J}_{\mathrm{esc}}(E)\approx \frac{3{\xi }_{\mathrm{CR}}}{2\mathrm{\Lambda }}{\dot{E}}_{\mathrm{kin}},$$(17)

where ξ_CR is the ratio between the cosmic ray pressure and ram pressure at the wind shock, for which we expect something on the order of 1−10%, and Λ≈ln(E_max/E_inj)≈14. One can see that Eq. (17) was obtained assuming an E⁻² spectral shape extending up to the maximum energy and neglecting the specific form of the high-energy cut-off.

Knowing that E²Φ_⊕(1PeV)≈10⁻⁴ GeV cm⁻² s⁻¹ sr⁻¹ (see, e.g., Recchia & Gabici 2024, and references therein), it is possible to invert Eq. (16) and write the total number of ULXs as

$${\mathcal{N}}_{\mathrm{ULX}}\approx 20·[E^2{\mathrm{\Phi }}_{\oplus }{]}_{-4}{D}_{28.5}{H}_{\mathrm{Gal},0.7}^{-1}{R}_{\mathrm{Gal},1.2}^2{\xi }_{\mathrm{CR},-1}^{-1}{\mathrm{\Lambda }}_{1.2}{\mathrm{\Delta }}_G{\dot{E}}_{\mathrm{kin},39}^{-1},$$(18)

where the kinetic power is computed assuming the benchmark values of $\dot{M}$ and V_w (see Table 1), and Δ_G=10^6δ_G−2. This suggests that, while remaining consistent with the order-of-magnitude prediction of their observable number in the Galaxy (see also Mapelli et al. 2010 and Prestwich et al. 2013), Galactic ULXs might meet the energy requirements to substantially contribute to the observed CR flux at the knee. Their relative reduced number would indeed be consistent with their nondetection, ${\mathcal{N}}_{\mathrm{ULX}}^{(\mathrm{obs})}\approx \mathcal{F}{\mathcal{N}}_{\mathrm{ULX}}$, in that all sources might simply be misaligned with favorable lines of sight or highly obscured. Assuming that ULXs are uniformly distributed in the circular Galactic disk, the average distance between them would be $d_{\mathrm{ULX}}\approx 128R_{\mathrm{Gal}}/[45\pi \sqrt{{\mathcal{N}}_{\mathrm{ULX}}}]\approx 3R_{\mathrm{Gal},1.2}{\mathcal{N}}_{\mathrm{ULX},1.3}^{-1/2}\mathrm{kpc},$ making our estimate plausible but somehow close to a critical threshold where its applicability may become questionable as d_ULX≲H_Gal. In fact, as the average distance between ULXs is expected to be on the order of a couple of kiloparsecs, it is reasonable to expect that the flux at Earth might be dominated by one, or at most two, nearby sources. This leads us to consider the second scenario.

Even though there are no ULXs detected in the Milky Way according to the observational definition of L_X≳10³⁹erg s⁻¹, there may be super-Eddington accreting sources with their funnels pointing away from our line of of sight, making them hidden in X-rays. In fact, according to their current state of accretion, SS 433 (King et al. 2023) and Cygnus-X3 (Veledina et al. 2024a) are considered (hidden) Galactic ULXs, and other objects, such as 4FGL J1405.1-6119 (Saavedra et al. 2023), might soon be added to the list. As currently there are no other closer ULX candidates in our Galaxy, it is reasonable to assume that SS 433 is the closest one. However, its distance (∼5.5 kpc), which is larger than the inferred height of the Galactic magnetized halo, implies that for isotropic diffusion, most of its CR flux would diffuse away from the Galaxy before reaching us. Therefore, unless peculiar transport conditions enhance the escaping flux of CRs from SS 433, this source seems unlikely to dominate the flux at the knee. On the other hand, if a highly obscured and undetected ULX were relatively close to Earth, its contribution to the CR spectrum at the knee could be substantial.

6. Conclusions

In this work, we explored the power of diffusive shock acceleration in the wind bubbles excavated by ULXs. We did so by solving the space and energy-dependent cosmic-ray transport equation, and particularly by focusing on the particle acceleration at the wind termination shocks of these bubbles. We studied the maximum energy of these accelerators and discussed the associated multi-messenger radiation in terms of gamma rays, high-energy neutrinos, and radio. We finally explored the possible role of ULXs as Galactic cosmic-ray sources.

We found that for a wide range of typical parameters, ULXs can accelerate protons up to several peta-electronvolts, while the most powerful sources could energize up to tens of peta-electronvolts. This suggests that ULX wind bubbles can be considered as a new class of PeVatrons.

We applied our model to the super-Eddington accreting Galactic source SS 433 where we found that the maximum energy can be as high as 3-8 PeV and the associated wind bubble could be already detectable in gamma rays as a diffuse GeV component. The bubble emission might also explain the flux of gamma-rays observed by LHAASO at energies >100 TeV as well as the change in the gamma-ray source morphology taking place at around that energy. In particular, we argued that below 100 TeV the collimated jets dominate the emission, whereas at higher energies, where the inverse Compton is expected to fade down, the hadronic bubble emission can be observed.

We noticed that in LHAASO Collaboration (2024), the hadronic interpretation is supported by a spatial coincidence of the emission (>100 TeV) with a giant molecular cloud. As the real position of such a cloud with respect to the geometrical boundaries of the bubble is uncertain, we leave a detailed study of the bubble-cloud modeling to a follow-up investigation. However, we point out that if the cloud was embedded inside the forward shock radius, the system might already be qualitatively described by our scenario S2. On the other hand, if the cloud was situated at larger distances than the forward shock, a proper diffusion model would be applied to the escaping particles to infer the actual cloud luminosity.

We estimated the high-energy neutrino flux from the wind bubble of SS 433, and we concluded that a neutrino observatory such as KM3NeT could detect a flux or constrain some realizations of our model with about ten years of observations. Moreover, our order-of-magnitude estimate of the radio flux at 1 GHz (due to synchrotron radiation of accelerated electrons in the bubble) is comparable with the total radio flux in the same energy band as measured from the nebula W50 that surrounds SS 433. This result led us to speculate on the possible ULX-bubble nature of the spherical part of W50, instead of interpreting it as a supernova remnant. In particular, the presence of a mildly relativistic disk wind (Middleton et al. 2021) might substantially increase former estimates (Chi et al. 2024) of its impact on the overall structure of the nebula.

We finally considered the possible role of ULXs as Galactic PeVatrons in contributing to the cosmic-ray spectrum at the knee. According to an observational scaling based on the star-formation rate, the number of observable ULXs expected in our Galaxy should be on the order of unity, while, according to our estimates, about 20 of these objects would be required to saturate the CR flux at the knee. Since ULXs are geometrically beamed, the two numbers are consistent with each other as all ULXs might be not favorably aligned with our line of sight. However, such a small number of ULXs in the Galaxy would result in one, or at most two, nearby objects dominating the observed flux. With our current knowledge, given the distance of the nearest Galactic ULX candidate, SS 433, the conclusion shall be that ULXs are unlikely to dominate the CR flux at the knee. However, we argue that in addition to Cygnus X-3 and SS 433, in agreement with our estimate, there could be other obscured ULXs in the Galaxy. We conclude that ULXs cannot be discarded as dominant cosmic-ray sources at the knee, as some other undetected, or unidentified, ULXs could be hosted in our Galaxy within a few kiloparsecs of the Earth.

Ultra-luminous X-ray sources are powerful sources powered by super-Eddington accretion onto stellar-mass compact objects. Our investigation highlighted their multi-messenger aspects and showed that they can be sources of PeV cosmic rays, high-energy gamma rays, radio, and neutrinos in the Galaxy.

Acknowledgments

E.P. is grateful to R. Amato for insightful discussions and to G. Israel for useful comments on the manuscript. E.P. is also grateful to V. Bosh-Ramon, S. Recchia, C. Pinto and L. Crosato Menegazzi for valuable comments that allowed improving the first version of the manuscript, and E. Amato, G. Morlino, N. Bucciantini, D. Caprioli, C. Evoli, A. Condorelli and A.Ambrosone for insightful discussions and comments. M.P. and G.V. would like to thank Université Paris Cité, where this project was conceived, for their hospitality. E.P. and S.G. were supported by Agence Nationale de la Recherche (grant ANR-21-CE31-0028). M.P. acknowledges support from the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the “2nd call for H.F.R.I. Research Projects to support Faculty members and Researchers” through the project UNTRAPHOB (Project ID 3013). G.V. also acknowledges support by the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the “3rd Call for H.F.R.I. Research Projects to support Postdoctoral Researchers” through the project ASTRAPE (Project ID 7802).

References

Appendix A: Gamma rays from pp interactions

We compute the production spectra of high-energy gamma rays (E_γ≥100 GeV) following the approach described in Kelner et al. (2006):

$$\begin{array}{cc}\hfill \frac{{\mathit{dN}}_{\gamma }(E_{\gamma })}{{\mathit{dE}}_{\gamma }}& =c{\displaystyle {\int }_0^{R_{\mathrm{fs}}}}\mathit{dr}4\pi r^{2}n(r){\displaystyle {\int }_{p_{\min }}^{\infty }}\frac{\mathit{dp}}{E_p(p)}4\pi p^2\times \hfill \\ \hfill & f(r,p){\sigma }_{\mathrm{inel}}(p)F_{\gamma }\left(\frac{E_{\gamma }}{E_p(p)},E_p(p)\right),\hfill \end{array}$$(A.1)

where $p_{\min }=\sqrt{E_{\gamma }^2/c^2-m_p^2{c}^2}$, while σ_inel and F_γ are given in Eqs. (79) and (58) of Kelner et al. (2006). Low energy (E_γ<100 GeV) gamma rays are computed following the delta-function approximation described in Kelner et al. (2006):

$$\frac{{\mathit{dN}}_{\gamma }(E_{\gamma })}{{\mathit{dE}}_{\gamma }}=2{\displaystyle {\int }_0^{R_{\mathrm{fs}}}}\mathit{dr}4\pi r^2{\displaystyle {\int }_{E_{\min }}^{\infty }}\frac{q_{\pi }(r,E_{\pi })}{\sqrt{E_{\pi }^2-m_{\pi }^2{c}^4}}{\mathit{dE}}_{\pi },$$(A.2)

where $E_{\min }=E_{\gamma }+m_{\pi }^2{c}^4/4E_{\gamma }$, m_π is the pion mass and the production rate of π⁰ mesons reads

$$\begin{array}{cc}\hfill q_{\pi }(E_{\pi })=& \frac{\tilde{n}c}{K_{\pi }}n(r){\sigma }_{\mathrm{inel}}\left[{\left(m_{p}c+\frac{E_{\pi }}{c{K}_{\pi }}\right)}^2-m_p^2{c}^2\right]\times \hfill \\ \hfill & 4\pi p^2\frac{\mathit{dp}}{{\mathit{dE}}_p}f\left(r,\left[{\left(m_{p}c+\frac{E_{\pi }}{c{K}_{\pi }}\right)}^2-m_p^2{c}^2\right]\right)\hfill \end{array}$$(A.3)

where $\tilde{n}$ and K_π are free parameters guaranteeing continuity with the high-energy branch (Eq. (A.1)).

All Tables

All Figures

	Fig. 1. ULX sketch showing various components (not to scale) that are relevant to our model. Two scenarios for the density distribution in the shocked wind that are investigated are displayed.
In the text

	Fig. 2. Typical timescales regulating behavior of non-thermal particles in ULX wind bubble. The acceleration timescale (blue dashed) is compared with the timescales of the main escape processes: advection (yellow dot-dashed) and diffusion (green dotted). The age of the system is also shown as a solid black line.
In the text

	Fig. 3. Proton spectra at different radii. The spectrum at the TS (solid black) is compared with the spectra at progressively larger radii approaching the FS (from the red dotted to the blue dashed).
In the text

	Fig. 4. $\dot{M}-V_w$ parameter space of typical Galactic wind bubbles and achievable maximum energy computed according to Eq. (9).
In the text

Fig. 5. Spectral energy distribution of gamma-ray emission from SS 433. Colored markers indicate flux measurements obtained with H.E.S.S. and HAWC. More recent data obtained by LHAASO are shown with black markers. The upper limits imposed by Fermi-LAT on the GeV emission from the eastern (E) and western (W) jets are overplotted with colored downward-pointing arrows. The predicted (hadronic) gamma-ray spectrum from the bubble is overplotted for the two scenarios presented in Fig. 1. The shaded region shows the expected range of gamma-ray emission when varying the terminal wind speed by a factor of two.

In the text