# Acceleration and transport of relativistic electrons in the jets of the microquasar SS 433

F Aharonian, F Ait Benkhali, J Aschersleben, H Ashkar, M Backes, V.

Barbosa Martins, R Batzofin, Y Becherini, D Berge, K Bernlohr, et al.

► To cite this version:

F Aharonian, F Ait Benkhali, J Aschersleben, H Ashkar, M Backes, et al. Acceleration and transport of relativistic electrons in the jets of the microquasar SS 433. Science, 2024, 383 (6681), pp.402-406. 10.1126/science.adi2048. hal-04425512

HAL Id: hal-04425512

https://hal.science/hal-04425512v1

Submitted on 17 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Acceleration and transport of relativistic electrons in
the parsec-scale jets of the microquasar SS 433

H.E.S.S. Collaboration*†

SS 433 is a microquasar, a stellar binary system with collimated relativistic jets. We observed SS 433 in gamma rays using the High Energy Stereoscopic System (H.E.S.S.), finding an energy-dependent shift in the apparent position of the gamma-ray emission of the parsec-scale jets. These observations trace the energetic electron population and indicate the gamma rays are produced by inverse-Compton scattering. Modelling of the energy-dependent gamma-ray morphology constrains the location of particle acceleration and requires an abrupt deceleration of the jet flow. We infer the presence of shocks on either side of the binary system at distances of 25 to 30 parsecs and conclude that self-collimation of the precessing jets forms the shocks, which then efficiently accelerate electrons.

<footnote>*laura.olivera-nieto@mpi-hd.mpg.de, brian.reville@mpi-hd.mpg.de, michelle.tsirou@desy.de, jim.hinton@mpi-hd.mpg.de, contact.hess@hess-experiment.eu</footnote>

<footnote>†H.E.S.S. Collaboration authors and affiliations are listed in the supplementary materials</footnote>
SS 433 (V1343 Aql) is a binary system comprising a compact object, likely a black hole (1–3), and a type A supergiant star (4). Accretion onto the black hole causes it to launch a pair of jets moving in opposite directions at approximately one quarter of the speed of light *c* (5–7), almost perpendicular to our line of sight (8). The jets precess with a half-opening angle of 20° and a period of 162 days (9–12). Adopting the distance measurement of 5.5 kiloparsecs (kpc) (7), optical and radio observations have shown that the precessing jets extend to distances of ~10⁻³ pc (13) and ~0.1 pc (7) from the black hole, respectively. X-ray emission reappears 25 pc from the binary (Figure 1), indicating collimated flows (the outer jets) on larger scales, which emit X-ray photons via non-thermal processes (14–17).

The outer jets terminate ~100 pc from the black hole (14), where they deform the surrounding radio nebula (known as W 50 or SNR G039.7-02.0) which is thought to be the supernova remnant associated with the formation of the compact object in SS 433 (18). The morphology of W 50 indicates that the opening angle of the outer jets is considerably smaller than the 20° precession angle of the inner jets (19); the origin of this discrepancy is unknown (20). The lack of apparent change in the measured positions of radio filaments in the jet termination regions over a 33 year period provides an upper limit on their velocity of < 0.023c (21, 22), though it is unclear whether the radio filaments trace the jets’ flow. Bright X-ray synchrotron knots have been observed in the outer jets but the temporal baseline and angular resolution was insufficient to determine their velocity. The dynamics of the outer jets and their termination process are poorly understood.

Since the initial X-ray detection of the outer jets (23), several attempts have been made to probe the non-thermal aspects and internal dynamics of the eastern (16, 24) and western (17) outer jets. However, observations of the X-ray synchrotron emission alone cannot resolve variations in the distribution of accelerated particles. The intensity of synchrotron emission is approximately proportional to the number density of accelerated electrons and the energy density of the magnetic field; the latter is poorly constrained. X-ray emitting electrons can also up-scatter low energy photons to the gamma-ray regime via the Inverse Compton scattering process. This process directly traces the population of high-energy electrons, since the diffuse low-energy photon distribution in the Galaxy is expected to be smooth on the spatial scale of the outer jets (25, 26). Previous observations of TeV gamma-rays emitted by the outer jets of SS 433 (27) indicate that the same energetic electrons responsible for the X-ray emission also produce gamma rays via inverse Compton scattering (28, 29). However, the angular resolution was insufficient to determine the emission regions and therefore the source of the energetic particles.
<figure><img src="image_2.png" /><figcaption>Figure 1: **Gamma-ray observations of SS 433.** A: Significance map of the H.E.S.S. observations at energies > 0.8 TeV (color bar). Cyan contours show the X-ray emission (14, 15). White crosses indicate locations of X-ray regions discussed in the text, w1, w2, e1, e2 and e3. Significance is for the H.E.S.S. excess counts above the background before accounting for statistical trials and after subtraction of the extended source HESS J1908+063 (subtraction shown in Figure S2). The map has been smoothed with a top-hat function of radius 0.1°. The white circle indicates the 68% containment region of the H.E.S.S. point-spread function (PSF). The green cross indicates the position of Fermi J1913+0515 and the green circle is its uncertainty. B: Orange circular points show our observed spectral energy distribution of the gamma-ray emission from the western jet. The brown square point is from previous observations (27). Error bars indicate the combined statistical (1σ) and systematic uncertainties; downward arrows indicate upper limits at 95% confidence. The solid line is the best-fitting power-law function, with dark and light shaded regions indicating the statistical and systematic uncertainties, respectively. C: Same as panel B but for the eastern jet. The regions from which the spectra shown in panels B and C were extracted are shown in Figure S2B.</figcaption></figure>
## H.E.S.S. observations of SS 433

We imaged the outer jets of SS 433 at TeV energies using the H.E.S.S. array of imaging atmospheric Cherenkov telescopes. The observations totalled over 200 hours of exposure time. They were analysed using previously described methods that were optimised for faint sources and improved performance at the highest energies (30). The extended source HESS J1908+063 (MGRO J1908+06) contaminates part of the SS 433 jet, so was modelled then subtracted from the data (Figures S2 and S3). The resulting gamma-ray image (Figure 1A) shows two regions of gamma-ray emission at the known positions of the eastern and western jets, with peak statistical significances of $7.8\sigma$ and $6.8\sigma$, respectively. No significant ($>5\sigma$) emission is detected from the central binary or the eastern termination region (Figure 1A), as we expect since the X-ray emission from those regions is predominantly thermal (14, 31). Fermi J1913+0515 (right ascension = $288.28 \pm 0.04^\circ$, declination = $5.27 \pm 0.04^\circ$) is a GeV gamma-ray source found to pulsate with a period consistent with the jet precession (32), suggesting a connection with the SS 433 system. No significant TeV emission is detected from this source (26). Figure 1B-C shows the measured spectral energy distributions of each of the jets.

To investigate the energy-dependence of the gamma-ray emission, we split the full energy range into three bands (0.8 to 2.5, 2.5 to 10 and $>10$ TeV), which were selected to have approximately the same gamma-ray excess counts over the background in each band. Figure 2 shows the significance maps for each band. We detect significant ($>5\sigma$) gamma-ray emission along both jets for the two highest energy bands. In the lowest energy band we find lower-significance evidence of emission at $4.4 \pm 0.7\sigma$ for the eastern and western jets, respectively. Gamma-ray emission $>10$ TeV appears only at the base of the outer jets visible in X-ray for both the eastern and western jets. In contrast, lower energy gamma rays have their peak surface brightnesses at locations further along each jet, except for the lowest-energy band on the eastern side. In the latter case, no significant emission is detected inside the X-ray jet region and evidence for emission appears close to the outer jet base (Figure 2A). In the western jet, the best-fit positions of the gamma-ray emission in each energy band have distances from the central binary (Table S4) that differ from each other by $0.97\sigma$ and $2.6\sigma$ when comparing adjacent energy bands, and by $5.3\sigma$ when comparing the lowest and highest energy band. The equivalent values for the eastern jet are $2.6\sigma$, $3.3\sigma$ and $0.1\sigma$. Our significance calculations include both systematic and statistical sources of uncertainty and a trials factor correction (26).

## Location of the particle acceleration

We interpret the offsets between the emission in different energy bands as indicating that transport of particles in the outer jets is dominated by the bulk jet flow (advection) and not the random scattering of the particles on magnetic field fluctuations (diffusion). The energy-dependent morphology then reflects an energy-dependent particle energy loss timescale. We infer that
<figure><img src="image_3.png" /><figcaption>Figure 2: Gamma-ray observations in different energy bands. Same as Figure 1A, but split into gamma-ray energy bands of (A) 0.8 to 2.5 TeV, (B) 2.5 to 10 TeV and (C) >10 TeV.</figcaption></figure>

the emission arises from relativistic electrons, and not hadrons, because the loss timescale for hadronic processes depends only very weakly on particle energy (33). The dominant energy loss mechanisms for high-energy electrons is likely to be synchrotron cooling. We conclude that the observed gamma-ray emission is the result of inverse Compton scattering (33, 34) of photons by high-energy electrons. Iron and other heavy nuclei are known to be present in the jet (35), so they might also be accelerated in the same region, but our observations cannot be used to constrain their presence (see Supplementary Text).

The shorter cooling time of higher-energy electrons limits the distance from the acceleration site within which they can radiate, because they are transported away by either diffusion or advection. The absence of emission above 10 TeV at the location of the X-ray knots (e2 and w2, Figure 1A) indicates that they cannot be sites of particle acceleration to TeV energies, contradicting previous interpretations (27, 36). Instead, the concentration of emission above 10 TeV at the base of the X-ray emission from the outer jets indicates this region is the site of particle acceleration to very high energies. We interpret the energy-dependent position of the gamma-ray emission in the jets of SS 433 as a consequence of the cooling and transport of particles that are accelerated at the base of the outer jets. Figure 3 shows a schematic diagram of our proposed interpretation.

## Modelling the outer jet dynamics

Previous studies have shown that the jets are launched from the black hole with initial velocities of $v_1 \approx 0.26c$ (5–7). We combined the distances between the gamma-ray excess regions in different energy bands with the electron cooling timescales (26) to determine the velocity $v_0$ of the outer jets at their base, approximately 25 pc away from the central binary. This calculation
<figure><img src="image_4.png" /><figcaption>**Figure 3: Schematic diagram of our model.** Thick black lines roughly outline the X-ray emission (grey shading) from the central region and the outer jets on the plane of the sky (rotated from the orientation in Figure 1). The precessing jet is launched with velocity $u_1 \approx 0.26c$ (purple spirals) and travels until it encounters a shock discontinuity (cyan bars), at the base of the outer jets. Our 1D model injects electrons continuously at the outer jet base, with an energy spectrum derived from a fit to multi-wavelength observations from the outer jets (Table S5). The injected electrons lose energy due to radiative losses, which affects their spectrum (indicated by the insets above the diagram). Particles are transported along the jets by the combination of diffusion (not depicted) and advection along the jet axis coordinate $z$ with the jet flow at velocity $v_{jet}(z)$ (purple arrows). The velocity at the base of the jets behind the shock $v_0$ can be determined by fitting the model prediction to the H.E.S.S. data. We assume the jet flow decelerates as it moves away the jet base, indicated by the purple curves below the diagram. We also considered the alternative case of a constant velocity jet (see Supplementary Text).</figcaption></figure>
requires us to assume a spatial dependence of the deceleration of the jets as a function of the dis-
tance from the central binary. We used the observed opening angle of the jets in X-ray images to
determine the deceleration profile by assuming the jet flow is incompressible (Figure S14). We
also considered a jet propagating with constant velocity, under different energy loss assump-
tions, which leads to consistent values of $v_0$ (see Supplementary Text). Our observations cannot
distinguish between the different jet propagation scenarios considered.

We model the energy-dependent morphology of the gamma-ray emission using a one-
dimensional Monte Carlo simulation which includes radiation and cooling of particles as they
are transported along the jet (26). The model injects electrons at the base of the outer jet with an
energy spectrum assumed to be of the form $dN/dE \propto E^{-\Gamma_e} \exp(-\frac{E}{E_{cut}})$ (Figure 3), where N
is the number of electrons, E their energy and $\Gamma_e$ and $E_{cut}$ the spectral index and cutoff energy,
respectively. We determine the best-fitting parameters of the injected electron spectrum and the
average local magnetic field strength from the multi-wavelength spectral energy distribution
of each outer jet separately. The value of $E_{cut}$ is not constrained by the data; we find only a
lower limit of > 200 TeV at 68% confidence level (C.L.). The model assumes the injection is
continuous for 10 000 years, this timescale being constrained by the combination of existing
GeV gamma-ray flux upper limits (37) and the measured TeV gamma-ray flux (Figure S11).
This electron injection timescale is consistent with previous dynamical estimates for the age of
the W50/SS433 complex, which range between 10 000 and 100 000 yr (19, 20). The simulation
evolves the electron population numerically in discrete time steps. In each step, electrons
are advected with the local jet velocity, then diffuse along the jet axis (neglecting transverse
diffusion) and cool radiatively. This leads to an energy- and spatially-dependent electron
distribution, from which we calculate one-dimensional profiles of the resulting non-thermal
emission. We find that the resulting spatial distribution in the gamma-ray range only weakly
depends on the parameters of the injected particle distribution (26).

