A rapidly-changing jet orientation in the stellar-mass black hole V404 Cygni

TL;DR

This study observes rapid jet orientation changes in the black hole V404 Cygni, suggesting that Lense-Thirring precession of the accretion disk influences jet dynamics during super-Eddington accretion.

Contribution

It provides observational evidence linking jet precession to accretion disk dynamics in a stellar-mass black hole, highlighting the role of Lense-Thirring precession.

Findings

Jet orientation changes on minutes to hours timescale.

Precession modeled as Lense-Thirring effect on a slim disk.

Implications for jet launching and feedback in accreting black holes.

Abstract

Powerful relativistic jets are one of the main ways in which accreting black holes provide kinetic feedback to their surroundings. Jets launched from or redirected by the accretion flow that powers them should be affected by the dynamics of the flow, which in accreting stellar-mass black holes has shown increasing evidence for precession due to frame dragging effects that occur when the black hole spin axis is misaligned with the orbital plane of its companion star. Recently, theoretical simulations have suggested that the jets can exert an additional torque on the accretion flow, although the full interplay between the dynamics of the accretion flow and the launching of the jets is not yet understood. Here we report a rapidly changing jet orientation on a timescale of minutes to hours in the black hole X-ray binary V404 Cygni, detected with very long baseline interferometry during the peak of its 2015 outburst. We show that this can be modelled as Lense-Thirring precession of a vertically-extended slim disk that arises from the super-Eddington accretion rate. Our findings suggest that the dynamics of the precessing inner accretion disk could play a role in either directly launching or redirecting the jets within the inner few hundred gravitational radii. Similar dynamics should be expected in any strongly-accreting black hole whose spin is misaligned with the inflowing gas, both affecting the observational characteristics of the jets, and distributing the black hole feedback more uniformly over the surrounding environment.

Figures

Equations

$${\mathrm{RA}}_{ik}={\mu }_{\mathrm{ra},k}(t_i-t_{\mathrm{ej},k})+{\mathrm{J}}_{\mathrm{ra},i},\mathrm{and}$$

$${\mathrm{Dec}}_{ik}={\mu }_{\mathrm{dec},k}(t_i-t_{\mathrm{ej},k})+{\mathrm{J}}_{\mathrm{dec},i},$$

$${\mu }_{\begin{array}{c}\mathrm{app}\\ \mathrm{rec}\end{array}}$$

$$\beta \cos \theta$$

$$\tan \theta$$

$$P=\frac{\pi }{3a_{\ast }}\frac{GM}{c^3}{r}_{\mathrm{sph}}^3\left[\frac{1-{\left(r_{\mathrm{in}}/r_{\mathrm{sph}}\right)}^3}{\ln \left(r_{\mathrm{sph}}/r_{\mathrm{in}}\right)}\right],$$

$$\frac{r_{\mathrm{sph}}/r_{\mathrm{in}}}{\dot{m}}\approx 1.34-0.4{ϵ}_{\mathrm{w}}+0.1{ϵ}_{\mathrm{w}}^2-(1.1-0.7{ϵ}_{\mathrm{w}}){\dot{m}}^{-2/3},$$

$$E_{\min }\approx 8\times 10^6{\eta }^{4/7}{\left(\frac{V}{{\mathrm{cm}}^3}\right)}^{3/7}{\left(\frac{\nu }{\mathrm{Hz}}\right)}^{2/7}{\left(\frac{L_{\nu }}{{\mathrm{erg}\mathrm{s}}^{-1}{\mathrm{Hz}}^{-1}}\right)}^{4/7}\mathrm{erg},$$

Full text

{affiliations}

International Centre for Radio Astronomy Research – Curtin University, GPO Box U1987, Perth, WA 6845, Australia

Department of Physics, University of Alberta, 4-181 CCIS, Edmonton, AB T6G 2E1, Canada

East Asian Observatory, 660 N. A’ohoku Place, University Park, Hilo, Hawaii 96720, USA

School of Physics & Astronomy, University of Southampton, Southampton SO17 1BJ, United Kingdom

INAF – Osservatorio Astronomico di Brera, Via E. Bianchi 46, I-23807 Merate (LC), Italy

Astrophysics, Department of Physics, University of Oxford, Keble Road, Oxford OX1 3RH, UK

