A jet bent by a stellar wind in the black hole X-ray binary Cygnus X-1

Understanding the impact of accreting supermassive black holes on the evolution of galaxies and cosmic structures is one of the key motivators for studying relativistic jets (for example, refs. ^1,2,3). Jets are observed to drive large-scale shocks, pollute interstellar gas with magnetic fields⁴ and cosmic rays, generate large-scale turbulence, and evacuate large-scale cavities of gas on scales of galaxy groups and clusters^5,6. A fundamental difficulty in modelling the kinetic feedback comes from the lack of an instantaneous jet power measurement, which would inform us of the fraction of the accreted energy that is converted into the kinetic energy of the jet. Hence, we have had to rely on time-averaged jet power estimates from the calorimetry of jet-inflated bubbles to constrain this feedback, which involves averaging the total kinetic power output by a jet over its lifetime, which far exceeds the variability timescale of the accretion flow. Owing to this mismatch in timescales, calorimetric measurements cannot be used to accurately calibrate the instantaneous kinetic feedback efficiency of accreting black holes, which is a key input for models of large-scale structure formation (Methods). These reasons underscore the need to develop a method that measures the instantaneous jet power of an accreting black hole.

High-mass black hole X-ray binaries provide a unique opportunity to measure the instantaneous jet power in an individual system via the theoretically predicted interactions between the stellar wind and the jets^7,8,9, although these have not to date been observationally confirmed. In these systems, the black hole accretes from the strong stellar wind of its massive companion star. The jets launched by the black hole must then propagate outwards through that wind. The impact of the wind can bend the jet away from the companion star^7,9, which, when combined with the orbital motion of the black hole, results in a helical jet outflow⁸. The overall jet trajectory is set by the relative strengths of the wind momentum flux and the momentum flux of the jet. Hence, if the wind parameters are known, instantaneous measurements of the jet power, speed, geometry and any misalignment between the jet and the binary can be made.

At a distance of 2.22 kpc (ref. ¹⁰), the high-mass black hole X-ray binary Cygnus X-1 contains a 21.2 ± 2.2 M_⊙ black hole¹⁰ in a 5.6-day binary orbit¹¹ with a supergiant companion of spectral type O and mass \(40.{6}_{-7.1}^{+7.7}\,{M}_{\odot }\) (ref. ¹⁰). The black hole accretes from the stellar wind of the donor star, whose mass-loss rate is (2.57 ± 0.05) × 10⁻⁶ M_⊙ yr⁻¹ (ref. ¹²). In the hard X-ray spectral state, a steady jet is launched from close to the black hole, which can be resolved by high-angular-resolution radio observations with very long baseline interferometry (VLBI)¹³. Early VLBI observations of Cygnus X-1 performed in the late 1990s detected only the Doppler-boosted approaching jet extending 10–15 mas away from the bright core at a position angle of approximately −26° east of north^13,14. A more sensitive VLBI campaign performed in 2016 detected the corresponding Doppler-deboosted receding jet¹⁰ for the first time. The receding jet was found to share a similar position angle to the approaching jet on the plane of the sky. The small, measured eccentricity of the binary (e = 0.019 ± 0.003; ref. ¹⁰) and low space velocity with respect to its suggested birthplace in the Cygnus OB3 association¹⁵ imply a relatively small natal kick at black hole formation, such that the jet axis should be relatively well aligned with the binary orbit. However, recent works have suggested a substantial misalignment between the jet axis and the orbital angular momentum vector^16,17. Moreover, the calorimetry of Cygnus X-1 (ref. ^18,19,20) has been used to anchor many observationally derived black hole scaling relations (such as jet power versus radio luminosity²¹) across the entire mass range of black holes (using the scale-invariant nature of black holes^22,23,24), thus further highlighting the need to measure the instantaneous jet power of Cygnus X-1.

The 2016 VLBI campaign undertook six 8.4-GHz observations with the Very Long Baseline Array (VLBA) and three 5-GHz observations with the European VLBI Network (EVN) to monitor the hard-state jet of Cygnus X-1 around an entire 5.6-day binary orbit. We reanalysed these data, detecting both the approaching and receding jets in each observation. Using a stack of the VLBA images, we found a brightness ratio in the range 1.5–5.5 (Methods and Extended Data Fig. 1). Stacking archival VLBA observations made in 1998 and 2009 also recovered the receding jet, at similar brightness ratios, thus showing that the jets are intrinsically two sided.

Our analysis of individual VLBA and EVN images from the 2016 campaign showed that the observed position angles of the approaching and receding jets varied with orbital phase (Fig. 1). We verified with archival data (Extended Data Table 1 and Supplementary Tables 1 and 2) that the variation in the position angle of the approaching jet was reproducible, with orbital-phase-dependent deviations from a median position angle that was constant over an 18-year period (Extended Data Figs. 2, 3 and 4). Moreover, a model-independent analysis of the jet (Methods) showed that the approaching and receding jets bend in different directions, being deflected away from the position of the donor star (Fig. 2).