Using the H.E.S.S. data, we derive spatial profiles of the gamma-ray flux along the axis
joining both outer jets through the central binary in the same three energy bands used in
Figure 2. We fitted the resulting model emission profile to the data with $v_0$ and the diffusion
coefficient, the latter assumed to be spatially uniform, as free parameters. The injected electron
spectrum parameters are fixed to the values obtained from the fit to the multi-wavelength
spectral energy distributions described above. We assume the same starting velocity for both
the eastern and western jet. The best-fitting value is $v_0 = (0.083 \pm 0.026_{stat} \pm 0.010_{syst})c$.
The systematic uncertainty is derived from the choice of parameters for the injected electron
spectrum (26). Figure 4 shows the gamma-ray spatial profiles and the best-fitting model.
<figure><img src="image_5.png" /></figure>

Figure 4: **Gamma-ray flux profiles along the jets compared with the model prediction.** Data points indicate the measured gamma-ray flux in spatial bins of 0.14° along the axis joining both jets through the central binary in the same three energy bands (panels A, B and C) as Figure 2. Error bars indicate the combined statistical (1σ) and systematic uncertainties. Solid lines show the prediction of our best-fitting 1D model. The shaded areas represent the combined statistical uncertainty of the best-fitting parameters. Dashed grey vertical lines represent the positions of the X-ray regions e1, e2, w1 and w2 (see Figure 1), which are labelled in panel C. The top axis assumes a distance to the system of 5.5 kpc (7).
Interpretation as a standing shock

Our modelling shows the data are consistent with the presence of a particle accelerator, likely a shock, at the base of the SS 433 outer jets which is capable of accelerating particles to very high energies. Our lower limit on $E_{cut}$ indicates the acceleration of electrons to energies >200 TeV (68% C.L.). At the inferred magnetic field strength of approximately 20 µG (Table S5), to keep up with cooling the acceleration rate must be close to the theoretical maximum, assuming diffusive shock acceleration (26, 38). Therefore the jet flow cannot have decelerated much from its inferred launch velocity of $0.26c$ prior to reaching the shock, because if it had, the particle acceleration could not compete with radiative losses at electron energies above several hundred TeV (Figure S13). The velocity we infer at the base of the outer jets $v_0$ is a fraction $\chi = 0.319 \pm 0.10_{stat} \pm 0.039_{syst}$ of the jet launch velocity. This is compatible with the velocity ratio expected for a sub-relativistic shock, which is $\chi \approx 0.25$ (39, 40). A shock at this location is consistent with the spatial coincidence between the position of the highest-energy gamma-ray emission and the location of the recently reported sharp X-ray reappearance of the X-ray emission (16, 17). This region has previously been interpreted as the acceleration site, but without involving shocks (16). Here we have shown that if the advection in the jet flow is taken into account, the observations are consistent with the shock acceleration scenario. Our observations also constrain the velocity of the shock, which would have needed to advance a small distance ($\ll 10$ pc) in the lifetime of the TeV gamma-ray emitting electrons (Figure S13).

There is no single model that has yet reproduced all the observational features of SS 433 (19, 41). While simulations can account for the observed difference in opening angle between the inner and outer jets due to the action of the ambient medium (20, 42, 43), this process would take place near the binary and does not result in the observed sharp transitions or shocks at 25 to 30 pc. The mirrored reappearance of the jets at this distance implies a physical significance to this radius, though there is no further observational evidence to indicate that this location is special. Radio observations of the jet launch region reveal the presence of an equatorial outflow perpendicular to the axis of the jets (44). Such an outflow has previously been argued to result in the formation of a quasi-spherical shock at distances of tens of pc from the binary (45). However, the X-ray shell that would be produced by such a shock has not yet been detected.

The proximity of SS 433 to Earth allows us to investigate shock physics and associated non-thermal processes in mildly relativistic jets. These insights can be applied both to other microquasars (46) as well as to the larger and more distant jets launched from the centres of other galaxies, in which jet sub-structure cannot be resolved at high energies (47). Our results imply that shocks forming within jets accelerate particles at close to the theoretical maximum energy (33, 48). Thus microquasars could be major contributors to the measured Galactic cosmic-ray flux at PeV energies, while extra-galactic jets could reach the EeV regime of ultra-high-energy cosmic rays (see Supplementary Text).
References

1. E. Seifina, L. Titarchuk, Astrophys. J. 722, 586 (2010).

2. A. M. Cherepashchuk, K. A. Postnov, A. A. Belinski, Mon. Not. R. Astron. Soc. 485, 2638 (2019).

3. A. M. Cherepashchuk, A. A. Belinski, A. V. Dodin, K. A. Postnov, Mon. Not. R. Astron. Soc. 507, L19 (2021).

4. T. C. Hillwig, D. R. Gies, Astrophys. J. 676, L37 (2008).

5. B. Margon, H. C. Ford, S. A. Grandi, R. P. S. Stone, Astrophys. J. 233, L63 (1979).

6. S. S. Eikenberry, et al., Astrophys. J. 561, 1027 (2001).

7. K. M. Blundell, M. G. Bowler, Astrophys. J. 616, L159 (2004).

8. D. H. Roberts, et al., Astrophys. J. 719, 1918 (2010).

9. A. C. Fabian, M. J. Rees, Mon. Not. R. Astron. Soc. 187, 13P (1979).

10. M. Milgrom, Astron. & Astrophys. 76, L3 (1979).

11. B. Margon, et al., Astrophys. J. 230, L41 (1979).

12. S. Fabrika, Astrophys. Space Phys. Res. 12, 1 (2004).

13. K. Davidson, R. McCray, Astrophys. J. 241, 1082 (1980).

14. W. Brinkmann, B. Aschenbach, N. Kawai, Astron. & Astrophys. 312, 306 (1996).

15. S. Safi-Harb, H. Ögelman, Astrophys. J. 483, 868 (1997).

16. S. Safi-Harb, et al., Astrophys. J. 935, 163 (2022).

17. K. Kayama, et al., Publications of the Astronomical Society of Japan 74, 1143 (2022).

18. R. Elston, S. Baum, Astron. J. 94, 1633 (1987).

19. P. T. Goodall, F. Alouani-Bibi, K. M. Blundell, Mon. Not. R. Astron. Soc. 414, 2838 (2011).

20. M. G. Bowler, R. Keppens, Astron. & Astrophys. 617, A29 (2018).

21. P. T. Goodall, K. M. Blundell, S. J. Bell Burnell, Mon. Not. R. Astron. Soc. 414, 2828 (2011).

22. H. Sakemi, R. Omae, T. Ohmura, M. Machida, Publ. Astron. Soc. Jpn 73, 530 (2021).
23. M. G. Watson, R. Willingale, J. E. Grindlay, F. D. Seward, Astrophys. J. 273, 688 (1983).

24. W. Brinkmann, G. W. Pratt, S. Rohr, N. Kawai, V. Burwitz, Astron. & Astrophys. 463, 611 (2007).

25. C. C. Popescu, et al., Mon. Not. R. Astron. Soc. 470, 2539 (2017).

26. Materials and methods are available as supplementary materials.

27. A. U. Abeysekara, et al., Nature 562, 82 (2018).

28. P. Bordas, V. Bosch-Ramon, J. M. Paredes, International Journal of Modern Physics D 19, 749 (2010).

29. T. Sudoh, Y. Inoue, D. Khangulyan, Astrophys. J. 889, 146 (2020).

30. L. Olivera-Nieto, H. X. Ren, A. M. W. Mitchell, V. Marandon, J. A. Hinton, The European Physical Journal C 82 (2022).

31. H. L. Marshall, C. R. Canizares, N. S. Schulz, Astrophys. J. 564, 941 (2002).

32. J. Li, et al., Nature Astronomy 4, 1177 (2020).

33. F. A. Aharonian, Very high energy cosmic gamma radiation : a crucial window on the extreme Universe (World Scientific Book, 2004).

34. G. R. Blumenthal, R. J. Gould, Reviews of Modern Physics 42, 237 (1970).

35. S. Migliari, R. Fender, M. Méndez, Science 297, 1673 (2002).

36. S. S. Kimura, K. Murase, P. Mészáros, Astrophys. J. 904, 188 (2020).

37. K. Fang, E. Charles, R. D. Blandford, Astrophys. J. 889, L5 (2020).

38. L. O. Drury, Reports on Progress in Physics 46, 973 (1983).

39. L. D. Landau, E. M. Lifshitz, Fluid Mechanics (Pergamon Press, 1987).

40. J. G. Kirk, P. Duffy, Journal of Physics G Nuclear Physics 25, R163 (1999).

41. T. Ohmura, et al., Astrophys. J. 910, 149 (2021).

42. R. Monceau-Baroux, O. Porth, Z. Meliani, R. Keppens, Astron. & Astrophys. 561, A30 (2014).

43. R. Monceau-Baroux, O. Porth, Z. Meliani, R. Keppens, Astron. & Astrophys. 574, A143 (2015).
44. K. M. Blundell, A. J. Mioduszewski, T. W. B. Muxlow, P. Podsiadlowski, M. P. Rupen, *Astrophys. J.* **562**, L79 (2001).

45. A. Konigl, *Mon. Not. R. Astron. Soc.* **205**, 471 (1983).

46. D. Kantzas, et al., *Mon. Not. R. Astron. Soc.* **500**, 2112 (2021).

47. H. E. S. S. Collaboration, et al., *Nature* **582**, 356 (2020).

48. A. M. Hillas, *Annu. Rev. Astron. Astrophys.* **22**, 425 (1984).

49. particle-transport-1d code (2023). doi:10.5281/ZENODO.10229720.

50. H. E. S. S. Collaboration, et al., *Astron. & Astrophys.* **457**, 899 (2006).

51. H. E. S. S. Collaboration, et al., *Astron. & Astrophys.* **612**, A1 (2018).

52. MAGIC Collaboration, et al., *Astron. & Astrophys.* **612**, A14 (2018).

53. T. Ashton, et al., *Astroparticle Physics* **118**, 102425 (2020).

54. B. Bi, et al., *37th International Cosmic Ray Conference. 12-23 July 2021*. Berlin (2022), p. 743.

55. S. Ohm, C. van Eldik, K. Egberts, *Astroparticle Physics* **31**, 383 (2009).

56. R. D. Parsons, J. A. Hinton, *Astroparticle Physics* **56**, 26 (2014).

57. L. Olivera-Nieto, A. M. W. Mitchell, K. Bernlöhr, J. A. Hinton, *European Physical Journal C* **81**, 1101 (2021).

58. A. Donath, et al., *Astron. & Astrophys.* **678**, A157 (2023).

59. A. Donath, et al., *Gammapy v0.19* (2021). doi:10.5281/ZENODO.5721467.

60. L. Mohrmann, et al., *Astron. & Astrophys.* **632**, A72 (2019).

61. T. P. Li, Y. Q. Ma, *Astrophys. J.* **272**, 317 (1983).

62. M. de Naurois, L. Rolland, *Astroparticle Physics* **32**, 231 (2009).

63. W. Reich, E. Fuerst, C. G. T. Haslam, P. Steffen, K. Reif, *Astron. & Astrophys. Suppl.* **58**, 197 (1984).

64. W. Reich, E. Fuerst, P. Reich, K. Reif, *Astron. & Astrophys. Suppl.* **85**, 633 (1990).

65. E. Furst, W. Reich, P. Reich, K. Reif, *Astron. & Astrophys. Suppl.* **85**, 691 (1990).
66. V. V. Davydov, V. F. Esipov, A. M. Cherepashchuk, *Astronomy Reports* **52**, 487 (2008).

67. H. E. S. S. Collaboration, *et al.*, *Astron. & Astrophys.* **672**, A103 (2023).