SRON, Netherlands Institute for Space Research, Sorbonnelaan 2, 3584 CA Utrecht, the Netherlands

Department of Astrophysics/IMAPP, Radboud University, Nijmegen, PO Box 9010, 6500 GL Nijmegen, the Netherlands

Universities Space Research Association, 7178 Columbia Gateway Dr, Columbia, MD 21046, USA

National Science Foundation, 2415 Eisenhower Ave, Alexandria, VA 22314, USA

Department of Physics & Astronomy, Wheaton College, Norton, MA 02766, USA

Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, the Netherlands

Gravitation Astroparticle Physics Amsterdam (GRAPPA) Institute, Science Park 904, 1098 XH Amsterdam, the Netherlands

ESAC/ESA, XMM-Newton Science Operations Centre, Camino Bajo del Castillo s/n, Urb. Villafranca del Castillo, 28692, Villanueva de la Cañada, Madrid, Spain

Institute of Cosmos Sciences, University of Barcelona, Martí i Franquès 1, 08028 Barcelona, Spain

NRAO, P.O. Box O, Socorro, NM 87801, USA

Caltech, 1200 E. California Blvd., MC 249-17, Pasadena, CA 91125, USA

Herzberg Astronomy and Astrophysics Research Centre, 717 White Lake Road, Penticton, BC V2A 6J9, Canada

New York University Abu Dhabi, P.O. Box 129188, Abu Dhabi, United Arab Emirates

Department of Astronomy, University of Virginia, 530 McCormick Road, Charlottesville, VA, 22903, USA

School of Astronomy and Space Sciences, University of the Chinese Academy of Sciences, Beijing 100049, China

Sydney Institute for Astronomy, School of Physics A28, The University of Sydney, Sydney, NSW 2006, Australia

Institute for Space Sciences, Atomistilor 409, PO Box MG-23, 077125 Bucharest-Magurele, Romania

During the 2015 outburst\citemainrodriguez15 of the black hole X-ray binary system V404 Cygni\citemainshahbaz94, we conducted high-angular resolution radio monitoring with the Very Long Baseline Array (VLBA). Our observations (Extended Data Table 1) spatially resolved the jets in this system, on size scales of up to 5 milliarcseconds (12 a.u. at the known distance of $2.39\pm 0.14$ kpc\citemainmiller-jones09; see examples in Figure 1). These jets evolved in both morphology and brightness on timescales of minutes.

The orientation of the jets on the plane of the sky varied between epochs, ranging between $-30.6$° and $+5.6$° east of north (Figure 1, 2, and Extended Data Table 2). This range encompasses the orientation inferred from the position angle of the linearly-polarised radio emission\citemaincorbel00 measured during the 1989 outburst ($-16\pm 6$° east of north; we state all uncertainties at 68% confidence)\citemainhan92. Moreover, during a period of intense radio and sub-millimetre flaring on June 22nd\citemaintetarenko17, we observed multiple ejection events spanning a similar range of orientations over a single four-hour observation (Figure 1), implying extremely rapid changes in the jet axis.

The time-resolved images from June 22nd (see Supplementary Video) show a series of ballistically-moving ejecta that persist for tens of minutes before fading below the detection threshold of $\approx 10$ mJy. The radio emission is dominated by a stationary core that is always present, allowing us to perform relative astrometry on the ejecta. The ejecta appear on both sides of the core, with proper motions ranging from 4.3 to 46.2 milliarcseconds (mas) day*-1* (0.06–0.64$c$ in projection; Figure 3), at position angles between $-28.6$° and $-0.23$° east of north on the plane of the sky (Extended Data Figures 1–4; Extended Data Table 3).

Under the (standard) assumption of intrinsic symmetry, then with the known distance\citemainmiller-jones09 we can use the measured proper motions of corresponding pairs of approaching and receding ejecta to determine $\theta$, the inclination angle to the line of sight, as well as the dimensionless jet speed $\beta=v/c$ (see Methods). We identify three likely pairs of ejecta with consistent position angles and ejection times (denoted N2/S2, N3/S3 and N6/S6; see Figure 3 and Extended Data Figures 1–3), although since their flux density evolution cannot be fully explained by Doppler boosting of intrinsically symmetric jets (see Methods), the assumption of symmetry remains unverified. From these three pairs we determine ($\beta=0.32\pm 0.02$, $\theta=40.6\pm 2.4$°), ($\beta=0.35\pm 0.01$, $\theta=32.5\pm 1.6$°), and ($\beta=0.48\pm 0.01$, $\theta=14.0\pm 0.8$°), respectively (Figure 4). In all three cases the northern component is the faster-moving, and must therefore be the approaching component. For unpaired ejecta, we can use the known distance to solve for $\beta\cos\theta$, subject to an assumption on whether the components are approaching or receding (Figure 4). Again, we find that the jet speed or inclination angle, or both, must vary between ejection events.