Full size image

a–i, Individual VLBI images (with the six 8.4-GHz VLBA observations in blue and the three 5-GHz EVN observations in purple), rotated anticlockwise by 25°, and with an asymmetric axis ratio to help visualize the jet bending. The orbit of the donor star (scaled up by a factor of 30 for the VLBA and 50 for the EVN) is shown in each image in red. The intensity of the radio jet images is shown in a logarithmic scale to visualize a wide dynamic range of emissions. The solid red circle indicates the star’s position at the midpoint of the observation. Each panel is for a different orbital phase ϕ of the observation: ϕ = 0.09 (a); ϕ = 0.27 (b); ϕ = 0.45 (c); ϕ = 0.56 (d); ϕ = 0.63 (e); ϕ = 0.75 (f); ϕ = 0.80 (g); ϕ = 0.92 (h); ϕ = 0.98 (i). j,k, The median position angle of the approaching jet (j) and receding jet (k) as measured from the core along with their 1σ errors. The bending seems to lag by quarter of an orbit, due to the finite time taken by the bent jet to travel downstream. l, The black hole orbit projected onto the plane of the sky, along with the phase coverage of each observation as shaded arcs in blue (VLBA) and purple (EVN). Their radii have been increased with time for clarity. All errors are shown at the 1σ level. PA, position angle.

Full size image

a,b, VLBA images of Cygnus X-1 close to superior conjunction (a) and inferior conjunction (b), when the jets showed significant deviations from the median position angle. The grey ellipse in each image indicates the synthesized beam size. The jet brightness profile was measured along lines of constant declination spaced by 1 mas at the locations indicated by the horizontal lines for the core (dotted), approaching jet (solid) and receding jet (dashed). A linear fit to the jet ridge lines is shown in white, demonstrating that the approaching (solid) and receding (dashed) jets are not collinear. c,d, Normalized cross-correlations of the downstream jet brightness profiles indicated in a (c) and b (d), with the brightness profile at the declination of the core. The negative of the lags is shown for the receding jets (dashed curves) to aid comparison with their approaching counterparts (solid curves). As the jets are not oriented north–south, the peak cross-correlation lag increases on moving downstream. Dec., declination; RA, right ascension.

The observed difference in the position angles of the approaching and the receding jets cannot be explained by simple precession of the jet axis, which would predict that the corresponding approaching and receding jet segments are aligned. However, it naturally arises in a scenario in which the jets are bent away from the donor star by the impact of the stellar wind. As the wind density decreases with distance from the star, most of the bending should occur within about an orbital separation of the jet nozzle (~0.1 mas on the plane of the sky), so that we see only the asymptotic bending angle at VLBI scales⁷. This explains the relatively linear nature of the observed approaching and receding jets and indicates that the central binary system is close to where the radio emitting jets meet, at the bright core. This is in agreement with theoretical predictions of the photospheric distance of Cygnus X-1²⁵.

We, therefore, consider an analytical model of the wind-induced jet bending. The model balances the wind momentum flux with the lateral momentum flux of the jet^8,26,27 and accounts for orbital motion to predict the overall helical structure of the jet (Extended Data Figs. 5 and 6). This physically motivated model also accounts for non-ballistic effects as the helical jet pushes against the wind. To infer the jet properties at launch, we fitted this analytical, numerically evaluated model^8,27 to our measured jet structure (Extended Data Fig. 7). Simultaneously fitting all six VLBA epochs from 2016 gave a jet power \({\log }_{10}[{L}_{{\rm{j}}{\rm{e}}{\rm{t}}}\,\)\(({\mathrm{erg}}\,{\mathrm{s}}^{-1})]=36.{9}_{-0.5}^{+0.2}\), a jet speed at launch \(\beta =0.4{7}_{-0.05}^{+0.06}\) and a jet half-opening angle 0.8 ± 0.5° (Extended Data Table 2). The fitted trajectories for each epoch of the 2016 VLBA observations are shown in Fig. 3. Our fit provides a measurement of the instantaneous jet power in an accreting black hole together with a speed measurement for the compact jet, quantities that have previously been difficult to constrain reliably (for example, refs. ^13,28,29). We note that our fit provides the same jet power for the different wind loss rates of the Cygnus X-1 donor star as those in the literature^12,30,31 (Methods and Extended Data Fig. 8).

Full size image