68. S. Gillessen, Sub-bogenminuten-genaue positionen von tev-quellen mit h.e.s.s., Ph.D. thesis, Ruprecht-Karls-Universität Heidelberg (2004).

69. P. Bordas, R. Yang, E. Kafexhiu, F. Aharonian, *Astrophys. J.* **807**, L8 (2015).

70. X.-N. Sun, R.-Z. Yang, B. Liu, S.-Q. Xi, X.-Y. Wang, *Astron. & Astrophys.* **626**, A113 (2019).

71. Y. Xing, Z. Wang, X. Zhang, Y. Chen, V. Jithesh, *Astrophys. J.* **872**, 25 (2019).

72. K. Rasul, P. M. Chadwick, J. A. Graham, A. M. Brown, *Mon. Not. R. Astron. Soc.* **485**, 2970 (2019).

73. J. Hahn, 34th International Cosmic Ray Conference (ICRC2015) (2015), vol. 34 of *International Cosmic Ray Conference*, p. 917.

74. J. Hahn, C. Romoli, M. Breuhaus, Gamera: Source modeling in gamma astronomy, Astro-physics Source Code Library, record ascl:2203.007 (2022).

75. D. L. Band, *Publ. Astron. Soc. Pac.* **99**, 1269 (1987).

76. A. M. Atoyan, F. A. Aharonian, *Mon. Not. R. Astron. Soc.* **278**, 525 (1996).

77. P. O. Lagage, C. J. Cesarsky, *Astron. & Astrophys.* **125**, 249 (1983).

78. A. W. Strong, I. V. Moskalenko, V. S. Ptuskin, *Annual Review of Nuclear and Particle Science* **57**, 285 (2007).

79. Y. Su, *et al.*, *Astrophys. J.* **863**, 103 (2018).

80. D. L. Band, M. A. Gordon, *Astrophys. J.* **338**, 945 (1989).

81. Q.-C. Liu, Y. Chen, P. Zhou, X. Zhang, B. Jiang, *Astrophys. J.* **892**, 143 (2020).

82. H. Yamamoto, *et al.*, *Publ. Astron. Soc. Jpn* **74**, 493 (2022).

83. P. Boumis, *et al.*, *Mon. Not. R. Astron. Soc.* **381**, 308 (2007).

84. S. Heinz, R. Sunyaev, *Astron. & Astrophys.* **390**, 751 (2002).

85. R. P. Fender, T. J. Maccarone, Z. van Kesteren, *Mon. Not. R. Astron. Soc.* **360**, 1085 (2005).

86. A. J. Cooper, D. Gaggero, S. Markoff, S. Zhang, *Mon. Not. R. Astron. Soc.* **493**, 3212 (2020).
87. G. J. Escobar, L. J. Pellizza, G. E. Romero, Astron. & Astrophys. **665**, A145 (2022).

88. R. M. Arnason, H. Papei, P. Barmby, A. Bahramian, M. D. Gorski, Mon. Not. R. Astron. Soc. **502**, 5455 (2021).

**Acknowledgements**

**Acknowledgements** We thank Jian Li and Ke Fang for providing the reduced Fermi-LAT dataset and GeV flux points, respectively. We appreciate the excellent work of the technical support staff in Berlin, Zeuthen, Heidelberg, Palaiseau, Paris, Saclay, Tübingen and in Namibia in the construction and operation of the equipment. This work benefited from services provided by the H.E.S.S. Virtual Organisation, supported by the national resource providers of the EGI Federation.

**Funding** The support of the Namibian authorities and of the University of Namibia in facilitating the construction and operation of H.E.S.S. is gratefully acknowledged, as is the support by the German Ministry for Education and Research (BMBF), the Max Planck Society, the German Research Foundation (DFG), the Helmholtz Association, the Alexander von Humboldt Foundation, the French Ministry of Higher Education, Research and Innovation, the Centre National de la Recherche Scientifique (CNRS/IN2P3 and CNRS/INSU), the Commissariat à l'énergie atomique et aux énergies alternatives (CEA), the U.K. Science and Technology Facilities Council (STFC), the Irish Research Council (IRC) and the Science Foundation Ireland (SFI), the Knut and Alice Wallenberg Foundation, the Polish Ministry of Education and Science, agreement no. 2021/WK/06, the South African Department of Science and Technology and National Research Foundation, the University of Namibia, the National Commission on Research, Science & Technology of Namibia (NCRST), the Austrian Federal Ministry of Education, Science and Research and the Austrian Science Fund (FWF), the Australian Research Council (ARC), the Japan Society for the Promotion of Science, the University of Amsterdam and the Science Committee of Armenia grant 21AG-1C085.

**Author Contributions:** L. Olivera-Nieto performed the H.E.S.S. analysis and M. Tsirou performed the cross-check analysis. N. Tsuji analysed the X-ray data. B. Reville, L. Olivera-Nieto and J. Hinton performed the interpretation and modelling. The manuscript was prepared by L. Olivera-Nieto, B. Reville and M. Tsirou. S. Wagner is the collaboration spokesperson. All other H.E.S.S. collaboration authors contributed to the design, construction and operation of H.E.S.S., the development and maintenance of data handling, data reduction or data analysis software. All authors meet the journal's authorship criteria and have reviewed, discussed, and commented on the results and the manuscript.

**Competing Interests:** The authors declare that they have no competing interests.

**Data and Materials Availability:** The H.E.S.S. data are available at: https://www.mpi-hd.mpg.de/HESS/pages/publications/auxiliary/2023NSS433/index.html. This includes the counts and background maps, the subsequently derived significance and flux maps (Figures 1A and 2), the flux pro-
files (Figure 3) and the raw data and derived flux points from the spectral mea-
surements of the jets (Figure 1B-C). Our modelling code is available at https:
//github.com/LauraOlivera/particle-transport-1D (49).

**Supplementary Material:** Authors and affiliations, Materials and Methods, Supplementary Text, Figure S1-S16, Tables S1-S5, References (50-88)
Supplementary Material for

Acceleration and transport of relativistic electrons in the
parsec-scale jets of the microquasar SS 433

H.E.S.S. Collaboration*

*Correspondence to: contact.hess@hess-experiment.eu; Laura Olivera-Nieto (laura.olivera-nieto@mpi-hd.mpg.de), Brian Reville (brian.reville@mpi-hd.mpg.de), Jim Hinton (jim.hinton@mpi-hd.mpg.de), Michelle Tsirou (michelle.tsirou@desy.de)

This PDF file includes:

- Authors and affiliations
- Materials and Methods
- Supplementary Text
- Figs. S1 to S16
- Tables S1 to S5
H.E.S.S. Collaboration authors and affiliations

F. Aharonian<sup>1,2</sup>, F. Ait Benkhali<sup>3</sup>, J. Aschersleben<sup>4</sup>, H. Ashkar<sup>5</sup>, M. Backes<sup>6,7</sup>, V. Bar-bosa Martins<sup>8</sup>, R. Batzofin<sup>9</sup>, Y. Becherini<sup>10,11</sup>, D. Berge<sup>8,12</sup>, K. Bernlöhr<sup>2</sup>, B. Bi<sup>13</sup>, M. Böttcher<sup>7</sup>, C. Boisson<sup>14</sup>, J. Bolmont<sup>15</sup>, M. de Bony de Lavergne<sup>16</sup>, J. Borowska<sup>12</sup>, M. Bouyahiaoui<sup>2</sup>, M. Breuhaus<sup>2</sup>, R. Brose<sup>1</sup>, A.M. Brown<sup>17</sup>, F. Brun<sup>18</sup>, B. Bruno<sup>19</sup>, T. Bulik<sup>20</sup>, C. Burger-Scheidlin<sup>1</sup>, S. Caroff<sup>16</sup>, S. Casanova<sup>21</sup>, R. Cecil<sup>22</sup>, J. Celic<sup>19</sup>, M. Cerruti<sup>10</sup>, T. Chand<sup>7</sup>, S. Chandra<sup>7</sup>, A. Chen<sup>23</sup>, J. Chibueze<sup>7</sup>, O. Chibueze<sup>7</sup>, G. Cotter<sup>17</sup>, S. Dai<sup>24</sup>, J. Dam-ascene Mbarubucyeye<sup>8</sup>, A. Djannati-Ataï<sup>10</sup>, A. Dmytriiev<sup>7</sup>, V. Doroshenko<sup>13</sup>, K. Egberts<sup>9</sup>, S. Einecke<sup>25</sup>, J.-P. Ernenwein<sup>26</sup>, M. Filipovic<sup>24</sup>, G. Fontaine<sup>5</sup>, S. Funk<sup>19</sup>, S. Gabici<sup>10</sup>, S. Ghafourizadeh<sup>3</sup>, G. Giavitto<sup>8</sup>, D. Glawion<sup>19</sup>, J.F. Glicenstein<sup>18</sup>, G. Grolleron<sup>15</sup>, L. Haerer<sup>2</sup>, J.A. Hinton<sup>2*</sup>, W. Hofmann<sup>2</sup>, T. L. Holch<sup>8</sup>, M. Holler<sup>27</sup>, D. Horns<sup>22</sup>, M. Jamrozy<sup>28</sup>, F. Jankowsky<sup>3</sup>, A. Jardin-Blicq<sup>29</sup>, V. Joshi<sup>19</sup>, I. Jung-Richardt<sup>19</sup>, E. Kasai<sup>6</sup>, K. Katarzyński<sup>30</sup>, R. Khatoon<sup>7</sup>, B. Khélifi<sup>10</sup>, S. Klepser<sup>8</sup>, W. Kluźniak<sup>31</sup>, Nu. Komin<sup>23</sup>, K. Kosack<sup>18</sup>, D. Kostunin<sup>8</sup>, A. Kundu<sup>7</sup>, R.G. Lang<sup>19</sup>, S. Le Stum<sup>26</sup>, F. Leitl<sup>19</sup>, A. Lemière<sup>10</sup>, J.-P. Lenain<sup>15</sup>, F. Leuschner<sup>13</sup>, T. Lohse<sup>12</sup>, A. Luashvili<sup>14</sup>, J. Mackey<sup>1</sup>, D. Malyshev<sup>13</sup>, D. Malyshev<sup>19</sup>, V. Marandon<sup>18</sup>, P. Marchegiani<sup>23</sup>, A. Marcowith<sup>32</sup>, G. Martí-Devesa<sup>27</sup>, R. Marx<sup>3</sup>, A. Mehta<sup>8</sup>, A. Mitchell<sup>19</sup>, R. Moderski<sup>31</sup>, L. Mohrmann<sup>2</sup>, A. Montanari<sup>3</sup>, E. Moulin<sup>18</sup>, T. Murach<sup>8</sup>, K. Nakashima<sup>19</sup>, M. de Naurois<sup>5</sup>, J. Niemiec<sup>21</sup>, A. Priyana Noel<sup>28</sup>, S. Ohm<sup>8</sup>, L. Olivera-Nieto<sup>2*</sup>, E. de Ona Wilhelmi<sup>8</sup>, M. Ostrowski<sup>28</sup>, S. Panny<sup>27</sup>, M. Panter<sup>2</sup>, R.D. Parsons<sup>12</sup>, G. Peron<sup>10</sup>, D.A. Prokhorov<sup>33</sup>, G. Pühlhofer<sup>13</sup>, M. Punch<sup>10</sup>, A. Quirrenbach<sup>3</sup>, P. Reichherzer<sup>18</sup>, A. Reimer<sup>27</sup>, O. Reimer<sup>27</sup>, H. Ren<sup>2</sup>, M. Renaud<sup>32</sup>, B. Reville<sup>2*</sup>, F. Rieger<sup>2</sup>, G. Rowell<sup>25</sup>, B. Rudak<sup>31</sup>, H. Rueda Ricarte<sup>18</sup>, E. Ruiz-Velasco<sup>2</sup>, V. Sahakian<sup>34</sup>, H. Salzmann<sup>13</sup>, A. Santangelo<sup>13</sup>, M. Sasaki<sup>19</sup>, J. Schäfer<sup>19</sup>, F. Schüssler<sup>18</sup>, U. Schwanke<sup>12</sup>, J.N.S. Shapopi<sup>6</sup>, H. Sol<sup>14</sup>, A. Specovius<sup>19</sup>, S. Spencer<sup>19</sup>, Ł. Stawarz<sup>28</sup>, R. Steenkamp<sup>6</sup>, S. Steinmassl<sup>2</sup>, C. Steppa<sup>9</sup>, K. Streil<sup>19</sup>, I. Sushch<sup>7</sup>, H. Suzuki<sup>35</sup>, T. Takahashi<sup>36</sup>, T. Tanaka<sup>35</sup>, A.M. Taylor<sup>8</sup>,
R. Terrier<sup>10</sup>, M. Tsrhou<sup>8*</sup>, N. Tsuji<sup>37</sup>, T. Unbehaun<sup>19</sup>, C. van Eldik<sup>19</sup>, M. Vecchi<sup>4</sup>, J. Veh<sup>19</sup>, C. Venter<sup>7</sup>, J. Vink<sup>33</sup>, T. Wach<sup>19</sup>, S.J. Wagner<sup>3</sup>, F. Werner<sup>2</sup>, R. White<sup>2</sup>, A. Wierzcholska<sup>21</sup>, Yu Wun Wong<sup>19</sup>, M. Zacharias<sup>3,7</sup>, D. Zargaryan<sup>1</sup>, A.A. Zdziarski<sup>31</sup>, A. Zech<sup>14</sup>, S. Zouari<sup>10</sup>, N. Żywucka<sup>7</sup>