The most natural interpretation for changes in jet orientation is precession, as best studied in the persistent X-ray binary SS 433. However, each individual jet component only samples the orientation of the jet axis at the time of ejection. With only twelve discrete components on June 22, we do not have sufficient sampling to determine whether the precession is regular. Our best constraint on the precession period comes from the $\sim 30$° swing in position angle between ejecta pairs N2/S2 and N6/S6, which were ejected only 1.3 hours apart. This places an upper limit of 2.6 hours on the period, although the varying position angles of the intervening ejecta suggest that the true period is significantly shorter. The lower limit of order $\approx 1$ second is set by the lack of any blurring motion of the point source components over the timescale on which they are ejected ($>0.1$s; see Methods). Regardless, since the distribution of position angles for a precessing jet will peak at the two extremes, we can infer a precession cone half opening angle of $\sim 18$° (Figure 2).

Since V404 Cygni likely received a natal supernova kick\citemainmiller-jones09b, a misalignment between the binary orbital plane and the black hole spin is expected. Plasma out of the black hole equatorial plane should then undergo Lense-Thirring precession\citemainlense18, potentially affected by torques from strong magnetic fields and associated jets\citemainliska18. This phenomenon has been proposed to explain the low frequency quasi-periodic oscillations (QPOs) observed at sub-Eddington accretion rates in many X-ray binary systems\citemainstella98,ingram16. Regardless, both theoretical predictions and magnetohydrodynamic simulations\citemainfragile07 of tilted disks have shown that a sufficiently geometrically thick disk\citemainpapaloizou95 can precess as a solid body. To enable communication of the warp, the precession timescale must exceed the azimuthal sound crossing time of the disk. The viscosity and magnetic fields should also be sufficiently low that the disk will not realign within a precession cycle\citemainmotta18.

During its 2015 outburst, the X-ray behaviour of V404 Cygni could be explained by invoking a geometrically thick slim disk configuration\citemainmotta17. The mass accretion rate inferred from the peak X-ray luminosity implies a spherisation (outer) radius for the slim disk consistent with the maximum for solid body precession set by the viscous alignment timescale (see Methods). This makes Lense-Thirring precession a plausible scenario for varying the disk orientation. Precession of the inner slim disk would naturally result in precession of the jets, whether due to the magnetic field lines anchored in the precessing disk, or to realignment of spin-powered jets, either by powerful outflows from the inner disk\citemainbegelman06 or by the precessing slim disk itself\citemainliska18.

While the maximum radiative luminosity detected in the outburst was twice the Eddington luminosity\citemainmotta17, super-Eddington accretion flows are known to drive powerful winds that can carry away a large fraction of the mass flowing in from the outer disk\citemainpoutanen07, implying an outer accretion rate well above Eddington. For moderate spins, mass inflow rates up to a few tens of times the Eddington accretion rate would imply precession periods\citemainfragile07 of up to a few minutes and spherisation radii of a few tens to hundreds of gravitational radii (Extended Data Figure 5). While such short periods would require the jet ejecta to be launched on timescales no longer than a few seconds, they would not require the jets to exceed the Eddington luminosity over the launching timescale (see Methods). The precessing jets could also give rise to optical or infrared QPOs in the optically-thin synchrotron emission from the jet base.

A precessing accretion flow is also consistent with the marginal detections of short-lived low-frequency X-ray QPOs reported at 18 mHz on June 22nd\citemainhuppenkothen17. However, the link between the QPOs and the precessing disk is not clear and their short-lived nature would argue against long-term stable precession. In such a case, the changing mass accretion rate (and hence spherisation radius) would cause bursts of precession, subsequently damped by either disk alignment, or by changes in the sound speed\citemainliska18,motta18. However, Figure 2 shows that the jet axis continues to vary over our full 2-week VLBA campaign. This suggests that precession continues with a relatively consistent cone opening angle, even if the precession timescale varies.

