© The Authors 2024

1. Introduction

Microquasars (MQs) are binary systems in which a compact object (such as a black hole or neutron star) draws matter from its companion star, leading to the formation of an accretion disk. Bipolar and highly collimated outflows emerge from the central region of these sources. Since the observation of the first MQs three decades ago (Mirabel & Rodríguez 1994), many of these systems have been discovered and studied, exhibiting a wide variety of characteristics (e.g., Mirabel & Rodríguez 1999; Bosch-Ramon et al. 2005b; Romero & Vila 2008; Belloni 2010; Bordas et al. 2015; Basak et al. 2017; Romero et al. 2017; Zdziarski et al. 2017; Koljonen et al. 2018; Middleton et al. 2021; Dexter & Begelman 2024).

S26 is a unique MQ located in the Sculptor Galaxy NGC 7793 at a distance of 3.9 Mpc (Karachentsev et al. 2003), whose nebula has a size of ∼350 × 185 pc (Pannuti et al. 2002; Dopita et al. 2012). According to observations in the optical/X-ray band, and following the relation of Weaver et al. (1977) between the wind-driven stellar bubbles and the jet-driven radio lobes, Pakull et al. (2010) showed that this MQ has the most powerful jets in accreting binaries, with a mechanical luminosity of L_j ∼ 5 × 10⁴⁰ erg s⁻¹. These authors also estimated its age (t ∼ 2 × 10⁵ yr), proposed that the binary system at the core of S26 is composed of a black hole (BH) and an early Wolf-Rayet (WNE) star, and identified nonthermal X-rays produced near the BH. Soria et al. (2010) used the Australia Telescope Compact Array (ATCA) to resolve the radio lobe structure of S26, and suggested that the radio emission from the terminal region of the jets is consistent with synchrotron mechanisms. They also analyzed X-ray observations from the Chandra space telescope, identifying X-ray hot spots ≈20 pc farther out than the peak of radio intensity in the lobes, and argued that this emission is most likely of optically thin thermal origin. Other authors have constrained the properties of the expanding nebula of S26 by spectroscopic analysis (Dopita et al. 2012) and studied its high-energy and neutrino emission (Inoue et al. 2017).

A key parameter for characterizing MQs is the accretion rate of matter onto the compact object. Accretion proceeds in three basic regimes, depending on the ratio of the actual accretion rate to the Eddington rate. At very sub-Eddington rates, the inner disk is optically thin and geometrically thick, while at moderate rates the regime approaches the standard Shakura-Sunyaev disk, which is optically thick and geometrically thin. In super-accreting BHs (Abramowicz et al. 1980; Fukue 2004), the radiation pressure overcomes gravity in the innermost part of the disk, and its surface layers are ejected as a consequence of the self-regulation of accretion at the Eddington rate (Fukue 2009). In these systems, the accretion disk launches powerful winds with mass loss rates similar to the accretion rate (Abaroa et al. 2023; Abaroa & Romero 2024a). An example of a super-Eddington MQ in our Galaxy is SS433, whose jets have a kinetic power of ≈10³⁹ erg s⁻¹ (see Fabrika 2004, for a detailed review) and presents opaque equatorial outflows (e.g., Picchi et al. 2020). Following the jet-disk symbiosis model (Falcke & Biermann 1995), the kinetic luminosity of the jet of S26 estimated by Pakull et al. (2010) implies that the accretion luminosity to the stellar-mass BH should be super-Eddington, on the order of L_acc ∼ L_j. However, the observed X-ray flux from the core of the system, measured by the Chandra and XMM-Newton telescopes, indicates an apparent very sub-Eddington luminosity of L_X ∼ 10³⁷ erg s⁻¹, orders of magnitude smaller than the jet power.

In this paper we study the properties of S26 and propose an explanation for the apparent discrepancy between the disk and jet luminosities. We hypothesize that the compact object is a 10 solar mass BH accreting at super-Eddington rates, and that its actual accretion power is a few times 10⁴⁰ erg s⁻¹. We suggest that some energy should also be extracted from the ergosphere of the spinning BH to provide enough power for the jet (Blandford & Znajek 1977). In such a scenario, the efficiency of converting accretion power into jet power can exceed 100 percent in highly rotating BHs, reaching transient values of 300 percent in the most extreme cases (Tchekhovskoy et al. 2010, 2011; McKinney et al. 2012). On the other hand, as mentioned above, one of the main features of these supercritical systems is a very dense wind ejected from the accretion disk.

In our hypothesis, the opaque wind almost completely obscure the disk, making it impossible for an observer to see the X-rays produced in its innermost region. If this is the case, then the X-ray flux detected from the core would not be coming from the disk, but from some other region of the system. This emission could be produced at the base of the jet by internal shocks, close to the BH but above the photosphere of the disk-driven wind. The phenomenology of S26 is essentially similar to that of an ultraluminous X-ray source (ULX), but seen from a large inclination angle: if the viewing angle were close to 0°, the observer would probably be able to see the innermost disk and thus the thermal hard X-ray component with high luminosity, enhanced by the geometric beaming (King et al. 2001; King 2009, 2010; Kaaret et al. 2017; Fabrika et al. 2021; Abaroa et al. 2023; Combi et al. 2024; Veledina et al. 2024).

We present here a study of the thermal and nonthermal radiation of S26 and we explore the relationship between the jet and the disk in this source. We analyze multi-epoch X-ray observations of the system obtained with XMM-Newton and develop a jet model to explain the emission. We study particle acceleration and radiative processes that occur in adiabatic shocks generated at the base of the jet and in the lobes. We also examine the feasibility of MQs as potential PeVatron sources, exploring their ability to accelerate protons to energies of about 1 PeV or higher.

We first describe the observations made by the XMM-Newton space observatory and the results of our data analysis. Then, we describe in Sect. 3 our model for the accretion disk, its wind, and the jets. We detail all the thermal and nonthermal radiative processes of the system in Sect. 4. We present the results of our calculations in Sect. 5 and, after a discussion (Sect. 6), we conclude with a summary and our conclusions.

2. X-ray observations

2.1. Data reduction

The XMM-Newton Observatory (Jansen et al. 2001) consists of three X-ray telescopes, each equipped with a European Photon Imaging Camera (EPIC) at its focus. The EPIC instrument includes the PN detector (Strüder et al. 2001) and two MOS detectors (Turner et al. 2001) covering the 0.15–15 keV energy range with a spectral resolution of R = E/dE ≈ 20–50.

S26 was within the field of view of XMM-Newton on 22 occasions between May 2012 and December 2023, located approximately 3–4 arcmin from the center point (on axis). Since S26 is an extragalactic source, to maximize the number of photons and thus the signal-to-noise ratio, we specifically selected long-duration observations, excluding those with raw exposure times less than 40 ks from the initial dataset. Therefore, we focused our analysis on 11 observations meeting this criterion.¹

We reduced and analyzed the XMM-Newton data using the Science Analysis System (SAS) software package version 20.0.0 and the High Energy Astrophysics Software (HEASOFT). We calibrated, filtered, and cleaned the event lists by applying the epproc and emproc tasks, for PN and MOS respectively, using the most up-to-date calibration files. We then identified and removed periods dominated by high particle background, defined as count rate greater than 0.35 cts s⁻¹ for MOS and count rate greater than 0.4 cts s⁻¹ for PN.

To extract the source spectrum, we selected a circular region of 15 arcsec radius centered at α = 23^h58^m00^s, δ = −32°33′21″ (J2000; Dopita et al. 2012). To account for the background emission, photons were selected from a source-free 30 arcsec circular region on the same chip as the source. We used the evselect task to retain source events for spectra with FLAG == 0 and single and double events with PATTERN < = 4 for the PN camera, and PATTERN <  = 12 for the MOS cameras. We generated the response matrix file (rmf) and ancillary response file (arf) using the rmfgen and arfgen tasks, respectively.

Finally, we further refined our dataset by excluding from our analysis the observations that were dominated by background above 2 keV, as this provided insufficient data for a robust spectral analysis of the source. We then analyzed three observations with obsIDs 0748390901 (October 12, 2014), 0804670301 (May 5, 2017), and 0853981001 (November 22, 2019), all performed using the EPIC full frame mode. The average durations for these observations were 46.1, 46.9, and 37.3 ks, respectively. The list of observations used in this study, along with their corresponding net exposure times after filtering, is presented in Table 1.

2.2. Spectral analysis