Contour plots of each image are shown. Blue and red markers denote the locations of the point source components required to model the approaching and receding jets, respectively, as determined from fitting the visibility data in the u–v plane. The errors associated with these point source components are smaller than the size of the markers. The purple dashed line shows the median trajectory from the results of fitting the physical model to the point source locations. The cyan shaded area shows a set of 200 random draws from the posterior distribution. All images have been rotated anticlockwise by 25° to align the median jet axis vertically on the figure. An asymmetric axis ratio has been used to better visualize the jet bending. The date and orbital phase of each epoch are indicated in the bottom left corner. Dist., distance.

The black hole in Cygnus X-1 was originally proposed to have formed by direct collapse¹⁵. In the absence of a supernova kick, the black hole spin axis would, therefore, be expected to be well aligned with the orbital angular momentum. However, the recent detection of high fractional X-ray polarization from the system indicates a substantial misalignment of >18° along the line of sight¹⁶. By contrast, a study of the frequency-dependent phase lags observed in orbital-phase-folded radio light curves indicated a misalignment of 20–30° in the plane of the sky¹⁷. Given these contrasting claims, we explored the impact of misalignment on the calculated jet trajectories.

In the presence of a misaligned jet, we would expect highly asymmetric non-ballistic wind–jet interactions to occur near the base of the jet. Regardless of the geometry of the misalignment, this would cause the approaching jet to propagate towards the star (where the wind is denser) at a particular orbital phase and the receding jet to do so half an orbit later. Substantial misalignment (≳10°) would then lead to strong, highly asymmetric jet bending between the approaching and receding jets due to the different non-ballistic forces imparted on them by the wind, which we do not observe (Supplementary Videos 1 and 2).

We refitted our jet trajectories with more parameters to allow for a misalignment between the jet and orbital axes and found a jet power \({\log }_{10}[{L}_{{\rm{j}}{\rm{e}}{\rm{t}}}\,({{\rm{erg}}\,{\rm{s}}}^{-1})]=37.{3}_{-0.2}^{+0.1}\), a somewhat higher jet speed \(\beta =0.6{8}_{-0.13}^{+0.08}\) and a best-fitting misalignment \(5.{2}_{-1.1}^{+1.0}\)°. However, to account for possible systematics due to averaging of the jet trajectories over a 12-h observation (Methods), we adopted a conservative 3σ upper limit of 8.2° on the misalignment between the jet and the binary.

Relaxing the requirement for the jet to be conical (that is, allowing the jet radius to scale with distance downstream as r ∝ z^ϵ) did not change this conclusion, as it gave consistent values for the jet power, jet speed and misalignment angle and favoured a conical jet with \(\epsilon =1.0{2}_{-0.08}^{+0.08}\). Our fitted conical jet geometry for Cygnus X-1 is in agreement with the lack of a strong re-collimation shock expected for our estimated jet power³² (Methods). The absence of substantial jet–orbit misalignment is in agreement with a lack of any observed Lense–Thirring precession of the jet, as indicated by the stable mean position angle of the jet on the plane of the sky over the 18 years of VLBI observations used here (Extended Data Fig. 4).

Such a small jet–orbit misalignment implies that alternative explanations are required for the high fractional X-ray polarization observed by the Imaging X-ray Polarimetry Explorer, such as the presence of a relativistic outflow in the corona³³. The radio phase lags found by ref. ¹⁷ could be explained by the helical jet structure created by the bent jet and the orbital motion of the black hole. The low misalignment is also in agreement with the low eccentricity and peculiar velocity of the system^10,15 and with theoretical expectations for the formation of such massive black holes³⁴.

Although our jet power estimate from the three physical models (without misalignment, with misalignment, and with misalignment and a non-conical jet) are in agreement with each other at the 1σ level and are relatively insensitive to uncertainties in the mass-loss rate of the donor star wind (Methods), we favour the jet power estimate from the model that allows for both misalignment and a non-conical jet geometry, that is, \({\log }_{10}[{L}_{{\rm{j}}{\rm{e}}{\rm{t}}}\,({{\rm{e}}{\rm{r}}{\rm{g}}\,{\rm{s}}}^{-1})]=37.{3}_{-0.2}^{+0.1}\), as it makes the smallest number of assumptions. This implies that over the few million years lifetime of the black hole X-ray binary, the total kinetic feedback from the jets in Cygnus X-1 would be several times ~10⁵⁰ erg, comparable with the kinetic feedback from a supernova. Our measured instantaneous jet power is in excellent agreement with the time-averaged jet power of 4–14 × 10³⁶ ergs s⁻¹ derived for Cygnus X-1 through calorimetry (see Methods for an extended discussion of the calorimetry) and is shown in Fig. 4. The striking similarity between the bolometric X-ray luminosity of Cygnus X-1 in the hard state³⁵ (rescaled to an updated distance of 2.22 kpc; ref. ¹⁰) and our jet power measurement validates the commonly assumed accretion to jet energy conversion fractions assumed in typical galaxy formation simulations^2,3,36,37 (Methods).