We have observed short-timescale changes in jet orientation from a black hole accreting near the Eddington rate, likely from a reservoir whose angular momentum is misaligned with the black hole spin. This spin-orbit misalignment in a low-mass X-ray binary suggests that the impact of black hole natal kicks can persist even after an evolutionary phase of accretion, and could therefore affect the observed gravitational waveforms\citemainapostolatos94 during black hole merger events arising from the evolution of isolated binary systems.

Our findings are consistent with results from recent relativistic magnetohydrodynamic simulations, which demonstrated (albeit in the absence of radiation pressure) that the accretion flow and jets precess together, due to the combination of Lense-Thirring and pressure or magnetic torques from the inflow/outflow system\citemainliska18. The presence of a rapidly-precessing jet in a high-accretion rate source implies that varying jet inclination angles likely need to be accounted for when interpreting observations of systems such as ultraluminous X-ray sources\citemainmiddleton18, black hole-neutron star mergers\citemainstone13, gamma-ray bursts, tidal disruption events\citemainlei13, and rapidly-accreting quasars in the early Universe.

Kinetic feedback from precessing jets or uncollimated winds in AGN that distribute energy over large solid angles\citemainfalceta-goncalves10 has been invoked to prevent the onset of cooling flows in cool core clusters\citemainvernaleo06 and to solve discrepancies between observed galactic properties and cosmological simulations\citemainweinberger17. For some low-luminosity AGN, which should host geometrically thick accretion flows, light curve periodicities and helical trajectories of jet components have been suggested as direct evidence of jet precession, typically attributed to the presence of a binary supermassive black hole\citemaincaproni04. However, Lense-Thirring precession can also match the observed timescales\citemainnagai10,britzen18, and might be expected in chaotic accretion scenarios. Therefore, as demonstrated by our findings, precessing jets need not always signify binary black holes.

\bibliographystylemain

naturemag \bibliographymainjmj_nature

{addendum}

is linked to the online version of the paper at www.nature.com/nature.

The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. JCAM-J is the recipient of an Australian Research Council Future Fellowship (FT140101082). AJT is supported by an Natural Sciences and Engineering Research Council of Canada (NSERC) Post-Graduate Doctoral Scholarship (PGSD2-490318-2016). AJT and GRS acknowledge support from NSERC Discovery Grants (RGPIN-402752-2011 & RGPIN-06569-2016). MJM appreciates support via an STFC Ernest Rutherford Fellowship. DA acknowledges support from the Royal Society. GEA is the recipient of an Australian Research Council Discovery Early Career Researcher Award (project number DE180100346) funded by the Australian Government. TMB acknowledges financial contribution from the agreement ASI-INAF n.2017-14-H.0. PGJ acknowledges funding from the European Research Council under ERC Consolidator Grant agreement no 647208. SM and TDR acknowledge support from a Netherlands Organisation for Scientific Research (NWO) Veni Fellowship and Vici Grant, respectively. KPM acknowledges support from the Oxford Centre for Astrophysical Surveys, which is funded through the Hintze Family Charitable Foundation. KPM is currently a Jansky Fellow of the National Radio Astronomy Observatory. This work profited from discussions carried out during a meeting on multi-wavelength rapid variability organised at the International Space Science Institute (ISSI) Beijing by T. Belloni and D. Bhattacharya. The authors acknowledge the worldwide effort in observing this outburst, and the planning tools (created by Tom Marsh and coordinated by Christian Knigge) that enabled these observations.

JCAM-J wrote the manuscript with input from all authors. JCAM-J wrote the observing proposal BM421 with help from all authors. GRS wrote the observing proposal BS249 with help from JCAM-J, AJT, RPF, PGJ, GEA and KPM. JCAM-J designed and processed the VLBA observations. AJT performed the Monte Carlo modelling. JCAM-J, AJT and GRS analysed the data. MJM led the development of the Lense-Thirring precession scenario, with help from S. Markoff.

Reprints and permissions information is available at www.nature.com/reprints.

The authors declare that they have no competing financial interests.

Correspondence and requests for materials should be addressed to J.C.A.M.-J. (email: james.miller-jones@curtin.edu.au).

{methods}