We performed the X-ray spectral analysis using the XSPEC software v.20.0.0 (Arnaud 1996) in the 0.5–5.0 keV energy range for obsIDs 0748390901 and 0804670301, and 0.5–4.0 keV for obsID 0853981001. We attempted different binning for the spectra using the specgroup tool. We applied a minimum of 16 counts per bin for obsIDs 0748390901 and 0853981001, and 20 counts per bin for obsID 0804670301 to allow for Gaussian statistics when fitting, and with an oversampling factor of 3. We began by fitting the spectra of each observation with two different single-component models. First, we fit the three XMM-Newton/EPIC spectra simultaneously with both an absorbed power-law (const*tbabs*powerlaw in XSPEC) and an absorbed blackbody model (const*tbabs*bbody in XSPEC).

Because of the lack of sufficient detected photons, we fixed the value of the hydrogen column density in all cases to N_H = 1.2 × 10²⁰ cm⁻² (Kalberla et al. 2005). We adopted the Wilms et al. (2000) abundance tables and cross sections from Verner et al. (1996). We fixed the constant factor to unity for the PN camera, but left it free to vary for the two MOS cameras to account for uncertainties in the cross-calibration. We also constrained the model parameters between the three detectors. Neither model provided a formally acceptable fit; the blackbody model did not adequately account for the observed emission above 2 keV, while the power-law model resulted in a poorly constrained spectral index.

We then applied a model consisting of a power-law and a blackbody component affected by photoelectric absorption (const*tbabs*(powerlaw + bbody) in XSPEC). We again fixed the model parameters, except for N_H, which was frozen, and the constant factor, which was fixed at one for PN but allowed to vary freely for the MOS cameras. This model gave a good fit to the observed emission, with χ² = 16.44/16.38/16.45 for 19/18/13 degrees of freedom (dof).

Overall, the parameters obtained with the PN and MOS cameras are consistent with each other within errors. The best-fitting parameters are reported in Table 2 (the errors indicate 1σ confidence levels), while the corresponding spectra and fitted model are shown in Fig. 1.

Fig. 1. X-ray multi-epoch observations. XSPEC CONSTANT*TBABS*(BBODY+POWERLAW) fitting model of the spectrum of S26 using observations 0748390901 (top), 0804670301 (middle) and 0853981001 (bottom) taken by XMM-Newton. The black, red, and green colors represent the PN, MOS1, and MOS2 cameras, respectively.

The model fit yields a power-law photon index Γ and a blackbody temperature kT_bb which are consistent with the values reported in previous X-ray studies (Soria et al. 2010). We obtain a 0.5–5.0 keV unabsorbed flux of $1.8_{-0.3}^{+0.2}\times {10}^{-14}$ erg cm⁻² s⁻¹ and $2.8_{-0.4}^{+0.3}\times {10}^{-14}$ erg cm⁻² s⁻¹ for obsIDs 0748390901 and 0804670301, respectively, and for obsID 0853981001, a 0.5–4.0 keV unabsorbed flux of $1.7_{-0.7}^{+0.1}\times {10}^{-14}$ erg cm⁻² s⁻¹. Assuming a distance of d = 3.9 Mpc (Karachentsev et al. 2003) these fluxes correspond to an unabsorbed luminosity of the order of ∼10³⁷ erg s⁻¹. The X-ray spectra of S26 obtained with XMM-Newton cannot be spatially resolved due to the limited resolution of the instrument. Therefore, we used the smallest possible extraction region to obtain the spectrum that includes the entire source. Thus, our model suggests that the emission from S26 could have two possible origins, one thermal and one nonthermal. This is consistent with the results reported by Soria et al. (2010), who show that the X-ray spectrum of S26 can be modeled by a combination of a thermal component and a nonthermal power law, the former associated with the lobes and the latter associated with the emission from the core. These observations were made with the Chandra space observatory, which has a better spatial resolution.

The temperature of the thermal component has negligible errors and remains essentially constant on timescales of years between observations. Conversely, the photon index of the nonthermal component has such large errors that we cannot rule out variability in the X-ray emission from the core.

3. Physical model of S26

We assume that the X-ray binary consists of a BH with mass M_BH = 10 M_⊙ and an early-type donor star. The star overflows its Roche lobe, transfers mass to the BH through the Lagrange point, and an accretion disk is formed due to the angular momentum of the system. In the following, we describe the semi-analytical models we use to study the accretion disk of the BH, the wind ejected from its surface, and the jet. We show a schematic picture of S26 according to our model in Fig. 2 and list all the parameters in Table 3.

Fig. 2. Conceptual scheme of the MQ S26, not to scale. The stellar-mass BH and the Wolf-Rayet star are in the core. The extended lobes are at the ends of the two opposite jets, where shocks form at the hotspots. There, the backward shock accelerates particles to relativistic energies, producing synchrotron radio emission, and the forward radiative shock produces optically thin thermal X-ray emission. Nonthermal X-ray emission is produced by the electron-synchrotron mechanism at the base of the jet in the core, close to the BH but above the photosphere of the disk-driven wind.

3.1. Accretion disk

We assume that the BH accretes matter at a super-Eddington rate at the outer part of the disk, $\dot{m}={\dot{M}}_{\mathrm{input}}/{\dot{M}}_{\mathrm{Edd}}\gg 1$, where Ṁ_input is the input of mass per unit of time. The Eddington rate is given by

$$\begin{array}{c}\hfill M_{\mathrm{Edd}}=\frac{L_{\mathrm{Edd}}}{\eta c^2}\approx 2.2\times {10}^{-8}{M}_{\mathrm{BH}}{\mathrm{yr}}^{-1}=1.4\times {10}^{18}\frac{M}{M_{\odot }}\mathrm{g}{\mathrm{s}}^{-1},\end{array}$$(1)

with L_Edd the Eddington luminosity (defined as the luminosity required to balance the attractive gravitational pull of the BH by radiation pressure). The parameter η ≈ 0.1 is the accretion efficiency (Fukue 2004), and c is, as usual, the speed of light.