<sup>1</sup>Dublin Institute for Advanced Studies, Dublin D02 XF86, Ireland

<sup>2</sup>Max-Planck-Institut für Kernphysik, Heidelberg D-69117, Germany

<sup>3</sup>Landessternwarte, Universität Heidelberg, Heidelberg D-69117, Germany

<sup>4</sup>Kapteyn Astronomical Institute, University of Groningen, Groningen 9747 AD, The Netherlands

<sup>5</sup>Laboratoire Leprince-Ringuet, École Polytechnique, Centre national de la recherche scientifique (CNRS), Institut Polytechnique de Paris, Palaiseau F-91128, France

<sup>6</sup>Department of Physics, University of Namibia, Windhoek 10005, Namibia

<sup>7</sup>Centre for Space Research, North-West University, Potchefstroom 2520, South Africa

<sup>8</sup>Deutsches Elektronen-Synchrotron (DESY), Zeuthen D-15738, Germany

<sup>9</sup>Institut für Physik und Astronomie, Universität Potsdam, Potsdam 14476, Germany

<sup>10</sup>Laboratoire Astroparticule et Cosmologie, Université de Paris, CNRS, Paris F-75013, France

<sup>11</sup>Department of Physics and Electrical Engineering, Linnaeus University, Växjö 351 95, Sweden

<sup>12</sup>Institut für Physik, Humboldt-Universität zu Berlin, Berlin D-12489, Germany

<sup>13</sup>Institut für Astronomie und Astrophysik, Universität Tübingen, Tübingen D-72076, Germany

<sup>14</sup>Laboratoire Univers et Théories (LUT), Observatoire de Paris, Université Paris Sciences et Lettres, CNRS, Université de Paris, Meudon 92190, France

<sup>15</sup>Laboratoire de Physique Nucléaire et de Hautes Energies (LPNHE), Sorbonne Université, Université Paris Diderot, Université Paris Cité, Institut national de physique nucléaire et de physique des particules (IN2P3), CNRS, Paris F-75252, France
<sup>16</sup>Laboratoire d'Annecy de Physique des Particules, CNRS, IN2P3, Université Savoie Mont Blanc, Annecy 74000, France

<sup>17</sup>Department of Physics, University of Oxford, Oxford OX1 3RH, UK

<sup>18</sup>Institute for Research on the Fundamental Laws of the Universe (IRFU), Commissariat à l'énergie atomique et aux énergies alternatives (CEA), Université Paris-Saclay, Gif-sur-Yvette F-91191, France

<sup>19</sup>Erlangen Centre for Astroparticle Physics, Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen D-91058, Germany

<sup>20</sup>Astronomical Observatory, The University of Warsaw, Warsaw 00-478, Poland

<sup>21</sup>Instytut Fizyki Jądrowej, Polska Akademia Nauk (PAN), Kraków 31-342, Poland

<sup>22</sup>Institut für Experimentalphysik, Universität Hamburg, Hamburg D-22761, Germany

<sup>23</sup>School of Physics, University of the Witwatersrand, Johannesburg 2050, South Africa

<sup>24</sup>School of Science, Western Sydney University, Penrith NSW 2751, Australia

<sup>25</sup>School of Physical Sciences, University of Adelaide, Adelaide 5005, Australia

<sup>26</sup> Centre de Physique des Particules de Marseille (CPPM), Aix Marseille Université, CNRS/IN2P3, Marseille 13288, France

<sup>27</sup>Institut für Astro- und Teilchenphysik, Leopold-Franzens-Universität Innsbruck, Innsbruck A-6020, Austria

<sup>28</sup>Obserwatorium Astronomiczne, Uniwersytet Jagielloński, Kraków 30-244, Poland

<sup>29</sup>Laboratoire de de Physique des deux Infinis (LP2I), Université Bordeaux, CNRS, Gradignan F-33170, France

<sup>30</sup>Institute of Astronomy, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University, Torun 87-100, Poland

<sup>31</sup>Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, Warsaw 00-716, Poland
<sup>32</sup>Laboratoire Univers et Particules de Montpellier, Université Montpellier, CNRS/IN2P3,
Montpellier F-34095, France

<sup>33</sup>Gravitation and Astroparticle Physics Amsterdam (GRAPPA), Anton Pannekoek Institute for Astronomy, University of Amsterdam, Amsterdam 1098 XH, The Netherlands

<sup>34</sup>Yerevan Physics Institute, Yerevan 375036, Armenia

<sup>35</sup>Department of Physics, Konan University, Higashinada-ku Kobe 658-8501, Japan, Japan

<sup>36</sup>Kavli Institute for the Physics and Mathematics of the Universe, The University of Tokyo
Institutes for Advanced Study, The University of Tokyo, Kashiwa Chiba 277-8583, Japan

<sup>37</sup>The Institute of Physical and Chemical Research (RIKEN), Wako Saitama 351-0198, Japan
*Corresponding author. Emails: Laura Olivera-Nieto (laura.olivera-nieto@mpi-hd.mpg.de),
Brian Reville (brian.reville@mpi-hd.mpg), Jim Hinton (jim.hinton@mpi-hd.mpg.de),
Michelle Tsirou (michelle.tsirou@desy.de)
# Materials and Methods

## The H.E.S.S. array

The H.E.S.S. array of imaging atmospheric Cherenkov telescopes is located in the Khomas Highland of Namibia, at an altitude of 1,835 m. H.E.S.S. is sensitive to gamma rays ranging from tens of GeV to tens of TeV. The array consists of five Cherenkov telescopes: four with mirror diameters of 12 m placed in a square configuration (designated CT1 to CT4) and a single telescope at the centre (CT5) with a mirror diameter of 28 m. The 12-m telescope array is sensitive to gamma rays of energies above several hundreds of GeV, while the central, large telescope is able to detect fainter Cherenkov emission, which in turn translates to lower gamma-ray energies (50).

## Observations of SS 433

The SS 433 region was observed by the H.E.S.S. array of telescopes as part of a Galactic plane survey (51), and again during two dedicated campaigns. The first dedicated campaign, between 2009 and 2011 did not detect SS 433 and has been described previously (52). A second campaign, between 2018 and 2021, collected around 150 h of data. We combined all three datasets, resulting in around 200 h of data. Table S1 lists the properties of each dataset we used. Some of this data was taken before an upgrade of the CT1-4 cameras (53) which we refer to as the HESS-I dataset, and the data taken after the upgrade is labelled HESS-IU. For HESS-IU, we make a further distinction between the observations taken with (HESS-IU-CT5) and without (HESS-IU) CT5. All the observations involving CT5 were performed after an upgrade of the CT5 camera in 2019 (54). The HESS-IU and HESS-IU-CT5 datasets overlap in time. Data-taking was split into observation runs usually spanning 28 minutes. The pointing positions of the observations comprising the HESS-I dataset are concentrated around HESS J1908+063 and the location of the western jet of SS 433. The pointing positions of the later observations were
chosen to achieve a uniform exposure over the entire field.

<table><caption>Table S1: **H.E.S.S. datasets.** Listed are the exposure time, mean zenith angle, start and end dates and array configuration for each of the used datasets. The exposure quoted corresponds to the maximum value of the livetime when the exposure from the different pointing positions are combined in the field of view.</caption><thead><tr><th>dataset</th><th>exposure (h)</th><th>mean zenith (deg)</th><th>start date</th><th>end date</th><th>CT5 included</th></tr></thead><tbody><tr><td>HESS-I</td><td>71.1</td><td>38.8</td><td>2005 Jun 3</td><td>2013 Aug 20</td><td>no</td></tr><tr><td>HESS-IU</td><td>33.5</td><td>44.2</td><td>2018 Aug 29</td><td>2021 Nov 3</td><td>no</td></tr><tr><td>HESS-IU-CT5</td><td>111.1</td><td>48.9</td><td>2020 Jun 19</td><td>2021 Nov 5</td><td>yes</td></tr></tbody></table>

## Event reconstruction and background rejection

H.E.S.S. records stereoscopic images of atmospheric showers produced by high-energy particles as they travel in the Earth's atmosphere. Gamma-ray-like events are selected by a Boosted Decision Tree classifier (55). The selected events are reconstructed using the Image Pixel-wise fit for Atmospheric Cherenkov Telescopes (ImPACT) algorithm (56), which uses a maximum-likelihood framework to fit a library of simulated templates to the data images. This process results in an estimate for the gamma-ray energy and direction, among other parameters. We only use information from CT1-4 for these steps. An extra step of background rejection is applied to the observation runs that include the CT5 telescope, which is the majority of runs taken in the 2019 to 2021 campaign. This method exploits the reduced threshold of the central telescope, which is more efficient at detecting faint emission from other particles, such as muons. Muons are produced in large numbers in hadronic showers, the main source of background for Cherenkov gamma-ray detectors (57). In this step, the CT5 image of events selected as gamma-ray-like based on their CT1-4 image is compared to the associated ImPACT prediction image of that event. Events for which the CT5 image differs according to a set of criteria from the expected template are rejected, which leads to a factor 3 to 4 improvement in background rejection. The method criteria, performance and implementation have been described
elsewhere (30).

## Data reduction

Subsequent steps in the data reduction and analysis were carried out using the gammapy software (58), version 0.19 (59). The selected gamma-ray-like events were binned into a three-dimensional data cube, consisting of sky-map of 6° width and 0.01° bin size, centred at the position of SS 433 and an energy axis with 22 bins equally spaced in logarithmic energy between 0.63 and 100 TeV. For each observation, a safe energy range was derived, defined as the range in which the energy bias was less than 10% (50). This criterion rejected the energy bin between 0.63 and 0.8 TeV in all our observations, so 0.8 TeV is the lowest event energy in our H.E.S.S. data.

Each observation run has an associated model of expected background counts which depends on radial offset from the pointing position and on reconstructed energy (60). These background models were derived from observation runs with mostly extra-galactic pointing positions in which expected gamma-ray sources were masked. The background model varies between runs due to different pointing altitude position and hardware changes. For each individual run in the dataset, the counts predicted by the model outside of an exclusion mask are fitted to those measured in the same region using two free parameters that determine the overall background normalisation and spectral shape. This procedure corrects the background model for possible variations due to atmospheric conditions and instrumental degradation (60). The exclusion mask was defined to cover known and expected gamma-ray sources in the SS 433 field of view. It is composed of a band of 2° height around the Galactic plane ($|b| \le 1°$ where $b$ is the Galactic latitude), a circle of 1.33° radius for the nearby extended source HESS J1908+063 (Figure S2) and two circles of 0.57° radius for the jets of SS 433.

The instrument response functions (IRFs) that describe the precision of the energy and di-
<figure><img src="image_6.png" /><figcaption>Figure S1: **Spectral energy distribution of HESS J1908+063.** Purple circular points show our observed spectral energy distribution of the gamma-ray emission from HESS J1908+063. Error bars indicate statistical (1σ) uncertainties. The solid purple line is the best-fitting log-parabola function. The dashed black line is the best-fitting power-law function.</figcaption></figure>

rection reconstruction, as well as the effective area of the detector were projected into multi-
dimensional sky-maps with the corresponding energy axes. The set of maps corresponding to
each observation were stacked by adding the counts and background maps and combining the
IRFs weighted by the exposure of each run, considering only the energy range determined as
safe for each observation. Significance maps were computed using a maximum-likelihood ratio
test based on the comparison of the number of measured counts to the expected background
counts (61). Both the background and counts map were smoothed with a top-hat kernel of ra-
dius 0.1°, chosen to be roughly the size of the H.E.S.S. point-spread-function (PSF). All data
reduction was confirmed by an independent reduction, acting as a cross check. The cross-check
analysis employed independent calibration, reconstruction and background suppression (62).