Full size image

Probability density functions estimated from posterior samples for the instantaneous jet power for our three different models: model 1, with no jet–orbit misalignment and a conical jet geometry (blue); model 2, which allows for misalignment with a conical jet geometry (magenta); and model 3, which allows for both misalignment and a non-conical geometry (orange). The probability density function of the derived jet power is relatively similar across all three models, and it shows that the instantaneous jet power is consistent with the time-averaged power required to inflate the arcminute-scale radio nebula surrounding the source (shown by the grey shaded regions^18,19,20).

With a robust, instantaneous measurement of the jet power in Cygnus X-1, we validated the use of calorimetry to calibrate the fraction of accretion energy converted into the kinetic energy of the jet. The close agreement between the instantaneous and time-averaged jet powers implies the long-term stability of the jets produced by radiatively inefficient accretion flows (in the hard X-ray spectral state). Although Cygnus X-1 is such a radiatively inefficient system, this validation strengthens confidence in the broader applicability of calorimetric techniques for estimating jet power in other black hole systems, including active galactic nuclei, irrespective of their specific accretion regime. On larger scales, our measurement of the instantaneous accretion to jet energy conversion fraction lends strong support for the assumed energy budget of accreting black holes in large-scale cosmological simulations.

Data reduction

The primary observations considered in this work were obtained using the VLBA and the EVN and covered a single binary orbit of Cygnus X-1 in May–June 2016¹⁰. To probe the long-term evolution of the jets, we also analysed all existing archival data obtained at frequencies between 5 GHz and 8.4 GHz (in the hard state)^13,14,38, which are tabulated in Extended Data Table 1.