The critical radius, given by $r_{\mathrm{crit}}\sim 40\dot{m}{r}_{\mathrm{g}}$, separates the disk into two regions: a standard outer disk (Shakura & Sunyaev 1973) and a radiation-dominated inner disk with advection (Fukue 2004), where r_g = GM_BH/c² is the gravitational radius of the BH, with G the gravitational constant. The disk becomes geometrically thick in the inner region, where the ejection of winds by the radiation force helps to regulate the mass-accretion rate onto the BH (Ṁ_acc) at approximately the Eddington rate (Fukue 2004):

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{acc}}(r_{\mathrm{d}})=\{\begin{array}{c}{\dot{M}}_{\mathrm{input}},r_{\mathrm{d}}>r_{\mathrm{crit}}\hfill \\ \hfill \\ {\dot{M}}_{\mathrm{input}}\frac{r_{\mathrm{d}}}{r_{\mathrm{crit}}}=\frac{1}{40}\frac{r_{\mathrm{d}}}{r_{\mathrm{g}}}{\dot{M}}_{\mathrm{Edd}},r_{\mathrm{d}}\le r_{\mathrm{crit}},\hfill \end{array}\end{array}$$(2)

where r_d is the distance to the BH in the equatorial plane of the disk. Such regulation leads to a total mass loss rate in the winds Ṁ_w that is approximately equal to the accretion input, ${\dot{M}}_{\mathrm{w}}\approx {\dot{M}}_{\mathrm{input}}\gg {\dot{M}}_{\mathrm{Edd}}$. The accretion power at a radius r_d is given by $L_{\mathrm{acc}}(r_{\mathrm{d}})\approx {\dot{M}}_{\mathrm{acc}}(r_{\mathrm{d}})c^2$.

In this supercritical scenario, the X-ray emitting region of the disk is completely obscured by the opaque wind described below. The radiation from the disk is therefore absorbed and reprocessed by the wind and shifted to lower energies.

3.2. Absorption by the disk-driven wind

We assume that the wind ejected from the BH disk is symmetric, smooth, and expands at a constant velocity, estimated to be ${\beta }_{\mathrm{w}}=v_{\mathrm{w}}/c=1/\sqrt{40\dot{m}}$ (King 2010). Since β_w ≪ 1 for ṁ ≫ 1, the Lorentz factor of the wind is γ_w ≈ 1; hence, there is no dependence of the wind properties with the line of sight (Θ) due to relativistic effects, so 1 − β_w cosΘ ≈ 1 (see e.g., Abaroa & Romero 2024b). However, we note that there could be a dependence on the viewing angle because of the narrow conical low-density funnel (with semi-opening angles ≲15° in some cases, as recently shown by Veledina et al. 2024, for Cygnus X-3) that forms due to the magneto-centrifugal barrier around the z-axis above the BH and causes the geometric beaming (King et al. 2001; King 2009, 2010; Kaaret et al. 2017; Fabrika et al. 2021; Abaroa et al. 2023). A low tilt with respect to the line of sight would allow the hard X-ray thermal emission to escape from the inner accretion disk and reach the observer, but this is not the case for S26 where a high inclination is inferred from radio and X-ray observations (Pakull et al. 2010; Soria et al. 2010).

The density of this wind at a distance R from the BH is

$$\begin{array}{c}\hfill {\rho }_{\mathrm{w}}(r)={\dot{M}}_{\mathrm{w}}/4\pi R^2{v}_{\mathrm{w}}.\end{array}$$(3)

where $R=\sqrt{r^2+z^2}$, with r and z cylindrical coordinates. We assume a steady state so that v_w and Ṁ_w are constants and therefore the density depends only on the distance to the BH.

The wind itself is opaque until it reaches the photospheric radius, where it becomes transparent to its own radiation (Poutanen et al. 2007; Fukue 2009; Tomida et al. 2015; Zhou et al. 2019). The apparent photosphere is defined as the surface where the optical depth τ_photo is unity for an observer at infinity. Its height z_photo, measured from the equatorial plane, can be found by integrating over the wind density (Fukue 2009):

$$\begin{array}{c}\hfill {\tau }_{\mathrm{photo}}={\int }_{z_{\mathrm{photo}}}^{\infty }{\gamma }_{\mathrm{w}}(1-{\beta }_{\mathrm{w}}cos\mathrm{\Theta }){\kappa }_{\mathrm{co}}{\rho }_{\mathrm{co}}\mathrm{d}z\approx {\int }_{z_{\mathrm{photo}}}^{\infty }{\kappa }_{\mathrm{co}}{\rho }_{\mathrm{co}}\mathrm{d}z=1,\end{array}$$(4)

where κ_co ≈ κ is the opacity and ρ_co ≈ ρ_w (the sub-index co refers to the comoving frame). As we assume a fully ionized wind, the opacity is dominated by free electron scattering (κ = σ_T/m_p, with σ_T the Thomson scattering and m_p the proton mass). We note that the location of the wind photosphere is completely determined by the rate of mass loss, and thus by the accretion input of matter.

We assume an accretion rate input of ${\dot{M}}_{\mathrm{input}}=10{\dot{M}}_{\mathrm{Edd}}$, which leads to a photosphere with a height of z_photo ≈ 2.6 × 10⁹ cm above the BH. Fig. 3 shows the geometry of the photosphere in Cartesian coordinates, in units of the gravitational radius. The color map represents the temperature, which is ∼3.9 × 10⁵ K at the top. In the present model for high inclination angles, the photosphere completely obscures the surroundings of the BH, and thus the X-ray emission from the inner disk cannot reach the observer². We estimate the temperature and emission from the wind photosphere in Sect. 4.1, along with all the other radiative processes.

Fig. 3. Photosphere of the disk-driven wind for an accretion rate of 10 ṀEdd. The axes are in Cartesian coordinates and units of the gravitational radius. The height zphoto is 1777 rg ( ≈ 2.6 × 109 cm) above the BH and its top temperature is ∼3.9 × 105 K.

3.3. Jet

Two highly collimated jets are generated near the BH in opposite directions, expanding laterally at the speed of sound. The jets are perpendicular to the orbital plane of the binary system. In the following sections, we discuss the characterization of jet power and launch, the magnetic field and matter density, as well as the acceleration of particles at the base and in the terminal regions of the jet.

3.3.1. Jet: available power and launching

Using general relativistic magnetohydrodynamics simulations, Liska et al. (2018) shows that, under certain physical conditions, a large-scale poloidal magnetic field can be generated from a purely toroidal field in the disk. The poloidal magnetic flux accumulates around the BH until it becomes dynamically important, leading to a magnetically saturated – also called “arrested” in the literature – accretion disk (Narayan et al. 2003). This happens at a certain radius r_d ≤ r_ms.

According to the disk-jet coupling hypothesis of Falcke & Biermann (1995), the kinetic power of each jet is assumed to be proportional to the accretion power at the magnetically saturated radius (where the poloidal magnetic flux becomes dynamically important, see Begelman et al. 2022), L_j ∼ q_j L_acc(r_ms), where q_j is the efficiency in transferring power from the disk to the jet. We assume that the disk is saturated at a distance from the BH of r_ms = 100 r_g. At this radius, a BH of 10M_⊙ accretes matter at ${\dot{M}}_{\mathrm{acc}}(r_{\mathrm{ms}})=2.5{\dot{M}}_{\mathrm{Edd}}=3.5\times {10}^{19}\mathrm{g}{\mathrm{s}}^{-1}$, so the accretion power is L_acc(r_ms) = 3.15 × 10⁴⁰ erg s⁻¹. This luminosity is not sufficient to explain the jet power in S26: some energy should also be extracted from the ergosphere of the rotating BH, which spins down as the rotational energy is transferred to the jet (Blandford & Znajek 1977). In such cases, the efficiency of converting accretion power into jet power can exceed 100 percent in highly rotating BHs, reaching transient values of 300 percent in the most extreme cases (Tchekhovskoy et al. 2010, 2011; McKinney et al. 2012). Considering that the jet kinetic luminosity of S26 estimated by Pakull et al. (2010) is L_j = 5 × 10⁴⁰ erg s⁻¹, then the jet-disk coupling constant should be q_j = L_j/L_acc ≈ 1.60.

We assume that the jet is fully formed at a distance z₀ = 100 r_g from the BH and with a radius r₀ = 0.1 z₀ = 10 r_g. Afterward, the outflow expands as a cone of radius r_j(z) = r₀(z/z₀) = 0.1 z (that is, with a typical semi-opening angle of θ_j ≈ 5.7°). We adopt a moderate value for the Lorentz factor γ_j = (1 − v_j²/c²)^−1/2 = 3, so v_j = 0.94 c.

Although the jet starts as a Poynting flux outflow, a fraction q_rel ≈ 0.01 of the jet power is expected to be in the form of relativistic particles at distances greater than ∼1000 r_g (see Romero & Gutiérrez 2020, for the matter loading mechanisms). Then we can write L_rel = q_rel L_j. We include both the hadronic and leptonic fractions, L_rel = L_p + L_e. The energy distribution between hadrons and leptons is unknown; we assume an equipartition near the base of the jet (L_p =  L_e), and a hadronic-dominant scenario in the terminal regions (L_p = 100 L_e), since strong matter entrainment along the jet path is expected at distances of ∼175 pc from the BH (e.g., Bosch-Ramon et al. 2005a; Bordas et al. 2009). At such distances the lobes formed by the jet are exposed to the interstellar medium (ISM) with conditions similar to those found in supernova remnants, where the hadron-to-lepton ratio is 100, in agreement with rough estimates of cosmic ray (CR) composition determined in the vicinity of the Earth (e.g., Beresinsky & Zatsepin 1969).

3.3.2. Jet: magnetic field and matter density

The magnetic field B₀ near the BH can be estimated from the Blandford-Znajek power P_BZ of the jet (Blandford & Znajek 1977; Romero & Gutiérrez 2020):

$$\begin{array}{c}\hfill P_{\mathrm{BZ}}\sim {10}^{46}{\left(\frac{B_0}{{10}^4\mathrm{G}}\right)}^2{\left(\frac{M_{\mathrm{BH}}}{{10}^9{\mathrm{M}}_{\odot }}\right)}^2{a}^2\mathrm{erg}{\mathrm{s}}^{-1},\end{array}$$(5)

where 0 ≤ a < 1 is the normalized spin of the BH. This is a 1st-order approximation, valid for low spin values: if a > 0.5, then the above equation underestimates the jet power by a factor of ∼3. As we are dealing with a highly rotating BH to account for the observed jet power in S26, we apply this factor (L_j ∼ 3 P_BZ) to correct the deviation from the more accurate 6th-order formula derived by Tchekhovskoy et al. (2010, 2011) (see also Chapter 3 of Contopoulus et al. 2015, for details). We obtain a magnetic field of B₀ ≈ 1.61 × 10⁹ G near the BH.

The magnetic field in the jet decreases with distance z from the BH: if the flow expands adiabatically,

$$\begin{array}{c}\hfill B(z)=B_0{\left(\frac{z^{\prime }}{z}\right)}^m,\end{array}$$(6)

where z′≈10 r_g and 1 ≤ m ≤ 2 is the magnetic index, which depends on the geometry of B. We assumed a mixed geometry of the magnetic field – poloidal and toroidal – beyond z′, so m = 1.5. The magnetic energy density is given by

$$\begin{array}{c}\hfill e_{\mathrm{m}}(z)=\frac{B^2(z)}{8\pi }.\end{array}$$(7)

The force exerted by the toroidal field helps to collimate the jet, where collimation is the degree of parallelism of the streamlines. A cylindrical jet with all streamlines parallel to the symmetry axis is then perfectly collimated (Romero & Vila 2014). We note that collimation is different from confinement, which refers to the width of the outflow³.

The number density of cold protons at a distance z from the BH in the jet is (e.g., Sotomayor Checa & Romero 2019)

$$\begin{array}{c}\hfill n_{\mathrm{p}}(z)\sim \frac{{\dot{M}}_{\mathrm{j}}}{\pi r_{\mathrm{j}}^2(z)m_{\mathrm{p}}{v}_{\mathrm{j}}},\end{array}$$(8)

where the mass loss rate in the jet is given by ${\dot{M}}_{\mathrm{j}}=L_{\mathrm{j}}/({\gamma }_{\mathrm{j}}-1)c^2$. The energy density of matter reads

$$\begin{array}{c}\hfill e_{\mathrm{p}}(z)=\frac{{\dot{M}}_{\mathrm{j}}}{2\pi z^2}{v}_{\mathrm{j}}.\end{array}$$(9)

3.3.3. Jet: shocks

The X-ray and radio maps of S26 show nonthermal emission from the core and the lobes (Pakull et al. 2010; Soria et al. 2010). Therefore, we considered two acceleration regions for relativistic particles in the jet: one at the base of the jet, but above the photosphere of the disk-driven wind, and the other at the terminal region of the jet. The acceleration is assumed to be due to diffusive shock acceleration in the matter-dominated regions of the system. We neglected further shocks or other acceleration mechanisms (such as magnetic reconnection) that may occur, since the emission regions are located in shock-prone regions in the observations.

Internal shocks can be produced near the base of the jet (Kaiser et al. 2000). These shocks occur when the jet is matter-dominated, that is, when the energy density of the matter is greater than the magnetic energy density, e_p(z) > e_m(z). By comparing Eqs. (7) and (9), we find that this happens at a height z_b ∼ 5 × 10⁹ cm, where we assume that particle acceleration occurs. We note that this shock region is located above the wind photosphere derived in the previous section (z_b > z_photo).

At the head of the jet, two more shocks are generated: a forward radiative shock that propagates into the ISM, and a reverse adiabatic shock that propagates into the jet. Particles escaping from the reverse shock form the thermal cocoon of the source. Particle acceleration to relativistic energies occurs only in the reverse adiabatic shock. Conversely, in the radiative shocks, the gas emits large amounts of thermal, optically thin radiation, and no particle acceleration is produced. The calculations of the radiation and its partial absorption are detailed in the next section.

4. Radiative processes

In the following, we characterize the different radiative processes that we have calculated for S26. First, we describe the thermal emission (from the disk-driven wind and the bow shocks of the jets), and then the nonthermal radiation (inner jet and reverse shocks in the terminal region). The processes studied include the acceleration, escape, cooling, distribution, and radiation of different types of particles. At the end of the section, the absorption of the radiation in the source is described.

4.1. Wind emission

We calculate the thermal emission from the photosphere of the disk-driven wind assuming that it expands at a constant rate equal to its terminal velocity. Since the mass loss rate of the disk is much higher than the critical rate, the wind is optically thick, and the emission from the disk is thermalized and re-emitted as a blackbody spectrum at its photosphere (Zhou et al. 2019). We note, however, that in the central region, as we already have mentioned, the dense wind surrounds a low-density funnel, where the hard X-rays from the central accretion disk can escape to reach a face-on observer (see Sect. 3.2) (King 2010; Zhou et al. 2019).

The photosphere temperature measured by an observer at infinity is given by (Fukue 2009):

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{T}}{T}_{\mathrm{w}}^4=\frac{\dot{e}{L}_{\mathrm{Edd}}}{{(1-\beta cos\mathrm{\Theta })}^{4}4\pi R^2}\approx \frac{\dot{e}{L}_{\mathrm{Edd}}}{4\pi R^2},\end{array}$$(10)