V404 Cygni was observed over fifteen epochs with the VLBA, between 2015 June 17th and July 11th (Extended Data Table 1).

0.1 Observations and data reduction.

External gain calibration was performed using standard procedures within the Astronomical Image Processing System\citemethodsgreisen03 (AIPS). We used geodetic blocks to remove excess tropospheric delay and clock errors for all observations of duration $\geq 3$ hours. Our phase reference calibrator was the bright (1.8 Jy at 15 GHz), nearby (16.6 arcmin from V404 Cygni) extragalactic source J2025+3343\citemethodsma98.

The strong amplitude variability seen in both the VLBA data and the simultaneous VLA data from 2015 June 22nd\citemaintetarenko17 violates a fundamental assumption of aperture synthesis. We therefore broke the data down into short segments, within which the overall amplitude would not change by more than 10%. This equated to 103 scan-based (70-s) segments in the 15-GHz data from June 22nd, and two-scan (310-s) segments in the 8.4-GHz data from the other epochs. The sparse uv-coverage in each individual segment meant that we could not reliably image complex structures. We therefore minimised the number of degrees of freedom during deconvolution and self-calibration by performing uv-model fitting using the Difmap\citemethodsshepherd97 software package (v2.41), rather than the standard clean algorithm. With this approach, we found that the source could always be represented by a small number ($\leq 6$) of point source components. To create the final images, we performed multiple rounds of phase-only self-calibration, and a final single round of amplitude and phase self-calibration (leaving noise-like residuals in all cases).

Since this version of Difmap did not provide uncertainties on the fitted model parameters, we used the Common Astronomy Software Application\citemethodsmcmullin07 (CASA; v4.7.2) to fit the self-calibrated data with the software tool UVMULTIFIT\citemethodsmarti-vidal14. We used the Difmap model fit results to define both the number of point sources used for each snapshot and the initial guesses for their positions and flux densities.

Given the sparse uv-sampling, we took additional steps to ensure the fidelity of our final images, taking guidance from previous time-resolved VLBI studies\citemethodsfomalont01. We examined each snapshot image to check for consistency between adjacent frames. Only a small minority of frames showed inconsistent structure, and were therefore reprocessed using prior knowledge from the adjacent frames. In a few cases, we imaged longer chunks of data (10–15 min) to assess the fidelity of the structures with better uv-coverage. As seen in Extended Data Figures 3–4, the positions and flux densities of our final set of components evolve smoothly with time (other than occasional jumps when a new component appears or a blend of two components separates sufficiently to become resolved). This gives us confidence in the fidelity of our images.

0.2 Markov Chain Monte Carlo analysis.

Short-timescale tropospheric phase variations, particularly at 15.4 GHz, coupled with the propensity of self-calibration to shift source positions by a small fraction of a synthesised beam combine to introduce low-level positional offsets between individual snapshots. While these would be averaged out in longer data segments, they affected the fitted component positions in our snapshot images. Furthermore, in snapshots made with fewer than 10 antennas (e.g. due to the source having set), poor uv-coverage made it hard to distinguish the true source position from the high sidelobes, and the initial peak position selected to start the model-fitting process dictated the astrometric registration of the final image.

To fit for the proper motions of the individual point source components on June 22nd, we first had to determine the positional offsets in each snapshot. We assumed ballistic motion and constructed a set of linear equations with $k$ ejecta components and $i$ images, such that,

[TABLE]

[TABLE]

where $\mu_{{\rm ra},k}$ and $\mu_{{\rm dec},k}$ represent the proper motions of the $k$th component, and $t_{{\rm ej},k}$ its ejection time. The atmospheric jitter parameters ${\rm J}_{{\rm ra},i}$ and ${\rm J}_{{\rm dec},i}$ represent the offsets in position for the $i$th image, allowing us to correct the positional shifts.

With $k=10$ moving components (labelled by ejection time and direction of motion; see Extended Data Table 3), and $i=103$ images, we had 359 individual measurements in both right ascension and declination. This translates to 20 linear equations, and 236 free parameters. We took a Bayesian approach for parameter estimation, simultaneously solving equations (1) and (2) using a Markov-Chain Monte Carlo (MCMC) algorithm implemented with the emcee package\citemethodsforeman-mackey13. Prior distributions for all parameters are listed in Extended Data Table 4. Lastly, due to the large number of rapidly-moving ejecta and the blending of components close to the core, it was occasionally difficult to distinguish between components. We therefore assigned a confidence flag to each component for each image prior to the fitting (H = high, M = medium, L = low, and B = possible blended component) and weighted the data according to these flags (H=1, M=0.7, L=0.3, and B=0.1).