**Removal of contaminating source HESS J1908+063**

The bright extended source HESS J1908+063 is located less than 2° away from the posi-
tion of SS 433 and detected with significance of more than 10σ in the combined dataset.
<figure><img src="image_7.png" /><figcaption>Figure S2: **Subtraction of HESS J1908+063.** Same as Figure 1A, but (A) before and (B) after subtracting the emission from the nearby extended source HESS J1908+063. In panel B, the white circle indicates the 68% containment region of the model fitted to HESS J1908+063, and the white cross is its best-fitting position. In both panels, the solid white contours show radio emission from the W50 nebula (63–65). In panel B, the blue ellipses show the regions from where the spectral measurement of the jets is extracted (Figure 1B-C). The dashed line shows the axis across the jets used to derive the gamma-ray spatial profiles shown in Figures 4 and S3.</figcaption></figure>

To assess the degree of contamination into the SS 433 region, we fitted the emission from HESS J1908+063 with a combined spatial and spectral model using a maximum-likelihood framework in gammapy. The model of HESS J1908+063 has a single Gaussian component with curved spectrum described by a log-parabola function (Figure S1). The curvature is preferred to a simpler power-law at a significance of 6.4σ. Figure S2 shows the significance map of the full field of view, including the HESS J1908+063 region, before and after the best-fitting model is subtracted.

Figure S3 shows the gamma-ray flux profiles measured along the jets of SS 433 for different energy ranges, before and after the subtraction of HESS J1908+063. The contamination of HESS J1908+063 to the jets region is strongest at the tip of the western jet, located at around 0.5° from the binary. At that location, HESS J1908+063 contributes about ≈ 40% of the measured flux. This contribution quickly decreases towards the eastern jet, where it is negligible.
We include the subtraction of HESS J1908+063 in our estimation of systematic errors (see below).

## Background systematic effects

We estimate the systematic uncertainties in our background estimation by measuring the significance distribution outside of the exclusion mask. Figure S4 shows the distribution over the entire energy range, together with the best-fitting Gaussian function. The Gaussian has a width (1.0830 ± 0.0054)σ, higher than the expected value of exactly 1σ. Its mean is -(0.0421 ± 0.0054)σ, which is lower than the expected value of 0σ. After correcting for this effect, emission is detected at significances of 7.2σ and 6.3σ for the eastern and western outer jets, respectively.

For the significance maps in different energy (Figure 2), we derived equivalent distributions for each of the energy bands separately. The resulting histograms are shown in Figure S5, with widths of (1.0379 ± 0.0068)σ, (1.0718 ± 0.0062)σ and (1.018 ± 0.015)σ from low to high energy. The corrected maximum significance values when accounting for this effect are 4.2σ, 7σ and 5.8σ for the eastern jet from low to high energy and 4.5σ, 5.2σ and 6.5σ for the western jet in the same order.

## Spectra of the jets

Due to the energy-dependent morphology, assuming a single spatial model for the jets across all energies would introduce inaccuracies in the associated spectral model when fitting both components together. Therefore, we measured the spectra of the jets (after subtracting HESS J1908+063) without assuming a spatial model by extracting the spectral information inside two elliptical regions (Figure S2) large enough to completely contain the gamma-ray excess in each jet. Spectral models were fitted to the measured excesses in these regions using
<figure><img src="image_8.png" /></figure>

Figure S3: Gamma-ray flux profiles along the jets showing the contamination of HESS J1908+063. The data points represent the measured flux in spatial bins of 0.14° along the axis joining both jets through the central binary (Figure S2) for energies (A) above 0.8 TeV, (B) 0.8 to 2.5 TeV, (C) 2.5 to 10 TeV and (D) above 10 TeV. Squares and circles indicate the flux before and after subtracting HESS J1908+063. Error bars indicate the combined statistical (1σ) and systematic uncertainties. Circles in panels B-D are the same data as shown in Figure 4. The dashed lines show the flux of the HESS J1908+063 model at each location. The top axis assumes a distance to the system of 5.5 kpc (7).
<figure><img src="image_9.png" /><figcaption>Figure S4: **Background significance distribution.** The distribution of significance values (red histogram) is plotted for all locations outside of an exclusion mask covering all expected sources. The black dashed line is a Gaussian function fitted to the distribution (black line), Its mean and width are labelled. The deviation of this model from the expected mean of 0 and width of 1 is used to assess the systematic uncertainties in the significance values due to background subtraction.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure S5: **Background significance distribution in different energy bands.** Same as Figure S4 but split into the energy bands (A) 0.8 to 2.5 TeV, (B) 2.5 to 10 TeV and (C) above 10 TeV.</figcaption></figure>
<table><caption>Table S2: **Parameters of the power-law model fitted to the gamma-ray spectra of the jets.** Listed are the best-fitting and fixed parameters in the model: $E_0$ is the reference energy, $\phi_0$ is the amplitude at the reference energy and $\Gamma$ is the photon spectral index.</caption><thead><tr><th></th><th>φ<sub>0</sub><br>(10<sup>-13</sup> TeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup>)</th><th>E<sub>0</sub><br>(TeV)</th><th>Γ</th></tr></thead><tbody><tr><td>east</td><td>2.30 ± 0.58<sub>stat.</sub> ± 0.32<sub>syst.</sub></td><td>1</td><td>2.19 ± 0.12<sub>stat.</sub> ± 0.12<sub>syst.</sub></td></tr><tr><td>west</td><td>2.83 ± 0.70<sub>stat.</sub> ± 0.39<sub>syst.</sub></td><td>1</td><td>2.40 ± 0.15<sub>stat.</sub> ± 0.13<sub>syst.</sub></td></tr></tbody></table>

a maximum-likelihood framework. Models were compared via their test-statistic (TS) value, with the significance of the improvement in the description due to one additional parameter computed as $\sigma = \sqrt{\Delta TS}$. For both jets the differential photon spectrum $\frac{dN}{dE}$ is described by a power-law $\phi_0 (\frac{E}{E_0})^{-\Gamma}$ where $\Gamma$ the photon spectral index, $E_0$ the fixed reference energy and $\phi_0$ is the amplitude at the reference energy. We find no significant (> 3σ) evidence for curvature (1.5σ and 0.7σ for east and west) or an exponential cutoff (1.3σ and 0.15σ) in the spectrum. The best-fitting parameters are listed in Table S2; both jets have consistent parameters within the uncertainties. Flux points were derived by fitting a normalisation parameter in each energy bin, assuming the best-fitting spectral shape derived from the wider energy range. The resulting flux and best-fitting model spectra are shown in Figure 1.

## Spatial properties of the energy-integrated gamma-ray emission

The energy-integrated spatial properties of the gamma-ray emission were modelled using two elliptical Gaussian components. The elliptical description is preferred to a symmetrical Gaussian description by 5.8σ and 3.5σ for east and west, which in turn is preferred to a point-like description by 7.8σ and 4.7σ for east and west. The angle of the asymmetrical Gaussian model $\theta$ (degrees east from north) is fixed to that of the SS 433 X-ray jets ($\theta = -19^\circ$). Allowing the angle to vary during fitting of the model to the data results in $\theta = -16.2 \pm 3.5^\circ$ for the eastern excess, preferred by 0.74σ, and $\theta = -7.1 \pm 5.4^\circ$ for the western excess, preferred by 1.9σ. As neither is significant (> 3σ), we keep the angle fixed to the value from X-ray observations. The
<table><caption>Table S3: **Parameters of the elliptical model fitted to the spatial morphology.** Listed are the best-fitting and fixed parameters in the model: *l* and *b* are the Galactic longitude and latitude respectively of the central position, *r*<sub>maj</sub> and *r*<sub>min</sub> are the major and minor axis standard deviation sizes and *θ* is the angle of the jets, which was fixed in the fit. The physical sizes for *r*<sub>maj</sub> and *r*<sub>min</sub> are calculated for a distance of 5.5 kpc (7). Uncertainties are both systematic and statistical.</caption><thead><tr><th></th><th>unit</th><th><i>l</i></th><th><i>b</i></th><th><i>r</i><sub>maj</sub></th><th><i>r</i><sub>min</sub></th><th><i>θ</i></th></tr></thead><tbody><tr><td rowspan="2">east</td><td>deg</td><td>39.875 ± 0.018</td><td>-2.687 ± 0.027</td><td>0.205 ± 0.035</td><td>0.044 ± 0.014</td><td>-19</td></tr><tr><td>pc</td><td></td><td></td><td>19.68 ± 3.36</td><td>4.22 ± 1.35</td><td></td></tr><tr><td rowspan="2">west</td><td>deg</td><td>39.564 ± 0.013</td><td>-1.853 ± 0.027</td><td>0.134 ± 0.036</td><td>0.046 ± 0.015</td><td>-19</td></tr><tr><td>pc</td><td></td><td></td><td>12.86 ± 3.46</td><td>5.37 ± 1.44</td><td></td></tr></tbody></table>

parameters of the best-fitting elliptical Gaussian model for each jet are presented in Table S3. The spatial extension is significant (> 3σ) in both the major and minor axis directions of the ellipse.

## Morphology in different energy bins

Gamma-ray excesses are detected along both jets in all energy bins (Figure 2) with peak significance values (from low to high energy) of 4.4σ, 7.6σ, 5.9σ for the eastern jet, and 4.7σ, 5.6σ and 6.6σ for the western jet. The morphology of the measured gamma-ray excesses was fitted separately in each energy band, to compare the resulting best-fit positions. The gamma-ray excesses are found to be significantly (> 3σ) extended (not point sources) and are modelled with a symmetric Gaussian in most cases; however the western jet at >10 TeV is not significantly extended, so we used a point source model. The spectral parameters are fixed to those from the fitting to the full energy range (Table S2). We tested making these free parameters, which resulted in consistent spatial parameters in all cases [both with each other, and with the values from the full energy range (Table S2)]. The distance from the position of the SS 433 binary at Galactic coordinates *l*=39.694°, *b*=-2.245° for each energy range is shown in Figure S6. As a further check, the data were split into 3 bands with different energy boundaries and into different numbers (5, 10) of energy bands. In all cases, the resulting trend was consistent within the
<table><caption>Table S4: **Results of spatial model fitting different energy ranges.** Same as Table S3 but for three different energy ranges, each fitted separately. The distance from the model center to the nominal position of SS 433 $d_{SS\ 433}$ is listed both in degrees and in pc, the latter adopting a distance estimate of 5.5 kpc (7). Uncertainties are both systematic and statistical.</caption><thead><tr><th>side</th><th>energy<br>(TeV)</th><th><i>l</i><br>(deg)</th><th><i>b</i><br>(deg)</th><th><i>r</i><br>(deg)</th><th><i>d</i><sub>SS 433</sub><br>(deg)</th><th><i>d</i><sub>SS 433</sub><br>(pc)</th></tr></thead><tbody><tr><td rowspan="3">east</td><td>0.8 to 2.5</td><td>39.913 ± 0.044</td><td>-2.614 ± 0.047</td><td>0.125 ± 0.022</td><td>0.428 ± 0.046</td><td>41.148 ± 4.424</td></tr><tr><td>2.5 to 10</td><td>39.924 ± 0.018</td><td>-2.772 ± 0.021</td><td>0.085 ± 0.015</td><td>0.575 ± 0.021</td><td>55.212 ± 2.007</td></tr><tr><td>above 10</td><td>39.840 ± 0.031</td><td>-2.643 ± 0.038</td><td>0.013 ± 0.029</td><td>0.424 ± 0.037</td><td>40.693 ± 3.593</td></tr><tr><td rowspan="3">west</td><td>0.8 to 2.5</td><td>39.537 ± 0.024</td><td>-1.759 ± 0.033</td><td>0.080 ± 0.016</td><td>0.510 ± 0.032</td><td>48.946 ± 3.089</td></tr><tr><td>2.5 to 10</td><td>39.582 ± 0.024</td><td>-1.826 ± 0.037</td><td>0.095 ± 0.018</td><td>0.433 ± 0.037</td><td>41.590 ± 3.552</td></tr><tr><td>above 10</td><td>39.560 ± 0.010</td><td>-1.951 ± 0.011</td><td>-</td><td>0.323 ± 0.011</td><td>31.015 ± 1.038</td></tr></tbody></table>

statistical uncertainties.