where $R=\sqrt{r^2+z^2}$. The radiative power of the wind in supercritical sources is limited to the Eddington luminosity (Zhou et al. 2019), so we assume that the comoving luminosity of the wind is a fraction ė = 0.1 of L_Edd (Sotomayor & Romero 2022). The temperature of the wind at the photospheric height is then T_photo = 3.9 × 10⁵ K.

4.2. Thermal lobes

The radiative forward shocks generated at the jets’ terminal regions, as they interact with the ISM, produce thermal extended emission due to the free-free interaction of the particles (thermal Bremsstrahlung). The emissivity of a thermal, optically thin plasma at a temperature T for a frequency ν is given by (Longair 2011):

$$\begin{array}{c}\hfill {\kappa }_{\nu }=\mathbb{C}{Z}^2{T}^{-1/2}{N}_{\mathrm{i}}{N}_{\mathrm{e}}g(\nu ,T)exp\left(\frac{-h\nu }{k_{\mathrm{B}}T}\right)\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-3}{\mathrm{Hz}}^{-1},\end{array}$$(11)

where ℂ = 6.8 × 10⁻³⁸, Z is the atomic number, N_i and N_e are the number density of ions and electrons (in units of cm⁻³), h is the Planck constant, and k_B the Boltzmann constant. The Gaunt factor g(ν, T), for X-ray energies, is (Longair 2011):