The best fitting results (Extended Data Table 3) were taken as the median of the posterior distributions from the converged MCMC solution, with the $1\sigma$ uncertainties reported as the range between the median and the 15th/85th percentile. Two components, N8 and N9, did not appear to move away from the core. Given the faint nature of the components and the sparse uv-coverage, these could be artifacts arising from the difficulty of representing complex structures with a small number of unresolved point sources.

0.3 Jet dynamics and Doppler boosting.

From the similarities in ejection time and position angle, we identified three likely pairs of components (N2/S2, N3/S3, N6/S6). In all cases, the proper motion of the northern component exceeded that of its southern counterpart, implying that the northern jets are approaching and the southern jets receding. This identification is supported by the first six epochs of our 8.4-GHz VLBA data, which all showed extensions to the north (see Figure 1), consistent with the northern components being both faster-moving and more Doppler-boosted. Furthermore, only with approaching northern components do we get constraints on $\beta\cos\theta$ for the individual ejecta that are consistent with paired ejections (see Figure 4).

Assuming that our identification of pairs was correct, we then re-fit the proper motions of these three pairs, tying the ejection times of each component in a pair. We use the results of these tied fits in Figures 2–4, and Extended Data Figures 1–3, and to calculate the jet physical parameters in Extended Data Table 5.

Assuming intrinsically symmetric jets at a distance $d$, we can determine the jet speed and inclination angle from the proper motions of corresponding approaching and receding components via

[TABLE]

With a known distance, equations (4) and (5) can be uniquely solved, allowing us to derive the jet Lorentz factor, $\Gamma=\left(1-\beta^{2}\right)^{-1/2}$ and the Doppler factors $\delta_{\rm app,rec}=\Gamma^{-1}\left(1\mp\beta\cos\theta\right)^{-1}$ (see Extended Data Table 5). For unpaired ejecta, we can only solve equation (3) for $\beta\cos\theta$.

Given our estimated precession cone half-opening angle of $\approx 18$°, the N2/S2 and N3/S3 pairs have inclinations consistent with being on the surface of a precession cone centred on the binary orbital angular momentum vector. However, the N6/S6 pair has a very low inferred inclination of $14.0\pm 0.8$°. Either these two ejecta do not form a corresponding pair, or (more likely) the proper motion of N6 is affected by additional, unaccounted systematic uncertainties due to its slow motion and the short lever arm in time (it is based on only six points). Thus this last pair should be treated as less reliable than the other two. Even should N6 have been ejected slightly later, its observed angular separation suggests an ejection time prior to 13:40 UT, so our robust upper limit on the precession timescale remains a few hours.

0.4 Mass accretion rate.

The slim-disk geometry inferred from the X-ray emission implies an accretion rate at or above Eddington. Further, the walls of the slim disk are likely to obscure the hottest inner regions of the accretion flow, implying an intrinsic luminosity higher than the maximum observed value of twice the Eddington luminosity ($2L_{\rm Edd}$)\citemainmotta17. Furthermore, a supercritical accretion disk is expected to launch a powerful outflow, which can expel a significant fraction of the infalling mass\citemainpoutanen07. Recent X-ray studies of ultraluminous X-ray sources have suggested that the wind kinetic power could be a few tens of times the bolometric luminosity (albeit reduced by the covering factor and solid angle of the wind)\citemethodspinto16,pinto17. The mass accreted during the 2015 outburst was inferred to be a factor of three lower than the mass transferred from the secondary over the preceding 26-year quiescent period\citemethodsziolkowski18. This was attributed to substantial wind mass loss, either from the outer disk\citemainmunoz-darias16 or from the inner regions\citemainmotta17. A total outer mass accretion rate of order ten times the Eddington rate would therefore be plausible, and would be sufficient to give rise to a precession period of order a minute (Extended Data Figure 5a).