## Gamma-ray spatial flux profiles

Spatial flux profiles along the jets in different energy bands were derived from the H.E.S.S. data. We used them to visualise the contamination from HESS J1908+063 (Figure S3) and for comparison with one-dimensional theoretical models (Figure 4). We constructed the profiles by defining an axis passing through both jets and the central binary as the line between the Galactic coordinates $(l, b)=(40.246^\circ, -3.695^\circ)$ and $(39.340^\circ, -1.295^\circ)$ (Figure S2). We defined a rectangular box of full width $0.7^\circ$ along this axis and sliced it into 19 perpendicular boxes of height $0.14^\circ$. The height of the boxes is chosen as a compromise between having enough signal in each slice and sampling the spatial shape of the gamma-ray emission, for comparison the 68% containment radius of H.E.S.S. PSF is ~$0.07^\circ$. The definition of these boxes does not consider the positions of, or distances between the emission centroids in the different energy bands (Table S4). In each of these boxes the excess counts were computed and a normalization was fitted assuming a power-law spectral shape with index 2.3, the mean of the measured indices of the eastern and western jets. Using instead the best-fitting value for either of the jets, of 2.2 and 2.4 for east and west, respectively (Table S2), gives consistent results. The resulting integrated flux profiles are shown in Figures 4 and S3. The outlier flux point seen at around $0.7^\circ$
<figure><img src="image_11.png" /><figcaption>Figure S6: **Distance of the peak gamma-ray emission from the central binary.** The distance from the best-fitting position of the gamma-ray emission centroid to the central binary is shown for both the eastern (green squares) and the western (orange circles) jets. In the energy range between 0.8 and 2.5 TeV, each jet is detected with a significance only above 4σ, which is why the measurements in this range are represented with a transparent point. Error bars indicate the combined statistical (1σ) and systematic uncertainties. The location of the base of the outer jets as inferred from X-ray maps (14) is marked with a dashed line for each side. The distances in parsec are calculated adopting a distance of 5.5 kpc to the system. (7)</figcaption></figure>
towards the eastern side above 10 TeV is the result of a background fluctuation outside the jet region (Figure 2C). Reducing the width of the box could prevent this outlier, but also excludes some of the gamma-ray flux, especially at the lower energies. Therefore, we choose to keep the box wide enough to contain the gamma-ray emitting regions entirely.

## Search for periodic variability

Fermi J1913+0515 is a GeV gamma-ray source reported to pulsate with a period consistent with the precession of the SS 433 jets (32). The position of Fermi J1913+0515 and its uncertainty are shown in Figure 1, with no significant emission at that location. However, due to the periodic nature of Fermi J1913+0515, it is possible that any TeV emission would only be detectable around certain phase ranges, like in the GeV range. We investigated the presence of TeV emission from Fermi J1913+0515 by phase-folding the H.E.S.S. observations. We adopted the jet precession period to be 162.250 days (66) and a starting time $T_0$ of JD 2443508.4098 (32) and derived a corresponding phase for each H.E.S.S. observation. We separated the observations into 8 phase bins and repeated the data reduction process for each group, resulting in a full map of the entire region in each case. No significant emission was found in the vicinity of Fermi J1913+0515 in any of the phase ranges considered, both for the full energy range or the three energy ranges used in Figure 2. No significant emission was found either when considering 4 and 2 bins of phase. A non-detection in the TeV range is consistent with the observation of energies only up to 10 GeV from Fermi J1913+0515 (32). No significant phase trend was found for any other part of the SS 433 system in any of the phase or energy ranges considered.

While a sub-threshold excess in the region between Fermi J1913+0515 and the X-ray contours can be seen in the phase-integrated significance map shown in Figure 1, no phase trend or variation is observed at this location. When split into different bins of phase, the sub-threshold excess at this location disappears, indicating no apparent relation with the jet precession phase.
Additionally, when the jet emission is modelled with an elliptical Gaussian, no significant (or near-threshold) emission is left at that location, meaning that our measured TeV gamma-ray counts from the eastern jet region are compatible with a single component aligned with the jet axis.

## Upper limits from the central source and e3 regions

No significant (> 5σ) excess of TeV gamma rays is detected at either the position of the binary system or at the X-ray jet termination region e3 (Figure 1A). Differential upper limits for both regions were computed assuming a spectral index of 2.7 following previous methods (52). The radius of the circular region where the points are extracted is 0.07° and 0.2° for the central and e3 regions, respectively, motivated by the size of the X-ray emission in each region. The resulting limits are shown in Figure S7, together with the previously most constraining upper limits (52) for the central source. The non-detection of TeV emission confirms that the X-ray emission from these two regions is predominantly of thermal origin (14, 31).

## Systematic uncertainties in model parameters

The systematic uncertainties in the model parameters are calculated with a Monte Carlo-based approach, in which the IRFs are randomly varied to generate random pseudo-datasets based on the best-fitting spatial and spectral models described above. These datasets are then re-fitted using the original, unmodified IRFs. The obtained spread in the fitted parameters then reflects their combined statistical and systematic uncertainty. This procedure is extensively described elsewhere (67). The resulting systematic uncertainties are consistent with previous estimates of the H.E.S.S. systematic errors (50). The contamination of HESS J1908+063 into the jet regions was assessed by not including that model component in half of the generated pseudo-datasets. The same procedure was used when calculating systematic errors for flux points. In this case,
<figure><img src="image_12.png" /><figcaption>Figure S7: **Upper limits for undetected regions.** 95% confidence level flux upper limits (downward arrows) derived at the position of the central binary (teal circles) and the e3 region (yellow squares). Gray circles are limits on the central source from previous observations (52).</figcaption></figure>

the resulting systematic uncertainty is of the same magnitude as the statistical one at energies
around 1 TeV, and quickly becomes negligible at higher energies. For certain fitted parameters,
such as the eccentricity and width of the elliptical Gaussian components, the distribution is not
broadened by the systematic effects considered in this estimate. This is also the case for the
best-fitting Galactic coordinates of the source components. However, these parameters may be
affected by other systematic effects neglected in the Monte Carlo approach. In particular, the
source positions are subject to a systematic uncertainty of the pointing position of the H.E.S.S.
telescopes, which is of the order of 10" - 20" (68). This value is around a factor 10 lower than
the statistical uncertainty on the positional measurements (see Table S3). We accounted for this
source of systematic uncertainty with a 20" systematic uncertainty included in the best-fitting
coordinates presented in Tables S3 and S4.
<figure><img src="image_13.png" /><figcaption>Figure S8: **X-ray flux points.** The X-ray flux of the eastern jet derived from XMM-Newton data is shown before (light blue upside down triangles) and after (dark blue triangles) the corrections described in the text. Same for the Chandra data of the western jet but with orange circles and yellow squares respectively. Error bar indicate statistical (1σ) uncertainties.</figcaption></figure>

## Multi-wavelength data

We used multi-wavelength (MWL) observations of the outer jets of SS 433 to determine the spectral energy distribution (SED) of both jets.

## Radio

The outer jets have not been detected in the radio band (hundreds of km to mm wavelengths), a range in which instead the shell of W 50 is a bright source. We used data from the 11 cm Effelsberg radio telescope survey (63–65) are used to determine the total flux inside the region where the gamma-ray spectra are extracted. This provides an upper limit of $7.084 \cdot 10^{-14}$ TeV s$^{-1}$ cm$^{-2}$ and $1.367 \cdot 10^{-13}$ TeV s$^{-1}$ cm$^{-2}$ to the radio flux coming from the eastern and western jets respectively, as they must at least be fainter than W 50.
## X-ray

X-ray flux was determined from data taken by the X-ray Multi-Mirrors-Newton (XMM-Newton) (24) and Chandra (17) space telescopes for the eastern and western jet, respectively. The flux points were extracted from the same spatial regions used for the TeV measurements, with the limitation of the different field of view from the different instruments. Only data above 2 keV is used to minimise the contribution of thermal X-ray emission and the effect of absorption. In the case of the Chandra data the effect of the small field-of-view is taken into account with a correction factor of order 10%, derived from the surface brightness of a soft X-ray image of the jets (24). The X-ray emission from the eastern jet includes that of the bright region labelled e2 (Figure 1A). The relative brightness of this region with respect to the rest of the X-ray emission from the eastern jet suggests that it might be the result of a local enhancement of the magnetic field. Because we are interested in the average magnetic field of the jets, we subtract the spectrum of this region (16); we discuss the effects of this correction below. The resulting X-ray flux points for both jets before and after the corrections are shown in Figure S8.

## GeV gamma rays

The presence of a bright pulsar in the field of view and uncertainties on the Galactic diffuse emission make the SS 433 region difficult to study at GeV energies, with different studies reaching conflicting conclusions (69–72). Evidence for (~4σ) an emission excess in the *Fermi* Large Area Telescope (*Fermi*-LAT) data has been reported (37) near the outer western jet, but not spatially coincident with it. The 95% containment upper limits for both the western and eastern jets (37) are reproduced in Figures S10 and S11.
Multi-wavelength spectral model

We modelled the MWL spectral energy distribution of each jet using a single electron popula-
tion injected continuously at the base of the outer jets. This model is fitted to the MWL data
described above and the H.E.S.S. data assuming two one-zone (eastern and western) scenarios.
Transport effects within the jets are not included, because we assume particles cool at the same
rate everywhere in the jet. We use the gamera package (73, 74) to model the temporal evolu-
tion of the electron distribution and the resulting radiation. The injected electron spectrum is
parameterised as a power-law with index Γe and an exponential cutoff at Ecut. Synchrotron
radiation produced by electrons with a spectral index Γe is expected to follow a photon spectral
index of (Γe+1)/2. The measurement of an X-ray photon index of roughly 1.5 at the base of
both the western (17) and eastern (16) outer jets consequently implies Γe = 2. Therefore, we
fix the value of the electron spectral index to Γe = 2 in our fit. The power injected in electrons
above 10 MeV is parameterised as a fraction α of the jet kinetic power, the latter taken to be
10³⁹erg s⁻¹. The injected electron spectrum is assumed to have the same parameters for both
jets.

The injected electrons lose energy via two processes: synchrotron emission in a uniform
magnetic field *B* and IC scattering of photons from a uniform ambient radiation field. The
photon emission from the binary peaks in the ultra-violet (UV) range (12); IC scattering on UV
photons is severely suppressed due to the Klein-Nishina effect (34), is included in the gamera
package implementation of IC scattering. For TeV gamma-ray emission, far-infrared (FIR)
photons provide instead the dominant scattering field. At distances of more than 25 pc from
the binary, the dominant FIR photon field is not that of the binary or W50 (75) but the diffuse
IR background. We use the combination of an axisymmetric Galactic model (25) derived from
observations at similar distances from the Galactic Centre as the outer jets and the cosmic
microwave background as target field. The energy densities of the assumed ambient fields are
<figure><img src="image_14.png" /><figcaption>Figure S9: **Model of the ambient radiation fields.** A: Model energy density of the background radiation fields (25) at the position of the central binary (teal solid line), the eastern (yellow dash-dotted line) and the western (red dashed line) jets. B: Energy densities shown in panel A normalised by the model at the central binary highlighting the differences between the two jets.</figcaption></figure>

shown in Figure S9. The radiation fields differ by less than 10% at the locations of each jet with respect to the one at the location of the central binary. The synchrotron photons produced by the electrons would in principle contribute to the target field (Synchrotron-self-Compton (76)), but this contribution is negligible as inferred from the predicted synchrotron flux (Figure S10).

The models were fitted to the TeV and X-ray data for each jet, using a maximum likelihood approach, and considering the GeV and radio upper limits. For simplicity, the magnetic field strength is the only parameter allowed to differ between the eastern and western sides. The best-fitting values for the parameters are shown in Table S5. We do not find significant (> 3σ) evidence for an exponential cutoff in the electron spectrum, meaning that if such a feature is present, it must be at higher electron energies than those probed by H.E.S.S. observations. Therefore we derive a 68% C.L. lower limit for the cutoff energy in the electron spectrum of $E_{\text{cut}} > 200$ TeV. The power required to reproduce the observed emission is around 0.13% of the jet kinetic power of $10^{39}$ erg s⁻¹. The best-fitting magnetic field strengths are roughly 19
and 21 µG for the eastern and western jet, respectively, in agreement with estimates derived from X-ray observations (16). For reference, the Poynting flux of the jet assuming a velocity ~0.1c and these values of the magnetic field is ~10³⁷ erg s⁻¹, around one order of magnitude more than the power in the relativistic electrons. The resulting MWL spectral model is shown in Figure S10 for both jets.