$$\begin{array}{c}\hfill g(\nu ,T)\approx \{\begin{array}{c}\frac{\sqrt{3}}{\pi }ln\left(\frac{k_{\mathrm{B}}T}{h\nu }\right),h\nu /k_{\mathrm{B}}T\ll 1\hfill \\ \hfill \\ \sqrt{\frac{h\nu }{k_{\mathrm{B}}T}},h\nu /k_{\mathrm{B}}T\gg 1.\hfill \end{array}\end{array}$$(12)

The luminosity of a region of volume V as a function of the frequency is then L(ν) = V κ_ν ν. This thermal radiation from the lobes of S26 is identified with the observed soft X-ray photons detected by the Chandra and XMM-Newton telescopes.

4.3. Acceleration and escape timescales

The diffusive acceleration rate of the particles is given by (e.g., Protheroe 1999):

$$\begin{array}{c}\hfill t_{\mathrm{acc}}^{-1}={\eta }_{\mathrm{acc}}\frac{eZcB}{E},\end{array}$$(13)

where e is the electric charge, Z the atomic number, c the speed of light, B the magnetic field, and E is the energy of the particle. The acceleration efficiency η_acc of the process can be estimated with the velocity of the jet, assuming Bohm diffusion: η_acc = 3β_sh²/8, where β_sh = v_sh/c, with v_sh the shock velocity in the corresponding region: in the case of the base, v_sh ≡ v_j while in the reverse shock v_sh ≡ v_rs (see Table 3). These values lead to η_acc ≈ 0.1.

Escape is advective at the base of the jet (particles are removed from the acceleration region by the bulk motion of the jet) and diffusive in the lobes:

$$\begin{array}{c}\hfill \begin{array}{c}{t}_{\mathrm{adv}}^{-1}={\left(\frac{\mathrm{\Delta }{x}_{\mathrm{acc}}}{v_{\mathrm{j}}}\right)}^{-1}(\mathrm{advection})\hfill \\ \hfill \\ t_{\mathrm{diff}}^{-1}={\left(\frac{\mathrm{\Delta }{x}_{\mathrm{acc}}^2}{D(E)}\right)}^{-1}(\mathrm{diffusion}),\hfill \end{array}\end{array}$$(14)

where Δx_acc is the size of the corresponding acceleration region and D(E) is the diffusion coefficient, considered in the Bohm regime, so

$$\begin{array}{c}\hfill D(E)=D_{\mathrm{B}}(E)=\frac{Ec}{3Be}.\end{array}$$(15)

Accelerated particles can cool by several different processes, which are described in the next subsection. The maximum energy for each type of particle can be deduced by looking at the point where the acceleration rate equals the total cooling or escape rate. This energy cannot exceed the maximum energy locally imposed by the Hillas criterion, $E_{\mathrm{e},\mathrm{p}}^{\max }<E_{\mathrm{Hillas}}^{\max }=eZ\mathrm{\Delta }{x}_{\mathrm{acc}}B$.

4.4. Energy losses

The timescales associated with cooling are related to the total energy loss of the particles:

$$\begin{array}{c}\hfill \frac{\mathrm{d}E}{\mathrm{d}t}\approx \frac{-E}{t_{\mathrm{cool}}},\end{array}$$(16)

where the total cooling rate is

$$\begin{array}{c}\hfill t_{\mathrm{cool}}^{-1}=\sum _i{t}_i^{-1},\end{array}$$(17)

where t_i is each timescale of the cooling processes involved.

Adiabatic losses may also be relevant at the internal shock, because of the pressure exerted by the particles to expand the lateral walls of the jet,

$$\begin{array}{c}\hfill t_{\mathrm{ad}}^{-1}={\left(\frac{3\mathrm{\Delta }{x}_{\mathrm{acc}}}{2v_{\mathrm{j}}}\right)}^{-1}.\end{array}$$(18)

Radiative cooling is caused by various nonthermal processes as a consequence of the interaction of particles with radiation, magnetic field, and matter. The intensity of the magnetic field in the shocks at the base of the jet is calculated according to Eq. (6); at the reverse shock it is fixed to a typical value of the ISM as 2.5 μG. In addition, we considered an enhancement of ∼4 due to compression (Lamers & Cassinelli 1999). The dominant photon fields are provided by the photosphere of the disk-driven wind (for the base of the jet), which decays with 1/d², and the cosmic microwave background (for the terminal region). The cold matter density is calculated according to Eq. (8) for the base of the jet, and is set to 0.1 cm⁻³ for the terminal region (Soria et al. 2010).

Our model is lepto-hadronic, and so we calculate the following cooling processes numerically:

– Synchrotron: Interaction of protons and electrons with the ambient magnetic field.

– Inverse Compton: Collision of relativistic electrons with photons of the ambient radiation field.

– Relativistic Bremsstrahlung: Coulombian interactions between relativistic electrons and cold matter.

– Photo-hadronic interactions: Interaction of highly relativistic protons with photons of the ambient radiation field. This interaction produces pions, which decay into gamma rays, e⁺e⁻ pairs, and neutrinos.

– Proton-proton: Collision of relativistic protons with cold matter. This process is also a source of gamma rays, pairs, and neutrinos.

4.5. Particle distributions

The relativistic particles have a distribution given by dN = n(r,E,t)dEdV, where n is the number density of particles, t the time, r the position, V the volume, and E the energy. The evolution of this distribution is determined by the transport equation (see Ginzburg & Syrovatskii 1964). We solve this equation numerically in steady state and in the one-zone approximation:

$$\begin{array}{c}\hfill \frac{\partial }{\partial E}\left[\frac{\mathrm{d}E}{\mathrm{d}t}N(E)\right]+\frac{N(E)}{t_{\mathrm{esc}}}=Q(E),\end{array}$$(19)

where t_esc is the corresponding escape time, and the particle injection function,

$$\begin{array}{c}\hfill Q(E)=Q_0{E}^{-p}exp(-E/E_{\max }),\end{array}$$(20)

is a power-law in the energy with an exponential cutoff and a spectral index of p ∼ 2, which is characteristic of the Fermi first-order acceleration mechanism (see e.g., Drury 1983). The normalization constant Q₀ is obtained from the available power for the different particle species

$$\begin{array}{c}\hfill L_{\mathrm{rel}(\mathrm{e},\mathrm{p})}=V_{\mathrm{acc}}{\int }_{E_{(\mathrm{e},\mathrm{p})}^{\min }}^{E_{(\mathrm{e},\mathrm{p})}^{\max }}{E}_{(\mathrm{e},\mathrm{p})}{Q}_{(\mathrm{e},\mathrm{p})}(E_{(\mathrm{e},\mathrm{p})})\mathrm{d}{E}_{(\mathrm{e},\mathrm{p})},\end{array}$$(21)

where V_acc is the volume of the acceleration region, and $E_{(\mathrm{e},\mathrm{p})}^{\max }$ the maximum energy reached by protons and electrons.

4.6. Nonthermal emission

Once we have the particle distributions, we calculate the spectral energy distribution (SED) for each of the relevant processes involved in the cooling. This radiation corresponds in particular to the nonthermal emission of S26 detected from the core (X-rays) and the lobes (radio emission). The reader is referred to Romero & Vila (2008), Romero & Paredes (2011), Romero & Vila (2014), and references therein for additional details on the radiative processes.

Relativistic Doppler affects the emission produced at the internal shock due to the relativistic motion of the gas. The Doppler factor depends on the angle i of the jets with respect to the line of sight:

$$\begin{array}{c}\hfill D=\frac{1}{{\gamma }_{\mathrm{j}}(1\pm {\beta }_{\mathrm{j}}cosi)},\end{array}$$(22)

where the sign ± accounts for the jet ( − ) and the counterjet ( + ). The apparent luminosity is then D² times the comoving luminosity,

$$\begin{array}{c}\hfill L_{\mathrm{app}}=D^2{L}_{\mathrm{emitted}}.\end{array}$$(23)

We get D² ≈ 0.19 for the jet and ≈0.07 for the counterjet (cj), so the apparent luminosity is less than the emitted. Since L_cj ∼ 0.4 L_j, then L_j + cj ∼ 1.4 L_j. We note also that since the lobes are static, their emission is not affected by relativistic effects (D ≡ 1).

4.7. Absorption

Finally, we calculate the absorption of gamma rays by pair creation from photon-photon annihilation, γ + γ → e⁺ + e⁻. The nonthermal photons on their way out of the emission region can find thermal ambient photons and annihilate. High-energy radiation emitted in the jets is also attenuated by γγ–annihilation with synchrotron photons produced internally. The absorption is quantified by the optical depth of the medium, τ_γγ, which depends on the energy. The attenuated, observed, luminosity reads:

$$\begin{array}{c}\hfill L_{\mathrm{obs}}(E_{\gamma })=L_{\mathrm{app}}(E_{\gamma })·{\mathrm{e}}^{-\tau },\end{array}$$(24)

where e^−τ is the attenuation factor. We refer to Sect. 4.5 of Abaroa et al. (2023) for details in the calculation of the absorption.

5. Results

5.1. Energy gain and losses

Figure 4 shows the timescales associated with particle acceleration, escape, and cooling. We calculate these processes for the base of the jet (left) and the reverse shocks in the lobes (right), in each case for both leptons (top) and hadrons (bottom). The efficiency of the acceleration in all cases is η_acc ∼ 0.1.

Fig. 4. Cooling timescales of particles accelerated at the base of the jet (left) and at the terminal region (right). We show the energy gains and losses for electrons in the top panel and the same for protons in the bottom panel. In the internal shock at the base of the jet, the particles cool mainly by the synchrotron mechanism due to the strong magnetic field, reaching energies of 10 GeV (electrons) and 50 PeV (protons). At the reverse shock in the terminal region of the jet, electrons also cool mainly by synchrotron radiation, reaching energies of 5 PeV, while protons reach energies of ∼10 PeV. The efficiency of the acceleration in all cases is ηacc ∼ 0.1.

In the internal shock at the base of the jet, electrons cool mainly by synchrotron (as expected, because of the strong magnetic field), reaching a maximum energy of ∼10 GeV, whereas protons reach a maximum energy of ∼50 PeV, being synchrotron and advective escape the dominant timescales for them. The Hillas maximum allowed energy is $E_{\mathrm{Hillas}}^{\max }=2\times {10}^{17}\mathrm{eV}$ in this region.

Electrons accelerated in the reverse shock in the lobes also cool mainly by synchrotron radiation, reaching a maximum energy of ∼5 PeV. Protons are accelerated to energies of ∼10 PeV, and their dominant loss is diffusive escape. These ultra-relativistic protons, once they have escaped, will in turn propagate through the ISM and contribute to the CR population in the host galaxy of the MQ.

5.2. Nonthermal radiation

Figure 5 shows the nonthermal SEDs of the jet base (left) and the lobes (right). In the case of the base, the electron-synchrotron emission extends from 10⁻¹ to ∼10⁶ eV with an almost constant luminosity of 5 × 10³⁶ erg s⁻¹, while the proton-synchrotron mechanism provides emission up to ∼1 TeV, although internal photon-photon annihilation produces absorption for energies higher than ∼1 GeV. The beaming effects are already included. The electron-synchrotron mechanism is then responsible for the X-ray emission produced at the core, detected by Chandra and XMM-Newton. For the lobes, synchrotron, Bremsstrahlung, inverse Compton, and proton-proton interactions are all relevant. The luminosity produced by synchrotron is above 10³⁴ erg s⁻¹ from very low (10⁻⁹ eV) to high (10⁷ eV) energies, while the emission due to π⁰ decays through the pp channel is above 10³⁵ erg s⁻¹ in the 10⁸ − 10¹⁴ eV energy range. The electron-synchrotron mechanism is responsible for the radio emission in the hotspots detected with ATCA by Soria et al. (2010), and pp and the subsequent decays are responsible for the high-energy gamma-ray emission.

Fig. 5. Nonthermal emission from the jet. Left: Nonthermal SEDs of radiation produced at the internal shock in the jet, close to the BH but above the disk-driven wind photosphere. The dominant mechanism is electron synchrotron. For energies > 1 GeV the emission is completely attenuated due to internal absorption by pair creation. Dashed lines represent the unabsorbed emission. Right: Nonthermal SEDs of the radiation produced at the reverse shock in the terminal region of the jet. Electron-synchrotron radiation is the dominant process at low energies while proton-proton dominates at high and very high energies. No attenuation affects the radiation produced at the lobes.

5.3. Total SED

Figure 6 shows the total SED of S26, already corrected for absorption and in the observer’s frame (on Earth). We can see the thermal bump caused by the wind photosphere and the forward shocks in the lobes. We also distinguish several plateaux due to the nonthermal processes (the total contribution of the internal shock in the jet and the backward shocks in the lobes). The emission from the photosphere of the disk-driven wind is calculated with an accretion input of 10 Ṁ_Edd. The optically thin emission from the radiative shocks in the terminal region of the jets fits the detected soft X-ray emission with a temperature of ∼4 × 10⁶ K, while the nonthermal emission from the inner jet (plus the small contribution from the counterjet) fits well the hard X-ray emission (Chandra from Pakull et al. 2010, and XMM-Newton from this paper). The combined nonthermal radiation from both lobes provides a good fit to the radio emission detected with ATCA (Pakull et al. 2010; Soria et al. 2010). The gamma rays from the lobes reach energies of ∼10¹⁴ eV. As the plot shows, the source cannot be detected at gamma-ray energies with current instruments, so it should be studied with radio interferometric arrays in the southern hemisphere (such as MeerKAT -originally known as the Karoo Array Telescope- and the future Square Kilometer Array -SKA-) and with space-based X-ray observatories.

Fig. 6. Luminosity of S26 as a function of the energy for our physical model with the parameters listed in Table 3. We plot the thermal and nonthermal SEDs of the system. The thermal contributions are the photosphere of the disk-driven wind and the optically thin emission in the lobes. Nonthermal emission is produced in the internal shocks at the base of the jet (orange solid line, jet + counterjet) and in the reverse shock of the lobes (purple solid line, both lobes). Individual nonthermal processes are shown with dashed and dotted-dashed lines (for details see Fig. 5). We plot radio data from ATCA (red dots, Soria et al. 2010) corresponding to the emission of the two lobes combined, and X-ray data from XMM-Newton (blue dots, this work) corresponding to the core and lobes, for obsID 0748390901. Sensitivities of the Fermi, CTA, SWGO and LHAASO gamma-ray detectors are also shown for two different distances: that from S26 (3.9 Mpc, solid colored lines), and a galactic distance of ∼10 kpc (black dashed lines).

We note, however, that a similar MQ in our Galaxy (i.e., a super-Eddington MQ) could be detected at very high and ultra-high energies with Fermi, Cherenkov Telescope Array (CTA), Southern Wide-field Gamma-ray Observatory (SWGO), and Large High Altitude Air Shower Observation (LHAASO). In fact, this is the case of SS433, a weaker source that has already been detected in very-high-energy (VHE) gamma rays (Abeysekara et al. 2018; H.E.S.S. Collaboration 2024).

5.4. Exploration of the parameter space

We explore the parameter space of S26 to analyze the output spectra as we vary some of the most sensitive parameters. Fig. 7 shows the SEDs we obtain for different scenarios, including the one shown in Fig. 6 (in solid black lines) for comparison. We consider five different scenarios for the jet base (Sb1 − Sb5, dotted lines) and the nonthermal lobe (Sl1 − Sl5, dashed lines):

• Sb1 (cyan line): We only vary the location of the acceleration region, setting z_b = 10¹¹ cm, which results in a lower magnetic field but a more extended emission compared to our original scenario. The nonthermal X-ray emission is underestimated in this case.

• Sb2 (blue line): We vary the BH mass M_BH = 15M_⊙ and the hadron-to-lepton ratio, K_b = 10. The power results underestimated by almost an order of magnitude.

• Sb3 (pink line): We vary the jet Lorentz factor γ_j = 5 and the fraction of energy going to relativistic particles q_rel = 0.1. This produce an overluminous source by a factor of ∼5.

• Sb4 (red line): We vary the BH mass M_BH = 15M_⊙, the spin a = 0.9, and the jet Lorentz factor γ_j = 5. The luminosity results strongly underestimated.

• Sb5 (orange line): We again vary the jet Lorentz factor γ_j = 1.5 and the fraction of the power going to relativistic particles q_rel = 0.1. The result is a powerful source that surpasses the actual observed one by nearly two orders of magnitude.

• Sl1 (cyan line): We vary the magnetic field B_l = 1 μG, the spectral index p_l = 2, and the hadron-to-lepton ratio, K_l = 10. In this scenario the radio emission is overestimated by an order of magnitude.

• Sl2 (blue line): We vary the magnetic field B_l = 1 μG, the number density n_l = 1 cm⁻³, and the hadron-to-lepton ratio, K_l = 10. This scenario overestimates the radio emission by a factor of 2.

• Sl3 (pink line): We vary the spectral index p_l = 2, the number density n_l = 1 cm⁻³, and the hadron-to-lepton ratio, K_l = 10. This is the scenario that overestimates the overall emission the most.

• Sl4 (red line): We vary the hadron-to-lepton ratio K_l = 10 and the number density n_l = 1 cm⁻³. This scenario overestimates the emission by a factor of 2.

• Sl5 (orange line): We vary only the magnetic field, B_l = 1 μG. The emission is underestimated by a factor of ∼3.

Fig. 7. Exploration of the space parameter of the model. We plot SEDs for sets of input parameters different from those of Table 3 and Fig. 6. The dotted lines correspond to the SEDs of the different specific scenarios for the recollimation shock at the base (b) of the jet, (Sb1 − Sb5). The dashed lines correspond to the SEDs of the scenarios of the reverse shock at the lobes (l) of the jet, (Sl1 − Sl5). The corresponding colors for both shocks are cyan, blue, pink, red, and orange for models 1–5, respectively. The different model parameterizations are described in Sect. 5.4. Solid black lines represent our initial scenario. The observational data are plotted in the same way as in Fig. 6.

As can be seen from Fig. 7, no alternative choice of parameter values at the base of the jet fits well with the nonthermal X-ray emission observed from the core of S26. The model shows sensitivity to the BH mass, the hadron-to-lepton ratio, the fraction of power going to relativistic particles, the acceleration region, or the Lorentz factor, all of which we vary to explore the parameter space. The luminosity of the radiation produced in this region varies uniformly over the energy range ∼10⁻¹ − 10⁶ eV in all scenarios. The hadron-to-lepton ratio and the power going to the relativistic particles have a linear behavior with respect to the luminosity: a higher q_rel and a lower K_b leads to a higher synchrotron emission and vice versa. The jet Lorentz factor is one of the most sensitive parameters, since the luminosity in the observer’s frame is corrected by Doppler effects: the higher γ_j, the lower the emission in the direction of Earth, because of the high inclination of the source (see Eq. (22)). The location of the acceleration region, the mass of the BH, and its spin have little effect on these alternative scenarios.

In the case of the lobes, we varied the BH spin and mass, the number density, the magnetic field, and the hadron-to-lepton ratio. Again, the model is sensitive to changes in these parameters, showing in the most extreme case a difference of two orders of magnitude with respect to the observations. Furthermore, scenarios Sl2, Sl3, and Sl4 have nonthermal X-ray luminosities resolvable with Chandra, but this telescope has not detected any nonthermal X-ray emission from the lobes. Unlike the emission from the base of the jet, the lobes do not show uniform variations across the energy spectrum. The curves vary at energies > 10⁸ eV, essentially due to a softening of the particle distribution (scenarios Sl1 and Sl3, where we adopt a spectral index of p = 2), with the emission at 10⁹ eV being an order of magnitude higher than in the case of our initial scenario. The variation of the magnetic field and the hadron-to-lepton ratio does not lead to significant changes at these energies. On the contrary, the model is susceptible to these two parameters in the energy range 10⁻⁸ − 10⁸ eV, where the synchrotron process dominates: a higher B_l and a lower K_l lead to a higher luminosity and vice versa. The BH mass and its spin have little impact on these alternative scenarios of the lobes.

6. Discussion

Our model of S26, which assumes that this source is a super-Eddington accretor seen at a large angle of inclination, allows us to solve the puzzle posed by the large power inferred from the large scale jets and the paucity of X-rays from the central source. It is the dense wind ejected by the accretion disk that obscures most of the X-ray emission. Only nonthermal X-rays produced outside the photosphere of the wind can reach the observer.

Future observations with different detectors could help to further constrain the free parameters of the models, such as the ratio of relativistic protons to electrons, the spectral index of the particles at different distances from the core, the magnetization of the jet, the BH mass and spin, or the accretion rate. For example, the detection of radio emission from the core of the system could be used to better determine the location of the acceleration zone at the base of the jet and the strength of the magnetic field there. This would require interferometric observations with a higher resolution than the one that can be obtained with the ATCA. In addition, the detection of very high energy gamma rays from S26 could indicate a higher particle density at the end of the jets and place a lower limit on the energy of the CRs accelerated in the source. A flaring event caused by the interaction of the jet with overdensities of matter could also favor the detection at very high or ultra-high energies (Araudo et al. 2009; Owocki et al. 2009).

As we mentioned in Sects. 1 and 3, we note that a source similar to S26, but with a different orientation with respect to the line of sight, would be seen as an ULX. The jet emission would also be magnified by a factor of γ_j². This would increase the received flux from the inner jet by about an order of magnitude. For similar considerations about SS433 as a ULX when viewed from a different point of view, see Begelman et al. (2006).

The interaction of the strong wind ejected from the disk of S26 with nearby clouds can also lead to gamma-ray emission, as recently proposed for SS433 (Li et al. 2020). The equatorial outflow could propagate far beyond the system and be revealed if it collides with any clouds (see Picchi et al. 2020, for a characterization of the wind). The shock generated by the collision would convert the kinetic energy of the plasmoids into internal energy and relativistic particles, which could then be cooled by pp interactions with the cloud material. Such a scenario might explain the detection of a GeV source by the Fermi satellite on the side of SS433 (Bordas 2020; Li et al. 2020). Since S26 is at a distance of 3.9 Mpc, a luminosity L_γ > 10³⁹ erg s⁻¹ is required to be detected by Fermi at GeV energies. According to Li et al. (2020), to account for this luminosity via hadronic interactions, we need a total energy of protons interacting with a nearby cloud W_p = 2.5 × 10⁴⁸(L_γ/10³⁴ erg s⁻¹)(n/20 cm⁻³) erg = 2.5 × 10⁵³ erg, where n is the number density of the cloud. The kinetic power of the wind is $L_{\mathrm{w}}={\dot{M}}_{\mathrm{w}}{v}^2/2=1.4\times {10}^{38}\mathrm{erg}{\mathrm{s}}^{-1}$; therefore, the system should be injecting wind into the cloud for a time t_inj = W_p/L_w ∼ 3 × 10⁷ yr, which is much longer than the estimated age of the MQ (t = 2 × 10⁵ yr), so the emission due to wind-cloud interactions would be undetectable with any reasonable choice for the values of the parameters in the model.

Another interesting result of our models is that S26 can accelerate protons to PeV energies both in the inner jet and in the extended lobes. The possibility of powerful MQs producing high-energy CRs has been addressed in several papers (see e.g., Heinz & Sunyaev 2002; Bosch-Ramon et al. 2005a; Fender et al. 2005; Cooper et al. 2020). The protons can escape from the jet through p + p → p + n + π⁺ or p + γ → n + π⁺ channels: the neutrons are no longer confined by the strong magnetic field and can decay far away from the sources, releasing very relativistic protons into the ISM, where they will diffuse until they interact with the surrounding matter and produce gamma rays (Escobar et al. 2022). The source, however, remains mostly obscured by the wind of the accretion disk. This type of model could be relevant for the recently detected PeVatrons in our Galaxy.

Current observations at radio wavelengths and X-rays provide us with information to constrain the parameters of our physical model, both for the core of the system and for the thermal and nonthermal lobes. Although we cannot in principle guarantee a direct relation between the observations and the acceleration of CRs to PeV energies, the results we have obtained for the emission produced by the system (which fits the observational data well), correspond to a set of parameters (see Table 2) that lead to the acceleration of protons to ultra-high energies at the reverse shock in the lobes, where the size of the accelerator is one of the relevant quantities involved in the process.

The LHAASO catalog lists 43 galactic ultra-high energy (UHE, E > 0.1 PeV) gamma-ray sources (Cao et al. 2024). Observations of UHE gamma rays provide a key clue to the acceleration of PeV protons and thus to the origin of galactic CRs. The detection of gamma rays above 10–100 TeV is crucial for tracing the accelerators (Vernetto & L.H.A.A.S.O. Collaboration 2016; Cao & Lhaaso Collaboration 2022) and determining the highest energies of CRs in the Galaxy.

We propose that yet-undetected super-Eddington MQs could be some of these sources (whose exact nature remains unknown). In the case of S26, the power of the system may allow protons to reach energies greater than 1 PeV, as we have shown. A key feature of such powerful jets relies on the final size of the accelerator in the reverse shock at the lobes. In our model, the maximum energy of the relativistic protons from the terminal region is limited by the diffusive escape (Fig. 4, bottom right), which depends on the square of the size of the acceleration region. A powerful jet could allow the acceleration region to be larger, imposing a cutoff to the protons at quite high energies. However, this scenario is conditioned by the physics in the system, which is quantified by several not well-known parameters, such as the efficiency of the acceleration processes in the different environments and the strength of the magnetic field in the acceleration regions. Nevertheless, our results suggest that super-Eddington MQs could potentially serve as PeVatron sources as long as they maintain their super-Eddington regime (Bykov et al. 2024).

Furthermore, if a Galactic super-Eddington MQ is located in a star-forming region, it could very easily go undetected if its jets are misaligned with the line of sight. The central X-ray emission would be obscured by the presence of the opaque wind, and proton-dominated jets tend to be dark most of the way (see e.g., Reynoso et al. 2008; Romero & Vila 2008). Only the impact of the jets on the environment would produce thermal X-rays and extended radio emission, but in large stellar associations there are many extended sources of uncertain origin. The protons would diffuse until they reach a suitable target. Then a gamma-ray source would appear with a significantly higher emissivity than the average gamma-ray production rate in the Galactic disk (Aharonian & Atoyan 1996). However, the origin of these PeV protons could be significantly displaced from the gamma-ray source itself, making identification difficult (see Aharonian & Atoyan 1996, for a discussion of such a scenario). A similar scenario was explored by Bosch-Ramon et al. (2005a), where the authors investigated hadronic MQs located in dense regions of the ISM.

The ultra-relativistic protons (E_p > 1 PeV) accelerated by the reverse shock in the lobes of S26 are injected into the ISM with a total power of ∼10³⁶ erg s⁻¹. As we just mentioned, they then diffuse and can reach regions of higher density, such as nearby molecular clouds. The energy distribution of protons at a distance D and time t can be determined by solving the diffusion equation (Ginzburg & Syrovatskii 1964). In the spherically symmetric case, it reads

$$\begin{array}{c}\hfill \frac{\partial f}{\partial t}=\frac{D}{R^2}\frac{\partial }{\partial R}{R}^2\frac{\partial f}{\partial R}+\frac{\partial }{\partial E}(Pf)+Q,\end{array}$$(25)

where f ≡ f(E, R, t) is the distribution function of particles, P = −(dE/dt) the continuous energy loss rate, Q ≡ Q(E, R, t) the source function, and D ≡ D(E) = D₀E^δ is the diffusion coefficient. Atoyan et al. (1995) found a solution to this equation, as a function of the diffusion radius, that is, the radius of the sphere up to which the particles of energy E_p effectively propagate during the time t after they left the acceleration site (Bosch-Ramon et al. 2005a).

The fraction of CRs that penetrate a cloud depends on a number of factors, including the magnetic field, the diffusion coefficient, advection, and energy losses due to ionization of neutral hydrogen atoms (Cesarsky & Volk 1978; Morlino & Gabici 2015; Dogiel et al. 2018; Ivlev et al. 2018). Approximately ∼0.1 of the proton population can produce gamma-ray emission through neutral pion decay if they interact with a nearby target. This, in turn, results in the conversion of up to 10% of the injected power into gamma rays if the cloud is dense enough. This means an emission of ∼10³⁴ erg s⁻¹ of very energetic gamma rays. If one of these sources is located in our Galaxy, approximately 10 kpc away, the flux on Earth would be about ∼10⁻¹² erg cm⁻² s⁻¹, which is easily detectable by LHAASO. Thus, it is possible that undetected super-Eddington MQs could be PeVatrons behind some of the unidentified UHE sources. The problem of gamma rays from the interaction of CRs with molecular clouds was addressed by Aharonian (1991), Combi & Romero (1995), Aharonian (2001), Gabici et al. (2007), among other authors. We plan to investigate this issue in detail elsewhere, following the approach of these papers, adapted to the case of super-Eddington MQs.

Other tracers that could help characterize the source could be neutrinos. The dominant scenario in the production of these particles will be pp interactions. Although the total power in neutrinos will be greater than in the case of SS433 (see e.g., Reynoso et al. 2008), S26 is thousands of times farther away, so the number of neutrinos reaching Earth will be lower. Furthermore, the only neutrino observatory currently operating at these energies is IceCube, which is located at the South Pole but detects neutrinos coming from the North, and S26 is only observable from the Southern Hemisphere.

As for the total power in CRs injected by a MQ similar to S26 in the host galaxy, although it is large, of the order of 5 × 10³⁸ erg s⁻¹ for E_p > 1 GeV (i.e, ∼5% of the CR power in our Galaxy), the lifetime of the super-Eddington phase (∼10⁵ yr) is not long enough to make a strong contribution to the average CR energy density. The contribution of a single source would be diluted by diffusion and cooling in the CR sea of the galaxy. We plan to continue the detailed study of CR production in super-Eddington MQs and their propagation and interaction in the future.

7. Conclusions

We developed a model to explain the puzzling features of the extragalactic MQ S26, which has the most powerful jets yet observed in accreting binaries. We propose that the four-order-of-magnitude difference between the jet and apparent disk luminosities in S26 is caused by the complete absorption of the disk radiation by the wind ejected from the super-Eddington disk. We suggest that some of the energy driving the jet must be extracted from the ergosphere of a high-spin BH.

We also reduced and analyzed data from multi-epoch X-ray observations of the X-ray XMM-Newton space telescope that contain S26 in the field of view. The information obtained complements that already available from Chandra observations at higher resolution. We find that the XMM-Newton can be well fitted by a thermal blackbody spectrum with a peak temperature of ∼0.1 keV and a hard power-law tail of index Γ ∼ 1.5, with large uncertainties. We cannot rule out variability in this component. In our model, we propose that the power-law component is the result of particle acceleration by internal shocks formed near the base of the jet, just above the wind photosphere, in the region where the outflow begins to be matter-dominated. The X-rays are produced by the electron synchrotron mechanism. The blackbody component, on the other hand, is produced in the terminal region of the jets, as a thermal contribution from the forward shock formed at the hotspots. The radio emission observed with ATCA can be explained as electron-synchrotron radiation from electrons accelerated by the reverse adiabatic shock in the lobes. All this is in agreement with the available X-ray Chandra images and the ATCA radio maps obtained by Pakull et al. (2010) and Soria et al. (2010).

We calculate that the protons accelerated by the shocks at both acceleration sites reach energies greater than ∼1 PeV. Therefore, we propose that S26 can be considered a PeVatron MQ with a super-Eddington accreting BH. A similar but weaker MQ in our Galaxy could also produce very high energy particles. If it is hidden within a stellar association, it could inject significant power into CRs that could “illuminate” molecular clouds and other matter overdensities, giving rise to UHE gamma-ray sources such as those observed by LHAASO. The scenario is similar to that proposed long ago by Aharonian & Atoyan (1996).

Acknowledgments

The authors thank the anonymous referee for a careful and constructive review, and for his/her comments, which have improved this paper. This work was supported by grant PIP 0554 (CONICET). LA acknowledges the Universidad Nacional de La Plata for the education received, Pablo Sotomayor for useful comments, and the 8th HEPRO organizers. GER was funded by PID2022-136828NB-C41/AEI/10.13039/501100011033/ and through the “Unit of Excellence María de Maeztu 2020-2023” award to the Institute of Cosmos Sciences (CEX2019-000918-M).

References

All Tables

All Figures

	Fig. 1. X-ray multi-epoch observations. XSPEC CONSTANT*TBABS*(BBODY+POWERLAW) fitting model of the spectrum of S26 using observations 0748390901 (top), 0804670301 (middle) and 0853981001 (bottom) taken by XMM-Newton. The black, red, and green colors represent the PN, MOS1, and MOS2 cameras, respectively.
In the text

Fig. 2. Conceptual scheme of the MQ S26, not to scale. The stellar-mass BH and the Wolf-Rayet star are in the core. The extended lobes are at the ends of the two opposite jets, where shocks form at the hotspots. There, the backward shock accelerates particles to relativistic energies, producing synchrotron radio emission, and the forward radiative shock produces optically thin thermal X-ray emission. Nonthermal X-ray emission is produced by the electron-synchrotron mechanism at the base of the jet in the core, close to the BH but above the photosphere of the disk-driven wind.

In the text

	Fig. 3. Photosphere of the disk-driven wind for an accretion rate of 10 ṀEdd. The axes are in Cartesian coordinates and units of the gravitational radius. The height zphoto is 1777 rg ( ≈ 2.6 × 109 cm) above the BH and its top temperature is ∼3.9 × 105 K.
In the text

Fig. 4. Cooling timescales of particles accelerated at the base of the jet (left) and at the terminal region (right). We show the energy gains and losses for electrons in the top panel and the same for protons in the bottom panel. In the internal shock at the base of the jet, the particles cool mainly by the synchrotron mechanism due to the strong magnetic field, reaching energies of 10 GeV (electrons) and 50 PeV (protons). At the reverse shock in the terminal region of the jet, electrons also cool mainly by synchrotron radiation, reaching energies of 5 PeV, while protons reach energies of ∼10 PeV. The efficiency of the acceleration in all cases is ηacc ∼ 0.1.

In the text

Fig. 5. Nonthermal emission from the jet. Left: Nonthermal SEDs of radiation produced at the internal shock in the jet, close to the BH but above the disk-driven wind photosphere. The dominant mechanism is electron synchrotron. For energies > 1 GeV the emission is completely attenuated due to internal absorption by pair creation. Dashed lines represent the unabsorbed emission. Right: Nonthermal SEDs of the radiation produced at the reverse shock in the terminal region of the jet. Electron-synchrotron radiation is the dominant process at low energies while proton-proton dominates at high and very high energies. No attenuation affects the radiation produced at the lobes.

In the text

Fig. 6. Luminosity of S26 as a function of the energy for our physical model with the parameters listed in Table 3. We plot the thermal and nonthermal SEDs of the system. The thermal contributions are the photosphere of the disk-driven wind and the optically thin emission in the lobes. Nonthermal emission is produced in the internal shocks at the base of the jet (orange solid line, jet + counterjet) and in the reverse shock of the lobes (purple solid line, both lobes). Individual nonthermal processes are shown with dashed and dotted-dashed lines (for details see Fig. 5). We plot radio data from ATCA (red dots, Soria et al. 2010) corresponding to the emission of the two lobes combined, and X-ray data from XMM-Newton (blue dots, this work) corresponding to the core and lobes, for obsID 0748390901. Sensitivities of the Fermi, CTA, SWGO and LHAASO gamma-ray detectors are also shown for two different distances: that from S26 (3.9 Mpc, solid colored lines), and a galactic distance of ∼10 kpc (black dashed lines).

In the text

Fig. 7. Exploration of the space parameter of the model. We plot SEDs for sets of input parameters different from those of Table 3 and Fig. 6. The dotted lines correspond to the SEDs of the different specific scenarios for the recollimation shock at the base (b) of the jet, (Sb1 − Sb5). The dashed lines correspond to the SEDs of the scenarios of the reverse shock at the lobes (l) of the jet, (Sl1 − Sl5). The corresponding colors for both shocks are cyan, blue, pink, red, and orange for models 1–5, respectively. The different model parameterizations are described in Sect. 5.4. Solid black lines represent our initial scenario. The observational data are plotted in the same way as in Fig. 6.

In the text