All VLBA observations were processed using the Astronomical Image Processing System, version 31DEC22³⁹, following the standard procedures for phase-referenced experiments (http://www.AIPS.nrao.edu/cook.html). Having performed the recommended initial calibration steps, we then exported the Cygnus X-1 data to difmap⁴⁰ for further calibration (phase only followed by joint amplitude and phase self-calibration) and imaging. For the EVN observations, we followed the EVN guide (https://www.evlbi.org/evn-data-reduction-guide). We used the flagging and the delay, band-pass and a priori amplitude calibration solutions generated by the EVN pipeline, followed by global fringe fitting performed with the Astronomical Image Processing System. The Cygnus X-1 data were then exported to difmap for self-calibration and imaging as done for the VLBA data. A montage of all naturally weighted VLBA (EVN) images is shown in Extended Data Fig. 2 (3).

We then used difmap to model the jet as a series of individual components in the visibility domain (as opposed to an image-based fit), as that method is less susceptible to artefacts introduced by the incomplete sampling in the u–v plane. As we expected the core to be compact but resolved parallel to the jet (due to possible ejection of new components), we modelled it using an elliptical Gaussian. The downstream emission was unresolved in the transverse direction, and hence, we modelled it as a series of delta functions. The number of delta components was determined by the extent of the resolved jet, and we iteratively increased the number of delta components until the residual maps (an image of the difference between the modelled visibilities and the observed visibilities) showed no significant emission. The results were not highly sensitive to the precise number of components used, provided the convergence criteria for the fitting were satisfied. Having generated a preliminary model in difmap, we conducted our final model fitting using the UVMULTIFIT package⁴¹ within the Common Astronomy Software Application⁴² to determine the uncertainties associated with each of the fitted parameters.

Model-independent illustration of bent jets in Cygnus X-1

We selected two of our 2016 VLBA epochs taken at different points in the orbit that showed significant jet bending (Fig. 1a,c, respectively). We extracted intensity profiles at constant declination, along the core and every ±1 mas downstream (see top panels of Fig. 2). The downstream intensity profiles were cross-correlated with the intensity profile of the core, such that the lag represents the shift in right ascension of the peak intensity at the relevant declination, relative to the right ascension of the core. We show the lag measurements in the bottom panels of Fig. 2, with the lags measured along the receding jet inverted for ease of comparison. For a collinear jet (with no bending), the lags measured at equal angular separations from the core along the approaching and receding jets would be expected to be equal in magnitude, contrary to what is shown in Fig. 2. For a jet aligned purely north–south, the cross-correlation in every declination slice would peak at zero lag. For any other orientation, the cross-correlation lag would increase linearly on moving downstream. On one side of the orbit (right panels), the lags along the approaching jet are larger in amplitude compared with those measured at the same distance downstream along the receding jet. This implies that the approaching jet is rotated further clockwise on the plane of the sky than the receding jet, confirming that the jets are not collinear. On the other side of the orbit (left panels), the opposite lag behaviour is seen, demonstrating that this bending of the jets varies with orbital phase.

Physical model of the jet bending

The combination of orbital motion, stellar winds and magnetic field strengths can produce very complex jet outflow morphologies in high-mass X-ray binaries^27,43. However, given that our VLBI images lack the resolution to resolve these fine-scale features, we employed a simplified analytical model that captures the average jet flow under the assumption of momentum conservation in the wind–jet interaction^8,26. In the absence of observational evidence to the contrary, we did not assume a dynamically important magnetic field in our model.

Having confirmed the existence of jet bending via the model-independent cross-correlation, we modelled the observed bending by considering the momentum transferred from the wind to the jet^8,26,27. In a non-rotating system, the wind would bend the jet away from the donor star, eventually reaching an asymptotic bending angle⁷. However, in a rotating binary system, the combination of orbital motion and radial bending would result in a helical jet outflow^8,26,27.

We used a framework adapted from ref. ²⁶ to calculate the jet bending, using a co-rotating frame, as shown in Extended Data Fig. 5. The x axis points towards the black hole, the z axis is parallel to the orbital angular momentum vector of the binary and the y axis is orthogonal to both, thus forming a Cartesian coordinate system. We define the initial momentum of the jet along the x, y and z directions as

where L_jet is the jet power, β is the jet speed at launch, Γ is the corresponding Lorentz factor and c is the speed of light. We assume an isotropic wind. The trajectory of the jet downstream is numerically calculated using the equation \({{\bf{p}}}_{i+1}={{\bf{p}}}_{i}+{{\bf{F}}}_{i}\,dt\), where the thrust imparted by the wind (\({{\bf{F}}}_{i}\)) on the jet that is transferred to the (i + 1)th jet segment is calculated using the wind density and wind velocity calculated upstream. Very close to the base of the jet, non-ballistic jet–wind interactions determine the jet structure, and further downstream, due to reduced wind density, a ballistic propagation of the jet is assumed by our model. An example jet trajectory calculated at different orbital phases is shown in Extended Data Fig. 6.

We used dynamical nested sampling⁴⁴ implemented using the dynesty Python package⁴⁵ to fit our model to the individual delta components from modelling the 2016 VLBA observations of Cygnus X-1. As an initial fit, denoted as model 1, we used a minimal subset of free parameters, comprising the jet speed at launch β, the jet power L_jet, the jet half-opening angle near the base ϕ_1/2, the inclination angle of the jet at launch i, the donor star wind velocity v_wind and the jet position angle at launch projected onto the plane of the sky. The terminal velocity of the donor star wind was expected to be within the range 1,600–2,400 km s⁻¹ (refs. ^46,47), and hence, it was left as a free parameter in our fit, with 1,600 km s⁻¹ and 2,400 km s⁻¹ set as the bounds. The priors and posterior distributions of the model fit are tabulated in Extended Data Table 2. We found a jet power \({\log }_{10}[{L}_{{\rm{j}}{\rm{e}}{\rm{t}}}\,({{\rm{e}}{\rm{r}}{\rm{g}}\,{\rm{s}}}^{-1})]=36.9{4}_{-0.49}^{+0.22}\), a jet launch speed \(\beta =0.4{7}_{-0.05}^{+0.06}\), a half-opening angle of the jet 0.8 ± 0.5°, a jet launch position angle of −25.0 ± 0.1° and a wind velocity of ≤1,870 km s⁻¹.

Recent studies have raised the possibility of a misalignment between the jet and the binary orbit in Cygnus X-1, whether along the line of sight¹⁶ or in the plane of the sky^17,48. A misaligned jet would result in orbital-phase-dependent asymmetric bending of the approaching and receding jet, as one of the jets would be launched towards the star and undergo significantly greater bending due to a strong non-ballistic interaction with the denser wind. Hence, we modified our physical model to evaluate the possibility of a misalignment by using two more parameters to describe the initial orientation of the jet at launch, namely the misalignment angle θ and the binary orbital phase az when the misaligned jet points along the line joining the star and the black hole. For an az value of 0 (or 0.5), the jets would be pointed towards (away from) the star at both superior and inferior conjunction (a misalignment along the line of sight direction), and for a value of az = 0.25 (or 0.75), the jets would be pointed towards (away from) the star at quadrature phases (a misalignment in the plane of the sky).

Furthermore, as the jet is evaluated in a co-rotating frame, in the presence of a misalignment, each downstream segment would have been launched at a different orbital phase, resulting in different initial jet momenta along the x and y directions (in the case of no misalignment, the initial momenta of all downstream jet segments are zero along the x and the y directions). Hence, in our numerical model, each downstream jet segment was mapped to the orbital phase at which it was launched, which in turn determined the initial momentum of the jet segment along the (x, y, z) directions. For an assumed jet segment length dl, an orbital period T₀ and a jet speed β, then each downstream jet segment (indexed as n = 0, 1, 2, 3, …, where n is the segment number with n = 0 being the segment right next to the black hole) would have been launched at an earlier orbital phase n × Δϕ, where Δϕ = dl/(βT₀). This offset in phase was then used to determine the initial momentum of the jet segment at launch. If the black hole is, at present, at an orbital phase ϕ₀, then the initial jet momentum for the nth downstream segment in the co-rotating frame is given by

where ϕ = 2π(ϕ₀ − nΔϕ + az), and ϕ₀ is the binary phase of the epoch being considered. Animations of jet trajectories for a misaligned jet are provided in Supplementary Information. Fitting this misaligned jet model (denoted as model 2) to the data points from the 2016 VLBA observations, we found a comparable jet power \({\log }_{10}[{L}_{{\rm{j}}{\rm{e}}{\rm{t}}}\,{(\mathrm{erg}}\,{\mathrm{s}}^{-1})]\)\(=37.{28}_{-0.23}^{+0.13}\), a jet speed at launch \(\beta =0.6{8}_{-0.11}^{+0.07}\), a half-opening angle of 1.8° ± 0.7°, a misalignment of \(5.{2}_{-1.1}^{+1.0}\)° with a phase of \(az=0.3{4}_{-0.02}^{+0.03}\) (approximately halfway between an orientation along the line of sight and in the plane of the sky) and a wind velocity of ≤1,980 km s⁻¹. The adopted priors and resulting posteriors are tabulated as model 2 in Extended Data Table 2. We also note that model 2 finds a mean jet position angle of \(-27.{8}_{-0.4}^{+0.5}\)°. We expect this to be the position angle of the jet at launch, and when combined with the asymmetric jet bending from the 5.2° jet–orbit misalignment, results in the apparent large-scale VLBA jet appearing to precess about a mean position angle of −25°, in agreement with the value from the scenario with no misalignment considered in model 1.

Our final model (model 3) incorporated non-conical jet geometries into model 2. We allowed the jet cross section to scale with distance downstream as r ∝ z^ϵ and fitted for the parameter ϵ (where ϵ = 1 for a conical jet). The fitting found the jet power to be consistent with values from the previous models, and ϵ = 1.0 ± 0.1, which is consistent with a conical geometry within uncertainties. The fitting parameters are tabulated as model 3 in Extended Data Table 2. The median and 200 random draws from the posterior distributions from models 1, 2 and 3 are shown overlaid on top of the six VLBA images from 2016 in the top, middle and bottom rows of Extended Data Fig. 7.

Misalignment

The presence of a small but non-zero misalignment of \(5.{2}_{-1.1}^{+1.0}\)° is consistent with the maximum plausible misalignment of ~10° inferred from the ratio of the kick velocity of the system and its pre-supernova orbital velocity¹⁰. The small but non-zero binary eccentricity and space velocity would also be consistent with a small but non-zero misalignment. However, we note the 12-h duration of the individual VLBA epochs, which would have averaged the jet trace over a non-negligible fraction of the 5.6-day binary orbit, could have affected the observed trajectories of our VLBI images. We, therefore, retain a degree of caution in interpreting the misalignment fit and prefer to draw a more conservative conclusion, ruling out misalignments >8°. This is far smaller than the 20–30° misalignments inferred by refs. ^16,17, and in agreement with the stable mean position angle of the jet observed during the 18 years of archival data (Extended Data Fig. 4).

Re-collimation shocks

For sufficiently low jet power (below a critical value) in high-mass X-ray binaries, a lateral re-collimation shock is expected to be driven into the jet by the wind³², resulting in non-conical jet geometry. For a given mass-loss rate \({\dot{M}}_{{\rm{w}}}\), wind velocity v_w and jet velocity v_j of the donor star, the critical power is given by \({P}_{{\rm{c}}{\rm{r}}}={\dot{M}}_{w}{v}_{{\rm{w}}}{v}_{{\rm{j}}}/16\) (ref. ³²). Using the binary parameters adopted in this paper, we found the critical jet power to be (2.5–5.9) × 10³⁷ erg s⁻¹, comparable with our estimate of instantaneous jet power. Thus, our determined conical jet (lack of a strong re-collimation shock) from model 3 is in perfect agreement with theoretical expectations.

Brightness and spectral index of the jet

The brightness ratio of the approaching and receding jets, combined with the spectral index, constrains the intrinsic jet speed^49,50. From the VLBA and EVN montage (Extended Data Figs. 2 and 3), we note that the receding jet was detected only in epochs with noise levels below 0.3 mJy per beam. To accurately measure the brightness ratio and confirm that the non-detection of the receding jet is due solely to noise, we created stacked images using temporally close observations with similar set-ups (for example, bandwidth). After correcting for parallax and orbital and proper motion shifts (to align the centroid of the individual images), we produced three stacks: one from 1998 (BS060AX, BS060BX and BS060CX), one from 2009 (BR141A, BR141B, BR141C and BR141D), and one from 2016 (BM429B–BM429G), shown in Extended Data Fig. 1 (left panels). Stacking revealed the receding jet in epochs where it was previously undetected.

With the receding jet confirmed as a persistent feature over ~2 decades, we measured the brightness ratio in the three stacked images (1998, 2009 and 2016). We extracted spatially resolved intensity profiles every 0.5 mas along both jets—approximately half the restoring beam size—and fitted one-dimensional Gaussians to these profiles. The measured intensities and brightness ratio are presented in the right and bottom panels of Extended Data Fig. 1. Although stacking substantially improved the signal-to-noise ratio, which allowed the previously undetectable counter-jet to be recovered in several epochs, the flux densities of the downstream emission may still be slightly affected by the differences in orbital-phase-dependent bending across the individual epochs. To mitigate this, we adopted a conservative range 1.5–5.5 for the brightness ratio, which we used to constrain the jet speed in section ‘Hard-state jet speed of Cygnus X-1’. We note, however, that the jet speed derived from our physically motivated model is our primary jet speed measurement, and we used the brightness-ratio-based estimate only as a consistency check.

During the 2016 VLBI observations, which monitored a full binary orbit (May–June 2016), the Karl G. Jansky Very Large Array (VLA) simultaneously observed Cygnus X-1 under programme VLA/15B-236. We reduced the data following the VLA guide published by the National Radio Astronomy Observatory (https://casaguides.nrao.edu/index.php/VLA_Continuum_Tutorial_3C391-CASA5.4.0). The spectral index (defined as S_ν ∝ ν^α, where S_ν is the flux density measured at a given frequency ν), derived from flux densities at 5.25 GHz and 7.45 GHz for each epoch, is listed in Supplementary Table 1. As the spectral index varies with orbital phase due to changes in path length through the stellar wind, we used the measurement closest to inferior conjunction (1 June 2016) as indicative of the intrinsic spectral index of the jet.

Hard-state jet speed of Cygnus X-1

The brightness ratio of intrinsically symmetric jets at equal angular separations from the core constrains the product \(\beta \,\cos \,i\), where β is the jet speed and i is the inclination of the system^49,50. For a continuous jet, as seen in Cygnus X-1, the jet speed depends on the measured brightness ratio (S_ratio):

where α is the spectral index (S_ν ∝ ν^α). Given the binary inclination of 27.5° ± 5° (ref. ¹⁰) and accounting for a possible misalignment <8°, we used a brightness ratio of 1.5–5.5 and a spectral index of 0.055 ± 0.016 to constrain the hard-state jet speed from the brightness ratio to be β = 0.27–0.46 (68% high-density interval).

Irrespective of the model used, we found the hard-state jet speed of Cygnus X-1 to be only mildly relativistic, with β < 0.76 within the 68% credible intervals. The physically motivated model found a jet speed of \(\beta =0.4{7}_{-0.05}^{+0.06}\) for a jet aligned with the binary (model 1), and \(\beta =0.6{8}_{-0.13}^{+0.08}\) in the presence of a misalignment (model 2). A rerun of the radio timing study presented in ref. ²⁹, incorporating updated constraints on the photosphere distance¹⁷, found β = 0.42–0.64 (68% high-density interval), providing further confidence in the speeds obtained from our physical model.

The jet speed determined from the brightness ratios is in marginally better agreement with the fitting results from physical model 1 than from models 2 and 3. This could be taken to favour an aligned jet in Cygnus X-1, but uncertainties associated with the brightness-ratio measurement (such as image-based convolution effects in the VLBI images, physical effects such as mass loading or a radio emission contribution from a slower jet sheath where much of the wind–jet interaction occurs) limits its discriminative power. However, such topics remain beyond the scope of the current paper and are impossible to probe observationally with the current data.

Regardless of the presence or absence of a misalignment in Cygnus X-1, based on our physical modelling, we found the jet speed to lie in the range β = 0.42–0.76. A more precise determination of the jet speed in Cygnus X-1 would require higher-precision measurements of the jet trajectory, on both smaller and larger angular scales.

Measurements and applications of jet power

Measurements of jet power are key input for theories of accretion as well as models of galaxy formation and evolution, as both require knowledge of the partitioning of accretion energy between advection, radiation, and jet and kinetic power. However, due to the featureless non-thermal synchrotron spectra, direct measurements of jet power from observations of jets are difficult. Neither the jet velocity nor composition can be determined directly from observations of the synchrotron emission alone, and jet power estimates therefore rely on ad hoc assumptions about the composition and the equipartition fraction of the plasma, as well as the bulk Lorentz factor. These assumptions introduce uncertainties of at least an order of magnitude.

Observational constraints on the time-averaged jet power can be derived by modelling the interactions between the jets and the surrounding interstellar medium^51,52.

The gold standard for measuring jet power to date relies on calorimetric methods, whereby the interaction of jets with some known target material is used to determine the average jet power (for example, refs. ^51,52). The most common method uses observations of jet-inflated cavities in galaxies, groups and clusters of galaxies⁶ (this method was also used when deriving the jet power required to inflate the nebula observed around Cygnus X-1¹⁸). These estimates have been used to derive scaling relations between the radiative and the kinetic power of large samples of jets^53,54.

Because these relations rely on measurements of the kinetic power averaged over the dynamical timescale of the cavity, they introduce large uncertainties when used to estimate instantaneous jet powers, given the generally rapid variability in accretion and, therefore, jet power inferred by observations of X-ray binaries, blazars and nearby jets.

Because of these large uncertainties, numerical models of structure formation generally leave the jet power as a free parameter \({\eta }_{{\rm{j}}{\rm{e}}{\rm{t}}}={P}_{{\rm{j}}{\rm{e}}{\rm{t}}}/\dot{M}{c}^{2}\approx 0.1\) (refs. ^{3,36,37,55,56,57}) that is then tuned so that the simulations reproduce the observed properties of galaxies. Although this assumption appears reasonable (given its order-of-magnitude agreement with the maximum radiative luminosity achievable for radiatively efficient accretion⁵⁸), it could be substantially strengthened by direct observational verification of the value of η_jet, which underscores the importance of the measurement of the instantaneous jet power of Cygnus X-1 presented here, which is broadly consistent with the findings of the aggregate jet power needed in cosmological simulations.

Calorimetry of Cygnus X-1

Of the four calorimetric studies of Cygnus X-1 in the literature^18,19,20,59, the average jet power measurement in the most recent study using MeerKAT data⁵⁹ differs by over an order of magnitude. We attribute this discrepancy to the difference in shock velocities considered in the different studies. The MeerKAT study considers a wide range of shock velocities between 21 km s⁻¹ and 364 km s⁻¹ using constraints from the observed radio emission and upper limits from a lack of observed X-ray emission. Reference ¹⁹ further rules out shock velocities <100 km s⁻¹ using the ratios of the observed O [iii] and Hα emission lines. Using constraints from radio, upper limits from X-ray and a series of optical emission lines, ref. ²⁰ places a very tight constraint on the shock velocity to \(22{0}_{-30}^{+50}\) km s⁻¹. As the jet power in calorimetric models scales as the third power of the shock velocity, we attribute the large discrepancy between ref. ⁵⁹ and the other studies to be due to their consideration of very small shock velocities. Owing to this reason, the average jet power estimates from ref. ⁵⁹ are not discussed in the main text.

Donor star mass-loss rate

Recent results^30,31 have suggested a higher uncertainty in the donor star mass-loss rate than inferred from the originally accepted value¹². We, therefore, tested the sensitivity of our results to the assumed wind mass-loss rate of the donor star by rerunning our model 1 fit for mass-loss rates in the range 0.5–7 × 10⁻⁶ M_⊙ yr⁻¹ (refs. ^30,31). As shown in Extended Data Fig. 8, we found that only the jet half-opening angle varies significantly to match the momentum rate of the wind, with the fitted opening angle increasing at lower mass-loss rates and decreasing at higher rates. This demonstrates that our inferred jet power is relatively insensitive to the assumed wind mass-loss rate from the donor star.

Jet composition of Cygnus X-1

Our physical model of the jet bending makes no assumptions about the jet composition but rather infers the jet power from the jet momentum, which is set by the relative balance between the jet momentum and the wind momentum required to produce the observed jet bending. However, we found a leptonic jet to be highly unlikely in Cygnus X-1 due to the large number of pairs required to match the mass flux of the hadronic jet. Moreover, for a high-mass X-ray binary, a leptonic jet at the base becomes mass-loaded by the wind and becomes baryon-dominated on the observed VLBI scales.

Wind clumping

Although the degree of wind clumping in Cygnus X-1 has not been quantitatively constrained in the literature, we discuss its potential implications for our jet power estimate. A clumpy wind could lead to an overestimation of the average wind density, which implies that the true wind momentum flux may be lower than assumed in this study. Consequently, a smaller jet momentum would be required to reproduce the observed bending angles. Such a reduction in jet momentum could be achieved either by a lower jet power or by a higher jet velocity (equation (3)). Conversely, ref. ⁶⁰ suggests that wind clumping can enhance jet instabilities and potentially disrupt the flow near its base, although the jet may subsequently re-collimate further downstream, a scenario that cannot be ruled out for Cygnus X-1. Nevertheless, given the small observed bending angle, such effects from dense wind clumping are unlikely to play a dominant role in this system.