The average bolometric luminosity over the outburst has been estimated as $\approx 0.1L_{\rm Edd}$\citemethodsziolkowski18, suggesting that the outer mass accretion rate likely varied substantially. This would alter both the spherisation radius $r_{\rm sph}$ and the precession period, and is consistent with the sporadic nature of the marginally-detected X-ray QPOs\citemainhuppenkothen17. This could suggest sporadic episodes of precession set by the changing mass accretion rate through the disk, rather than a long-term, stable, phase-coherent precession. Assuming that the optical polarization (attributed to jet synchrotron emission) reflects the orientation of the jet axis, the slower inferred variation of the optical polarization position angle on June 24th ($4^{\circ}$ in $\sim 30$ min)\citemethodsshahbaz16 would support this scenario.

0.5 Precession mechanisms.

Various mechanisms have been put forward to explain X-ray binary jet precession. In the slaved disk model (as applied to SS 433), tidal forces on the equatorial bulge of a misaligned early-type donor star cause the star to precess, thereby inducing the disk and jets to precess likewise\citemethodsroberts74. However, the predicted precession period\citemethodshut81 for V404 Cygni is $\sim 100$ times the 6.5-day orbital period, and cannot explain the observed changes in the jet axis. Alternatively, massive outflows from a radiatively-warped, precessing outer disk could collimate and redirect the jets\citemainbegelman06. Existing treatments of radiatively-driven warping\citemethodswijers99,ogilvie01 again predict precession periods significantly longer than the orbital period, although they were restricted to standard thin accretion disks ($H/R<\alpha$). For more vertically-extended, super-critical disks, the outer disk (where the radiation warping instability acts most strongly) is shielded from the most luminous inner regions by the puffed up slim disk and the associated clumpy wind outflow, and radiation can be advected with the outflow, making radiative warps unlikely\citemainmiddleton18.

Resonances between the donor star orbit and the orbits of disk particles can also cause disk precession, giving rise to superhumps for systems with mass ratios $q\lesssim 0.3$\citemethodswhitehurst91. However, the predicted periods are a few per cent longer than the orbital period, and again insufficient to explain the rapid changes we observed. The tidal torque from the secondary is of order $10^{-9}$ times the Lense-Thirring torque at the spherisation radius, so cannot produce the required precession. Finally, since V404 Cygni is a dynamically-confirmed black hole, we can rule out precession driven by magnetic interactions between the compact object and the accretion disk\citemethodsmushtukov17.

0.6 Predicted precession period.

The expected Lense-Thirring precession period for an inner super-critical accretion disk rotating as a solid body is\citemainfragile07,middleton18

[TABLE]

where $M$ is the black hole mass, $a_{\ast}$ is the dimensionless black hole spin $Jc/GM^{2}$ (with $J$ being the spin angular momentum), $G$ is the gravitational constant, and $r_{\rm in}$ and $r_{\rm sph}$ are the inner and outer radii of the slim disk (the latter being the spherisation radius), with all radii given in units of the gravitational radius $r_{\rm g}=GM/c^{2}$. We assume that $r_{\rm in}$ is located at the innermost stable circular orbit. Since the structure of the outer part of a supercritical disk is set by the angular momentum carried away by the disk wind, $r_{\rm sph}$ depends on the fraction of the radiation energy $\epsilon_{\rm w}$ used to launch the wind, as\citemainpoutanen07

[TABLE]

where $\dot{m}$ is the mass accretion rate in units of the Eddington rate. The spin parameter of V404 Cygni was estimated\citemethodswalton17 as $a_{\ast}>0.92$, but without accounting for the slim disk geometry (which would require less light bending and hence a lower spin) and assumed the disk inclination to be that of the binary orbit, which our measurements show is not the case. The true spin could therefore be somewhat lower. With a black hole mass of $12^{+3}_{-2}M_{\odot}$\citemainshahbaz94, we can then estimate the precession timescale of the slim disk for a given wind efficiency $\epsilon_{\rm w}=(1+L_{\rm rad}/L_{\rm wind})^{-1}$, where $L_{\rm rad}$ and $L_{\rm wind}$ are the radiative luminosity and wind power, respectively.