<table><caption>Table S5: **Best-fitting parameters from the multi-wavelength emission model.** Best-fit values of the injected electron spectral index ($\Gamma_e$), exponential cutoff ($E_{cut}$), power as a fraction of the jet kinetic power ($\alpha$) and the magnetic field ($B$). Uncertainties are statistical only. The third column indicates whether each parameters was shared for both jets, or allowed to vary between them.</caption><thead><tr><th></th><th>east</th><th>west</th><th>shared</th><th>fixed</th></tr></thead><tbody><tr><td>Γ<sub>e</sub></td><td>2</td><td>2</td><td>yes</td><td>yes</td></tr><tr><td>E<sub>cut</sub> (TeV)</td><td>&gt;200</td><td></td><td>yes</td><td>no</td></tr><tr><td>α</td><td>(1.287 ± 0.029)·10<sup>-3</sup></td><td></td><td>yes</td><td>no</td></tr><tr><td>B (μG)</td><td>19.5±2.7 21.1±1.8</td><td></td><td>no</td><td>no</td></tr></tbody></table>

If the X-ray flux from the e2 region in the eastern jet is not subtracted (Figure S8), the parameter $\alpha$ varies between the eastern and western jets, resulting in $\alpha_{east} = (1.597 \pm 0.027) \cdot 10^{-3}$, $B_{east} = 20.0 \pm 2.6$ µG, $\alpha_{west} = (1.184 \pm 0.073) \cdot 10^{-3}$ and $B_{west} = 20.8 \pm 1.8$ µG. While these values of the magnetic field are consistent (within the uncertainties) with those reported in Table S5, including the emission from the e2 region requires more power to be injected into the electrons on the eastern side. Providing more power to the eastern side would not affect our conclusions, so for the subsequent discussions we refer to the model derived from the corrected X-ray data.

The combination of the injection index required by the X-ray data, the GeV upper limits and the steepness of the H.E.S.S. spectra constrains for how long the electrons are injected (Figure S11), which in turn constrain the age of the jets (or its recent activity). Assuming continuous injection with the parameters listed in Table S5, injection is constrained to last more than 1000 yr, but less than 30000 yr, with an age of between 3000 and 10000 yr yielding
<figure><img src="image_15.png" /><figcaption>Figure S10: Multi-wavelength spectral energy distribution of the SS 433 outer jets and best-fitting model. Flux points and upper limits (downward arrows) derived from the gamma-ray emitting regions in the (A) eastern and (B) western outer jets in the radio (stars), X-ray (triangles), GeV gamma-rays (upside-down triangles) and TeV gamma-rays including both this work (circles) and previous observations by the High-Altitude Water Cherenkov Observatory (HAWC, squares). The best-fitting model for each of the jets is shown with a solid line. The residual is shown in panels C and D.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure S11: **Model predictions for different values of the injection time.** Coloured lines show emission as a function of the assumed injection time (colour bar) for the eastern (A) and western (B) jets, assuming the parameters listed in Table S5. TeV measurements and GeV upper limits are the same as in Figure S10.</figcaption></figure>

better agreement with the observations (Figure S11). This result is consistent with previous theoretical and numerical studies, which placed the age of the W50/SS 433 complex, ranging between 10 000-100 000 yr, e.g. (19, 20).

**Comparison with X-ray morphology**

To compare the observed gamma-ray morphology to the X-ray morphology, we focus on the eastern jet using the X-ray images covering the entire eastern jet in a wide range of photon energies from the Roentgensatellit (ROSAT) (14) and XMM-Newton (16, 24) space telescopes. The morphology of the X-ray emission from the eastern jet is also energy-dependent (16), with predominantly hard (2 to 12 keV) X-ray emission detected from the e1 region, a mixture of soft (0.5 to 2 keV) and hard emission from the e2 region and predominantly soft thermal X-ray emission from the e3 region. Comparing the gamma-ray morphology to the X-ray morphology therefore requires a choice of which X-ray energies to consider. We do so based on the model for the emission described above (Table S5). Figure S12 shows the model contribution of electrons of different energies to the total SED. The electrons responsible for the X-ray emission in the
XMM-Newton energy range (2 to 7 keV) mostly produce gamma rays of energies above 10 TeV. The lower energy (0.5 to 2 keV) X-ray emission detected by ROSAT is mostly produced by lower energy electrons, which are responsible for the remaining gamma-ray emission measured by H.E.S.S. below approximately 10 TeV. We therefore constructed two spatial templates using flux maps measured by ROSAT in the 0.5-2 keV band and XMM-newton in the 2-7 keV band.

The spatial templates were used as a model that is fitted to the morphology of the gamma-ray excess in the relevant energy range. This was done by multiplying the templates by a flux normalisation factor, which is a free parameter. Besides this template, an additional model component for HESS J1908+063 was used, also with the flux normalisation as a free parameter. The ROSAT template was fitted to the H.E.S.S. data below 10 TeV and the XMM-Newton template was fitted to the H.E.S.S. data above 10 TeV. In both cases, the residual map contained no excesses with significance greater than $5\sigma$. However, there is evidence for ($\sim 4\sigma$) an excess in the residual map in the lower energy model around the e1 position, likely due to the e2 region being several times brighter than the surrounding emission in the ROSAT image. Less significant emission can be seen in the residual of the higher energy model, where sub-threshold ($\sim 2.5\sigma$) positive and negative excesses are found around e1 and e2 respectively. If the ROSAT template is fitted to the entire H.E.S.S. energy range, a significant ($> 5\sigma$) excess remains at the base of the outer jet around e1, because the emission from that region arises predominantly from the highest energy electrons. We conclude that the observed energy-dependent morphology in the X-ray and gamma-ray bands are consistent with a shared origin due to the acceleration of particles at the base of the outer jet and their transport in the jet flow. The observed sub-threshold discrepancies may be explained with a modest local increase of the magnetic field in the e2 region.
<figure><img src="image_17.png" /><figcaption>Figure S12: **Contribution from electrons of different energies to the eastern SED model.** Coloured solid lines show emission from electrons with different energies (colour bar). The range of electron energies depicted was chosen to overlap with the H.E.S.S. observations. The model parameters used are listed in Table S5. The ROSAT energy range is indicated with a blue band. Flux points, upper limits and the total SED (dashed blue line) are the same as in Figure S10A.</figcaption></figure>
<figure><img src="image_18.png" /></figure>

Figure S13: **Calculated cooling and acceleration timescales.** The cooling timescale $t_{\text{cool}}$ is plotted as a function of electron energy (black dashed line) assuming the best-fitting values for the magnetic field (Table S5) and the same ambient radiation field as used in the MWL modelling. The contributions from IC ($t_{\text{IC}}$) and synchrotron ($t_{\text{synch}}$) losses are shown in orange and blue, respectively. Plotted are the values for the eastern jet; the western jet is almost identical. Dotted lines display the acceleration timescales ($t_{\text{acc}}$) for different values of the efficiency $\eta$, assuming a jet velocity prior to the shock of $v_1 = 0.26c$. The maximum electron energy range allowed by the lower limit on the cutoff energy of the electron spectrum is indicated with a grey band. The region defined by the combination of this energy range, the cooling time and acceleration time is indicated with vertical hatching. Adiabatic loss time scales ($t_{\text{ad}}$) assuming a jet with constant velocity, are depicted with green dot-dashed lines for different values of $v_0$, a proxy for the advection time. For a decelerating constant-density jet or a jet moving with both constant section and velocity, adiabatic losses are negligible.
## Shock acceleration timescale

We consider the timescale requirements for electrons to be accelerated at a shock located at the base of the outer jet. The acceleration timescale as a function of energy $t_{acc}$ is plotted in Figure S13, computed assuming test-particle diffusive shock acceleration theory (38):

$$t_{acc} = \frac{3}{u_1 - u_2} \left( \frac{D_1}{u_1} + \frac{D_2}{u_2} \right) \approx \frac{8}{\eta} \frac{D_{Bohm}}{u_1^2}, \qquad (S1)$$

where $u_1$, $u_2$, $D_1$ and $D_2$ are the velocities in the shock frame and diffusion coefficients in the upstream and downstream media, respectively. The parameter $\eta$ is the ratio of an electron's mean free path to its gyroradius ($r_g$), a measure of the scattering efficiency with respect to idealised magnetised transport where $D = \frac{1}{\eta}D_{Bohm} = \frac{1}{3\eta}r_g c$. Equation S1 assumes approximately equal upstream and downstream residence times. The case with $\eta = 1$ corresponds to the limiting minimum diffusion coefficient for magnetised transport (i.e. for which particles can still be considered to undergo helical trajectories) and therefore sets a lower limit on the acceleration time (77). Our emission model (Table S5) set a lower limit for the energy cutoff of the electron spectrum of >200 TeV. This implies values of $\eta (u_1/0.26c)^2 \gg 0.01$ for electrons to compete with radiative cooling losses. Upstream velocities $u_1$ much less than $0.1c$ would require $\eta > 1$, implying non-magnetised transport. In this regime, additional care is needed in applying Equation S1, but in general would make the inferred energies unreachable.

## The effect of particle transport

We interpret the observed energy-dependent position of the gamma-ray emission in the jets of SS 433 as a consequence of the combination of particle cooling timescales (Figure S13) and advection with the jet flow, and thus can be used to constrain the internal dynamics of the outer jets. We test this assumption using a one-dimensional Monte Carlo simulation which models the transport of particles via the combined action of advection and diffusion. We used gamera
to calculate the radiation and cooling of particles as they are transported. The simulation was
run from some initial time $t = 0$ to an assumed duration of the electron injection $t = t_{\text{inj}}$.

We define the coordinate $z$ as the position along the jet axis. At each time-step $t_{\text{step}}$, elec-
trons are injected at the base of the outer jets, assumed to be $z = 0$ with a spectral shape
described by a power-law with an exponential cutoff. The parameters of the injected electron
distribution and the value of the magnetic field strengths are those obtained from the multi-
wavelength model (Table S5 and Figure S10). The injected particles undergo radiative losses in
the magnetic and soft radiation fields, both assumed to be uniform and isotropic. This leads to a
change in the electron spectrum at each $t_{\text{step}}$. Particles move along the $z$ axis by the combination
of diffusion and advection. Diffusion parallel to the jet direction is included as a random Gaus-
sian smearing with scale $\sqrt{2Dt_{\text{step}}}$, where $D$ is the spatially homogeneous diffusion coefficient.
Diffusion is parameterised as $D = D_{100} (\frac{E}{100\text{TeV}})^{1/3}$, where $D_{100}$ is the value of the diffusion
coefficient at 100 TeV. We do not consider diffusion orthogonal to the jet axis. The assumed
energy dependence of the diffusion coefficient follows Kolmogorov scaling $D \propto E^{1/3}$, con-
sistent with inferences for Galactic cosmic-ray transport (78). The advection flow is described
by a spatially-dependent velocity $v_{\text{jet}}(z)$. Diffusion moves particles in both the increasing and
decreasing $z$ directions but advection only transports particles towards increasing values of $z$,
away from the acceleration site. Figure 3 shows a schematic diagram of this model.

When the time reaches $t_{\text{inj}}$ we output the a two-dimensional distribution of the number of
electrons as a function of their energy and position along the $z$ direction. Electrons injected ear-
lier in the simulation time have suffered more cooling losses (Figure S13) and on average are
also farther away from the injection location than freshly injected electrons, due to advection.
With the computed electron distribution as input, the resulting radiation spectrum due to syn-
chrotron and IC emission was calculated using gamera, including Klein-Nishina corrections to
the IC spectrum. Simulated emission profiles for selected photon energies and positions along
<figure><img src="image_19.png" /></figure>

Figure S14: **Calculated deceleration profile for each of the jets.** Shape of the deceleration profile $\Lambda(z)$ as a function of distance from the central binary $z$. The profile is derived from the X-ray data (14, 15) for the western (red line) and eastern (blue line) jets. The curves have been smoothed with a Gaussian kernel with a width 0.07°, approximately the H.E.S.S. PSF.

the $z$ direction were determined. We fitted the predicted flux from the model to the measured gamma-ray profiles with two free parameters: the velocity at the base of the outer jets $v_0 = v_{jet}(0)$ and the diffusion coefficient $D_{100}$. The results are shown in Figure 4.