Based on the peak intrinsic luminosity\citemainmotta17, and with a wind power fraction $\epsilon_{\rm w}$ of 0.25–0.5 (as estimated from relativistic magnetohydrodynamic simulations\citemethodsjiang14), slim disk models imply $15<\dot{m}<150$\citemainpoutanen07. For moderate spins, we therefore predict precession timescales of order minutes and spherisation radii of tens to hundreds of $r_{\rm g}$ (see Extended Data Figure 5). The predicted spherisation radii are consistent with the maximum radius expected for rigid precession\citemainmotta18. While the 18 mHz QPO detected simultaneously with our observations (at 11:17 UT on June 22nd) was relatively low-significance at 3.5$\sigma$, it would imply a precession timescale of 56 s. Given the uncertainty in mass accretion rate and black hole spin, this timescale is roughly consistent with these predictions. Since the maximum radius for rigid precession implied by the disk alignment criterion sets a spin and aspect-ratio dependent lower limit on the precession frequency\citemainmotta18, then for an aspect ratio of $H/R=0.5$, this timescale would imply a spin of $a\lesssim 0.3$.

0.7 Jet energetics.

The minimum amount of energy required to produce a given synchrotron luminosity is\citemethodsfender06

[TABLE]

where $\eta=(1+\beta)$ and $\beta$ is the ratio of energy in protons to that in the radiating electrons, $L_{\nu}$ is the monochromatic radio luminosity (given by $L_{\nu}=4\pi d^{2}S_{\nu}$, where $S_{\nu}$ is the measured flux density), $\nu$ is the observing frequency and $V$ is the emitting volume. We make the standard assumption that there is no energy in protons ($\eta=1$). The brightest of our ejecta is knot S3, which at 12:07 UT has a flux density of 461 mJy at 15.26 GHz (Extended Data Figure 4), and is unresolved to the synthesised beam of $1.2\times 0.4$ mas2. Assuming a maximum knot radius of 0.4 mas at 2.39 kpc, we derive an upper limit on its minimum energy of $8\times 10^{38}$ erg.

While this knot would have been expanding adiabatically (with an expansion speed 0.01-0.15$c$\citemaintetarenko17), it never became significantly resolved to the VLBI beam, so should have been substantially smaller than 0.4 mas at 12:07 UT. Hence the minimum energy is likely to be significantly lower than derived above. On the other hand, if the magnetic field deviated significantly from equipartition, the energy could be somewhat higher than the minimum.

Should the precession period indeed be of order minutes, the knots would need to be launched over a timescale small enough that they were not significantly extended due to the precessional motion over the launching period. This would argue for ejection on timescales no longer than a few seconds. A lower limit on the timescale comes from the light crossing time of the jet acceleration zone, which was found to be 0.1 light seconds ($3\times 10^{9}$ cm)\citemethodsgandhi17. Alternatively, modelling the multi-frequency radio light curves gave fitted component radii of 0.6–1.3$\times 10^{12}$ cm at the peak of the sub-mm emission in each flare\citemaintetarenko17, corresponding to light crossing times of 20–40 s. Since the sub-mm emission does not come from the jet base itself, the timescale of ejection would likely be significantly shorter. In either case, our minimum energy synchrotron calculations above would not require the jets to exceed the Eddington luminosity. However, even this would not be a hard limit given recent jet power constraints from ultraluminous X-ray sources\citemethodspakull10,soria14.

0.8 Data availability

The raw VLBA data are publicly available from the National Radio Astronomy Observatory archive (https://archive.nrao.edu/archive/advquery.jsp). All software packages used in our analysis (AIPS, Difmap, CASA, UVMULTIFIT, emcee) are publicly available. The final calibrated images and uv data are available from the corresponding author upon reasonable request. The data underlying the figures are available as csv or xlsx files, and the measured positions and flux densities of all VLBA components from 2015 June 22nd are included with the MCMC fitting code (see below).

0.9 Code availability

The MCMC fitting code is available at https://github.com/tetarenk/jet-jitter.

\bibliographystylemethods

naturemag \bibliographymethodsjmj_nature

Extended Data

Supplementary Video: Movie showing the evolution of the jet morphology over four hours on 2015 June 22nd. Time (indicated in UT) has been sped up by a factor of 1000. In the 103 separate snapshot images, we identify twelve separate components, together with a persistent core. Ejected components appear to move ballistically outwards over time, with varying proper motions and position angles, implying precession of the jet axis. Images have been corrected for atmospheric jitter (see Methods). Contours are at $\pm\sqrt{2}^{n}$ times the rms noise level of 3 mJy beam*-1*, where $n=3,4,5,...$. Top color bar is in units of mJy beam*-1*.