The only deviation from the observed trend of lower energy gamma rays reaching peak surface brightness further away from the central binary is the lowest-energy bin on the eastern side. In this case, sub-threshold excess emission is only detected outside of the X-ray jets region and close to the outer jet base. However, the 1D model prediction for the spatial distribution of the gamma-ray emission is still compatible with the flux measured at that location within the $1\sigma$ uncertainty.

## Velocity profiles of the jets

The jets are expanding into slightly different ambient conditions: the western jet advances toward the direction of the Galactic plane [(19), their figure 1], leading to an increase in the ambient density. The termination regions of the eastern and western jets lie at a vertical distance from the Galactic plane of approximately 320 and 120 pc, respectively. This probably explains
the observed different jet length on each side (19). The flow velocity along the jet, $v_{jet}(z)$ was characterised as the product of the value of the velocity at the base, $v_0$ and the normalised shape of the velocity profile $\Lambda(z)$, such that $v_{jet}(z) = v_0 \cdot \Lambda(z)$. We constructed the profile $\Lambda(z)$ assuming the jet is cylindrically symmetric and approximately incompressible, meaning that as the jet decelerates, its cross-section increases to preserve the density. Assuming the jet axis to be perpendicular to the line of sight [the measured angle is $\approx 80^\circ$ (8)], the width of the jet can be estimated using the X-ray data (14, 15). The absolute value of the jet width is not required for the model, only its evolution as a function of distance from the base of the outer jet. The velocity profile $\Lambda(z)$ was computed as the inverse square of the width and normalised such that $\Lambda(0) = 1$. The resulting velocity profiles for each of the jets are shown in Figure S14, smoothed with the H.E.S.S. PSF. Variations of this profile on scales smaller than the H.E.S.S. resolution will not affect our results. The best-fitting value of $v_0 = (0.083 \pm 0.026_{\text{stat}})c$ results in a velocity profile $\Lambda(z)$ which is consistent with the upper limit of 0.023c derived around the termination regions furthest from the central binary.

## Systematic uncertainties on the velocity

To assess the impact of the choice of injected electron spectrum parameters and magnetic field on the derived value of $v_0$, we investigated alternative combinations of these parameters. We fixed $\Gamma_e$ and $\log_{10}(E_{cut})$ to values in the ranges 1.6 to 2.4 and 1.8 to 4, respectively, in steps of 0.05 for both parameters. The remaining free parameters (magnetic field and normalisation of the electron spectrum) were fitted to the multi-wavelength data (Figure S10). All the parameter sets defined by the interior of the $2\sigma$ likelihood surface in the $\Gamma_e$ and $\log_{10}(E_{cut})$ space were selected. For each of these electron spectrum parameters, the best-fitting value of $v_0$ was derived by fitting the spatial model prediction to the gamma-ray profiles shown in Figure 4 under the assumption of a decelerating jet. The diffusion coefficient was fixed, otherwise the fit was not
<figure><img src="image_20.png" /><figcaption>Figure S15: **Estimation of systematic uncertainties due to the model parameters.** A: Likelihood profiles of the velocity parameter for each of the considered combinations of magnetic field and electron spectral parameters (grey lines). The black solid line shows the likelihood profile corresponding to the parameters obtained from the fit to the multi-wavelength (indicated as MWL) SED. The black dotted line indicates the best-fitting value of the velocity from the multi-wavelength SED. Black dashed lines mark the likelihood values corresponding to 1, 2 and 3 σ significance levels. The red solid line shows the best-fitting model when only H.E.S.S. data is considered. B: Distribution of the $v_0$ values from all the models. The black dotted line and grey band indicate the best-fitting value and statistical uncertainty of the velocity from the multi-wavelength SED. The blue line shows a Gaussian function fitted to the histogram, where the fitted width is indicated by the blue dashed lines. This value is an estimate of the systematic uncertainty $\sigma_{v_0, syst}$.</figcaption></figure>

stable for some of the combinations. Figure S15A shows the likelihood profiles for each of these models. This approach only uses the H.E.S.S. data, so the parameters from the multi-wavelength model are not necessarily the best-fitting combination, although they lie within $1\sigma$ of the minimum. Figure S15B shows the distribution of the best-fitting values of $v_0$ for the considered models, compared to the result obtained with the parameters derived from the SED fit (Table S5). All the values of $v_0$ derived with different model parameters lie within the statistical uncertainty band, indicating that the choice of parameters does not introduce a large systematic bias.
# Supplementary Text

## Jet deceleration

In our model we have assumed a decelerating jet, as motivated by the apparent observed spread of the X-ray emission. Here we consider alternative possible jet dynamics to test the robustness of our conclusions. We explore the case of a jet moving with constant velocity $\Lambda(z) = 1$ for both constant section (no adiabatic losses) and expanding section (with adiabatic losses $\frac{dE}{dt} = -\frac{1}{3}(\vec{\nabla} \cdot \vec{v}_{jet})E$ (Figure S13). Fitting the value of $v_0$ to the gamma-ray profiles under these assumptions results in $v_0 = (0.045 \pm 0.014_{stat})c$ and $v_0 = (0.061 \pm 0.013_{stat})c$ for the cases with and without adiabatic losses, respectively. These values of $v_0$ are consistent (within the uncertainty) with the value obtained assuming deceleration, $v_0 = (0.083 \pm 0.026_{stat})c$. Any scenario in which $v_0 = 0$ is disfavoured by more than $5\sigma$, so regardless of the jet internal structure, advection of the particles is required to explain the observations. We cannot distinguish between the different jet propagation scenarios with the available observations.

## Diffusion

The value of $D_{100}$ is fitted at the same time as $v_0$ to the gamma-ray data shown in Figure 4. The resulting best-fitting values are $D_{100} = (2.3 \pm 1.4) \cdot 10^{28} \text{ cm}^2 \text{ s}^{-1}$ when assuming a decelerating jet and $D_{100} = (4.7 \pm 4.1) \cdot 10^{27} \text{ cm}^2 \text{ s}^{-1}$ when assuming the velocity to be constant (with constant section). Comparing the advection and diffusion lengths within a cooling time $t_{\text{cool}}$ (Figure S13) we find that $D \ll \frac{t_{\text{cool}}v_0^2}{2}$ at all relevant energies, which indicates that advection is the dominant particle transport process taking place in the outer jets of SS 433. The best-fitting value of $D_{100}$ is approximately an order of magnitude lower than the average Galactic diffusion coefficient (78), which we interpret is due to the stronger magnetic field in SS 433. There is no a priori reason to expect the diffusion coefficient to match the Galactic value, given that the properties of the medium inside W50 are likely determined by the jets themselves and are
unlikely to be similar to the average interstellar medium.

## Magnetic field

The magnetic field strengths $B$ reported in Table S5 were derived from a fit to the multi-wavelength SED of the jets, including both the X-ray and gamma-ray data (Figure S10). The energy-dependent morphology provides an alternative way to estimate the average magnetic field in the outer jets that is independent of the X-ray observations. The dominant processes responsible for the observed morphology are the advection of particles in the jet flow, parameterised by $v_0$ and energy losses due to synchrotron emission, parameterised by $B$. We fit the model gamma-ray profiles to the observed data with $v_0$, $B$ and the normalisation $\alpha$ as free parameters, with the rest of the parameters fixed to the values reported in Table S5. We do this for both the constant velocity and the decelerating flow case. We find that the best-fitting value for the average magnetic field strength in both outer jets is $21.0\pm1.7$ µG and $20.5\pm2.5$ µG for the constant and decelerating flows respectively. These magnetic field values are consistent with those obtained from the SED, implying that the estimate of the average magnetic field reported in Table S5 is consistent with the observed energy-dependent morphology. The values of $v_0$ are consistent with those found for the fixed magnetic field, but with larger statistical uncertainties due to the additional free parameter.

## Contribution from hadronic processes

The observed energy-dependent morphology requires the bulk of the gamma-ray emission to be the result of IC emission from electrons. The cooling time of protons via hadronic interactions with surrounding gas at high proton energies is nearly independent of energy and inversely proportional to the gas density (33). Consequently, protons lose energy slowly for any expected value of gas density in the outer jets (see below) and across all relevant proton energies. In
contrast, the electrons cool rapidly at high energies via inverse Compton and synchrotron losses. While the presence of accelerated hadrons in the jets cannot be ruled out by our observations, we rule out a large contribution of hadronic emission to the observed gamma-ray flux due to the lack of dense target material in the jets region. While the ambient density in the SS 433/W50 complex is unknown, observations (79) and simulations (19) have derived values in the $0.1 \text{ to } 2 \text{ cm}^{-3}$ range. If we attempt to model the H.E.S.S. observations as proton-proton (p-p) interactions using gamera, we find that for values of the ambient density below $3 \text{ cm}^{-3}$, more than 100% of the available jet power needs to be injected into protons to match the observed flux. To reach the observed fluxes while requiring lower fractions (~10%) of the available power, p-p interactions would require a denser ($>20 \text{ cm}^{-3}$) medium. However, the absence of spatial correlation between the gamma-ray emission, especially at the highest energies, and the presence of dense target material (Figure S16) provide an additional indication that the majority of the observed emission is due to relativistic electrons.

While the observed gamma-ray emission is dominated by the inverse Compton emission of electrons, we expect protons or heavier nuclei to be accelerated to similar or greater energies at the same acceleration site as the TeV emitting electrons. Radio observations have detected several molecular clouds in the vicinity of the system (80), some aligned with the SS 433 jet axis (79), mostly in the western extremity of W50 (22, 81, 82). Their kinematics indicate a possible connection with SS 433 (81, 82). However, there are no distance estimates to these clouds, so their position with respect to the jets remains unknown. We do not observe gamma-ray emission from any of these clouds.

## SS 433 as a cosmic-ray source

Microquasars have been proposed as candidate sources of the Galactic cosmic-rays (84–87). The hard injection spectrum and high maximum energy inferred from our and X-ray obser-
<figure><img src="image_21.png" /><figcaption>Figure S16: **Location of possible sites of hadronic interactions.** The H.E.S.S. significance map above 10 TeV (rotated from the orientation in Figure 2C) compared to observed gas locations, which are possible target material for hadronic interactions. Equatorial coordinates are shown for the J2000 equinox. Black contours indicate significances of 4, 5 and 6σ in the H.E.S.S. map. The pink contour indicates Hα + [NII] emission from ionised gas (83), green corresponds to CO observations revealing four molecular clouds N1 to N4 (82) and light blue to neutral hydrogen emission from diffuse neutral gas (79). The ROSAT X-ray contours (14) are shown for reference in white. There is no correlation between any of the potential targets and the H.E.S.S. emission above 10 TeV.</figcaption></figure>
vations are consistent with this proposal, though SS 433 is unlikely to contribute to the local cosmic-ray flux at Earth. Given an assumed upper limit on the system age of $t = 30\,000$ years and adopting the average Galactic diffusion coefficient $D_{gal} = 10^{28} (E/1\text{GeV})^{1/3} \text{cm}^2 \text{s}^{-1}$ (78), the distance that cosmic rays can traverse via diffusive propagation is approximately $r = \sqrt{4D_{gal}t}$. Even for proton energies of 1 PeV, this gives $r \approx 0.6 \cdot (t/30\,000\,\text{yr})^{1/2}$ kpc, much smaller than even the lowest estimates for the distance to SS 433 of 3.8 kpc (88). To reach that distance, the average Galactic diffusion coefficient would need to be several orders of magnitude larger than predicted, or the system would need to be around 40 times older, which is incompatible both with the GeV measurements (Figure S11) and with the highest estimates of the age of W50 of 100 000 yr. (20). Although we expect protons and other nuclei to be accelerated in the jets of SS 433, we conclude that they do not contribute to the cosmic-ray flux measured on Earth.

We nevertheless consider the potential contribution of microquasars to the average Galactic cosmic ray population. Taking our derived value of the magnetic field (B), and the extent of the jet base in X-ray observations (16) as a proxy for the shock width (R), we infer from the Hillas limit (48), a maximum energy $E_{Hillas}$ of

$$E_{Hillas} \approx 10Z \left(\frac{B}{20\mu\text{G}}\right) \left(\frac{u_1}{0.26c}\right) \left(\frac{R}{1.6\text{pc}}\right) \text{PeV}, \quad (S2)$$

where Z is the atomic number. While this exceeds the maximum electron energy discussed above (Figure S13), protons and other nuclei experience less radiative losses. Systems similar to SS 433, should they exist, provide intermittent contributions to the Galactic cosmic-ray budget at a few PeV. If similar acceleration occurs in extra-galactic jets on larger scales, they could reach the EeV regime of ultra-high-energy cosmic rays.