# A shared accretion instability for black holes and neutron stars

F.M Vincentelli, J Neilsen, A.J Tetarenko, Y Cavecchi, N Castro Segura, S del Palacio, J van den Eijnden, G Vasilopoulos, D Altamirano, M. Armas Padilla, et al.

► To cite this version:

F.M Vincentelli, J Neilsen, A.J Tetarenko, Y Cavecchi, N Castro Segura, et al. A shared accretion instability for black holes and neutron stars. Nature, 2023, 615 (7950), pp.45-49. 10.1038/s41586-022-05648-3. hal-04038940

HAL Id: hal-04038940

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

Submitted on 4 Oct 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
# A shared accretion instability for black holes and neutron stars

https://doi.org/10.1038/s41586-022-05648-3

Received: 17 May 2022

Accepted: 12 December 2022

Published online: 1 March 2023

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

Check for updates

F. M. Vincentelli<fnref n="1" />,²,³,⁴✉, J. Neilsen³, A. J. Tetarenko⁵, Y. Cavecci⁶,⁷, N. Castro Segura⁴, S. del Palacio⁸, J. van den Eijnden⁹, G. Vasilopoulos¹⁰,¹¹, D. Altamirano⁴, M. Armas Padilla<fnref n="1" />,², C. D. Bailyn¹⁰, T. Belloni¹², D. J. K. Buisson⁴, V. A. Cúneo<fnref n="1" />,², N. Degenaar¹³, C. Knigge⁴, K. S. Long¹⁴,¹⁵, F. Jiménez-Ibarra<fnref n="1" />,², J. Milburn¹⁶, T. Muñoz Darias<fnref n="1" />,², M. Özbey Arabacı⁴, R. Remillard¹⁷ & T. Russell¹⁸

Accretion disks around compact objects are expected to enter an unstable phase at high luminosity<fnref n="1" />. One instability may occur when the radiation pressure generated by accretion modifies the disk viscosity, resulting in the cyclic depletion and refilling of the inner disk on short timescales². Such a scenario, however, has only been quantitatively verified for a single stellar-mass black hole³⁻⁵. Although there are hints of these cycles in a few isolated cases⁶⁻¹⁰, their apparent absence in the variable emission of most bright accreting neutron stars and black holes has been a continuing puzzle¹¹. Here we report the presence of the same multiwavelength instability around an accreting neutron star. Moreover, we show that the variability across the electromagnetic spectrum—from radio to X-ray—of both black holes and neutron stars at high accretion rates can be explained consistently if the accretion disks are unstable, producing relativistic ejections during transitions that deplete or refill the inner disk. Such a new association allows us to identify the main physical components responsible for the fast multiwavelength variability of highly accreting compact objects.

Swift J1858.6–0814 (hereafter Swift J1858) is a low-mass X-ray binary (LMXB) that was first detected in November 2018 (ref. ¹²) and reached a maximum X-ray luminosity of roughly $10^{37}$ erg s⁻<fnref n="1" /> (0.6–79 keV)¹³. Spectral analysis showed peculiar properties, including substantial obscuration¹³,¹⁴ ($N_H \approx 10^{23}$ cm⁻²) and outflows in X-rays¹⁵, optical¹⁶ and ultraviolet (UV)¹⁷. Moreover, for more than a year after its discovery, the source showed remarkable flaring activity from radio to hard X-rays¹³,¹⁵,¹⁸,¹⁹. The source returned to quiescence in 2020, but not before exhibiting X-ray eclipses¹⁹ and type I X-ray bursts²⁰, indicating the presence of an accreting neutron star with an orbital inclination >70° at a distance of about 13 kpc.

On 6 August 2019, we coordinated a multiwavelength campaign to observe Swift J1858 simultaneously for about 4 h with high time resolution in five bands: X-rays (3–79 keV) with the Nuclear Spectroscopic Telescope Array (NuSTAR); UV (150 nm) with the Cosmic Origins Spectrograph aboard the Hubble Space Telescope (HST); optical (i+z sdss band, effective wavelength $\lambda_{eff} = 720$ nm) with the RISE at the Liverpool Telescope (LT); near-infrared (IR) ($K_s$ band, $\lambda_{eff} = 2.2$ µm) with HAWK-I on the Very Large Telescope; and radio (4.5 and 7.5 GHz) with the Karl G. Jansky Very Large Array (VLA). The source showed very strong variability with similar patterns in UV, optical, IR (UV/O/IR) and X-ray (see Fig. 1a,b). On long timescales, Swift J1858 exhibited a repetitive

behaviour, alternating between quiet and active/variable phases (Figs. 1 and 2). The active phases showed oscillatory behaviour on timescales of about 100 s; we refer to these as ‘beats’, given their visual similarity to the ‘heartbeat’ variability pattern in GRS 1915+105 (ref. ⁵). On timescales of seconds, the source showed episodic fast flaring events (seen only in IR), which we refer to as ‘flares’.

To explore the multiwavelength temporal behaviour, we computed the cross-correlation function (CCF) between NuSTAR and HAWK-I for all the simultaneous segments in our dataset (see Methods). We measured a clear correlation between the two bands, but the IR lags the X-ray variability with a delay that changes from about 2.5 s to about 5.5 s (see Fig. 1c). The magnitude and orbital-phase dependence of these lags are fully consistent with a model in which the UV/O/IR beats originate from the irradiation of X-ray beats on a disk and donor star with high orbital inclination (approximately 80°) and the orbital period of Swift J1858 (about 21.3 h (ref. ¹⁹)).

Simple mass-accretion-rate variations in a hot inflow are unlikely to explain the driving X-ray light curve². The X-ray variability observed in Swift J1858 shows marked spectral evolution not compatible with standard variability of accreting compact objects³,⁴,²¹. Also, similar variability has been seen in the archetypal high-accretion-rate stellar-mass black holes GRS 1915+105 and V404 Cygni¹³. These sources also share

<footnote n="1">¹Instituto de Astrofísica de Canarias, La Laguna, Tenerife, Spain. ²Departamento de Astrofísica, Universidad de La Laguna, La Laguna, Tenerife, Spain. ³Department of Physics, Villanova University, Villanova, PA, USA. ⁴Department of Physics and Astronomy, University of Southampton, Southampton, UK. ⁵Department of Physics and Astronomy, Texas Tech University, Lubbock, TX, USA. ⁶Instituto de Astronomía, Universidad Nacional Autónoma de México, Ciudad de México, Mexico. ⁷Departament de Física, EEBE, Universitat Politècnica de Catalunya, Barcelona, Spain. ⁸Department of Space, Earth and Environment, Chalmers University of Technology, Gothenburg, Sweden. ⁹Astrophysics, Department of Physics, University of Oxford, Oxford, UK. ¹⁰Department of Astronomy, Yale University, New Haven, CT, USA. ¹¹Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, Strasbourg, France. ¹²INAF - Osservatorio Astronomico di Brera, Merate, Italy. ¹³Anton Pannekoek Institute for Astronomy, University of Amsterdam, Amsterdam, The Netherlands. ¹⁴Space Telescope Science Institute, Baltimore, MD, USA. ¹⁵Eureka Scientific, Inc., Oakland, CA, USA. ¹⁶Cahill Center for Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA, USA. ¹⁷MIT Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA, USA. ¹⁸INAF, Istituto di Astrofisica Spaziale e Fisica Cosmica, Palermo, Italy. ²e-mail: fvincentelli@iac.es</footnote>
<figure><img src="image_3.png" /></figure>

**Fig. 1** | Multiwavelength variability of Swift J1858. a, Campaign of Swift J1858 with time-resolved observations in (starting from the lowest panel) radio (VLA), IR (HAWK-I), optical (LT), UV (HST) and X-rays (NuSTAR). The zero time is 06/08/2019 00:00:00 UTC. All light curves are plotted with their standard error. For the top four light curves, an interpolation line is also plotted for clarity. The inset in the IR panel shows an example of a fast flare from the grey-shaded area. b, Strictly simultaneous X-ray UV/O/IR segment using rebinned light curves: all bands follow the same repetitive pattern. The inset shows the X-ray/IR CCF with a peak at 2.8 ± 0.5 s (see Methods). c, X-ray versus IR

lag as a function of the orbital phase of Swift J1858. Errors are computed as the standard deviation of the lag distribution (see Methods). The solid line represents the expected best fitted delays from an irradiated disk and companion star in a neutron star LMXB with an orbital period of 21.3 h and an inclination of about 80°. The dot-dashed and dashed lines correspond to the 68% and 99% confidence intervals, respectively. d, The average light-curve profile of the fast IR flares over 400 s. Error bars show the standard error. The grey curves are the overlay of all of the individual light curves used for the stack.

other important properties with Swift J1858, such as high luminosity (40% of the Eddington luminosity for Swift J1858), obscuration and outflows<sup>13,14</sup>. This association is strengthened by the remarkable similarity of the IR light curve of Swift J1858 and the X-ray light curve of the so-called 'b' variability class of GRS 1915+105 (ref. 21) (Fig. 2). Even though the patterns are less discernible in the X-ray band for Swift J1858 (probably owing to variable line-of-sight obscuration, given its high inclination<sup>9,13,15,16</sup>), the irradiation origin of the UV/O/IR light curve strongly suggests a common physical mechanism for the driving variability in both sources.

From a physical point of view, it is commonly accepted that the recurrent behaviour of GRS 1915+105 (that is, heartbeats and other limit cycles) is due to a radiation-pressure instability in the disk at high accretion rates<sup>2–5</sup>. Although not fully confirmed by general relativistic magnetohydrodynamic simulations, this instability is believed to drive cyclic accretion or ejection and rebuilding of the inner disk, generating repeating patterns in X-rays on 10–1,000-s timescales<sup>3–5</sup>. If this

emission irradiates the disk and companion star, it will give rise to a delayed UV/O/IR light curve, such as the one observed in Swift J1858. The interpretation of beats as a disk instability can be tested: both models<sup>4</sup> and observations<sup>5</sup> of GRS 1915+105 need short-lived jet ejections near the peak luminosity (roughly coincident with the depletion of the disk).

The fast IR flares in Swift J1858 seem to verify this hypothesis, giving credence to the radiation-pressure instability interpretation of the limit cycles. Aligning and averaging the flares, including 200 s of data before and after each flare, shows that they take place after the peak of the slower IR beats (see Fig. 1d). But these flares are inconsistent with a thermal origin (see Methods) and, given their red colour, we interpret them as direct evidence of optically thin synchrotron emission from transient/short-lived relativistic jet ejections expected to occur<sup>4</sup> during these beat oscillations.

Swift J1858 also showed substantial radio variability throughout our campaign<sup>18</sup>, which requires explanation. The fast IR flares cannot be responsible for the observed low-frequency variability because

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

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

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 2 | The unified accretion instability for black holes and neutron stars.**
*a, HAWK-I light curve of the accreting neutron star (NS) Swift J1858 (red line) compared with the Chandra X-ray light curve of the accreting black hole (BH) GRS 1915+105 (blue line). The light curves were rebinned to 10 s and aligned to their minimum value (at approximately 11,000 s). The similarity of the limit-cycle timescales is evident and strongly indicates a common physical mechanism. **b, Drawing of the underlying physical mechanism that explains the*</figcaption></figure>

multiwavelength phenomenology of both sources. A limit cycle in the accretion rate *m* generates 'beats' (or 'heartbeats') when the disk reaches the luminous phase, along with radio ejecta ('baby jets'). Fast IR flares, instead, originate from optically thin synchrotron emission, which does not contribute to the overall radio emission. Depending on the contribution from obscuration, reprocessing and the jet, the beat oscillations can be detected in X-rays or at UV/O/IR wavelengths.

their amplitude and duration would naturally lead to their radio emis-
sion being completely self-absorbed (τ ≫ 1 at 10 GHz; see Methods).
However, observations of GRS 1915+105 also show 'baby jets': strong
radio flares (although their synchrotron emission can contribute sub-
stantially in the IR band<sup>22,23)</sup> that are consistent with emission from
adiabatically expanding blobs<sup>24</sup> (although their launching mechanism
is still not clear). To search for baby jets in Swift J1858 and make a com-
parison with GRS 1915+105, we modelled its variable radio emission as
the sum of several ejecta<sup>25</sup>, performing the same analysis on an archi-
val radio observation of GRS 1915+105 (coincident with the β-class
X-ray light curve shown in Fig. 2). The results presented in Fig. 3 show
that the radio variability of both sources is closely reproduced by our
modelling. For Swift J1858, the model suggests baby jet ejection times
(grey-shaded areas in Fig. 3) near quiet/active phase transitions; most
of the ejecta in GRS 1915+105 occur during quiet phases but several fall
close to quiet/active transitions as well.

For self-consistency, we then tested whether the baby jets of Swift
J1858 would be detectable in the IR as for GRS 1915+105. Previous
studies<sup>5,24</sup> show accretion instabilities in GRS 1915+105 when the X-ray
and radio luminosity are L<sub>BH</sub> ≈ 10<sup>38</sup> erg s<sup>-1</sup> and L<sub>BH,radio</sub> ≈ 10<sup>30</sup> erg s<sup>-1</sup>, res-
pectively. For Swift J1858, we find L<sub>NS</sub> ≈ 10<sup>37</sup> erg s<sup>-1</sup> and L<sub>NS,radio</sub> ≈ 10<sup>29</sup> erg s<sup>-1</sup>
(ref. <sup>18</sup>). Even under the conservative assumption that the ratio between
the IR and radio flux from the jet in Swift J1858 is the same as that
observed in GRS 1915+105 during the β-class instability (IR/radio ≈ 1.4)<sup>24</sup>,
then we expect an IR baby jet flux of only approximately 0.24 mJy. This
is almost a factor of two fainter than the reprocessed emission during
the beats (roughly 0.4 mJy). This indicates that the two sources share

the same disk-jet coupling, despite having qualitatively different radio
and IR light curves. More broadly, regardless of the jet-launching
mechanism, this shows how the presence of accretion instabilities can
depend not only on the accretion rate and disk-jet geometry but also
on the binary orbit and the mass of the compact object.

There is growing evidence that high-accretion-rate black hole
sources such as GRS 1915+105, V4641 Sagittarii, Cygnus X-3 and V404
Cygni all share common X-ray spectral variability properties¹⁴. How-
ever, multiwavelength parallels have proved more difficult owing to
their different extinctions, hampering efforts to produce a unified
physical scenario for this class of sources. Yet, as anticipated from our
conclusions, Swift J1858 shows clear analogies with all these objects.
Simultaneous multiwavelength observations of the 2015 outburst
of V404 Cygni showed repetitive optical/X-ray patterns with a lag
consistent with reprocessing¹⁰,²⁶,²⁷ and fast non-thermal flares²⁸. Fur-
thermore, its extreme radio variability is consistent with jet ejections
taking place during X-ray spectral transitions²⁵. Moreover, similar
O/IR recurrent patterns with comparable timescales have also been
observed in V4641 Sagittarii²⁹ and Cygnus X-3 (ref. ³⁰). Finally, we
note that X-ray heartbeats have also been detected in sources such
as the LMXB IGR J1709−3624 (ref. ⁷) and the ultraluminous X-ray
source (ULX) NGC 3621 (ref. ³¹), which also shows notable line-of-sight
obscuration despite having a lower inclination. Thus, the recent
association of Swift J1858 as a low-luminosity Z source³², and the
isolated presence of X-ray ‘GRS 1915-like’ patterns in other accreting
neutron stars such as the Rapid Burster⁶ and the Bursting Pulsar³³,
strongly indicate that Swift J1858 represents the missing link for
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 3 | Radio oscillations modelling.** a,b, Radio time series in the two frequency bands (green and yellow points), overplotted with the total jet model (solid black line) and the individual model contributions from the jet and the counterjet for each ejection (dashed blue and red lines) of Swift J1858 (a)</figcaption></figure>

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

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

and GRS 1915 (b). Residuals (data-model/uncertainties) are also shown for each fit. c,d, Modelled jet ejection times (dashed black lines; grey shading indicates 1σ confidence intervals) with respect to the IR for Swift J1858 (c) and X-ray for GRS 1915 (d).

multiwavelength variability in high-accretion-rate sources (Fig. 2 and Extended Data Fig. 1).

It was also noted during review that, although the limit-cycle timescale is similar in GRS 1915+105 and Swift J1858 (despite their very different masses; see Methods), the beat timescale is much shorter around the black hole in the example light curves shown in Fig. 2. In fact, GRS 1915+105 exhibits a wide range of beat durations in similar limit cycles²¹, which suggests that the beats may represent a second instability timescale⁴ or may be affected by other factors in the accretion flow. One possibility is the jet power, which is expected to have a marked impact on the disk structure and thus on the observed X-ray light curve³,⁴. A careful comparison of the time-dependent radio/O/IR properties in states or sources with different beat timescales³⁴ could further explain the role of jets in shaping these instabilities.

Our results draw a new coherent picture that links key aspects of the multiwavelength variability of both black holes and neutron stars at high accretion rate: recurrent repetitive patterns, radio oscillations and fast flaring. At these high accretion rates, the accretion disk becomes unstable, resulting in disk-jet cycles on timescales of about 10 s to about 1,000 s. These have historically been observed in X-rays, but our work shows that, given the right conditions (such as inclination, orbital period, obscuration and the relative brightness of the jet), accretion instabilities may in fact be more readily observable at UV/O/IR wavelengths. These instabilities are also observationally associated with radio-emitting discrete ejections; therefore—for the first time—we can define a consistent physical scenario that can quantitatively account for most of the multiwavelength variability observed from accreting compact objects at high luminosity. We argue that accretion instabilities, irradiation/obscuration and jet ejecta should be seen as three fundamental pillars that can be used to study other classes of objects accreting near the Eddington limit. With this insight, future time-resolved multiwavelength campaigns

on compact objects will lead to better constraints on the physics of these instabilities and their hosts, independently of the nature of the central object⁸.

## Online content

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-022-05648-3.

1. Lightman, A. P. & Eardley, D. M. Black holes in binary systems: instability of disk accretion. Astrophys. J. 187, L1 (1974).

2. Belloni, T., Méndez, M., King, A. R., van der Klis, M. & van Paradis, J. An unstable central disk in the superluminal black hole X-ray binary GRS 1915+105. Astrophys. J. 479, L145–L148 (1997).

3. Janiuk, A., Czerny, B. & Siemiginowska, A. Radiation pressure instability as a variability mechanism in the microquasar GRS 1915+105. Astrophys. J. 542, L33–L36 (2000).

4. Nayakshin, S., Rappaport, S. & Melia, F. Time-dependent disk models for the microquasar GRS 1915+105. Astrophys. J. 535, 798–814 (2000).

5. Neilsen, J., Remillard, R. A. & Lee, J. C. The physics of the “heartbeat” state of GRS 1915+105. Astrophys. J. 737, 69 (2011).

6. Bagnoli, T. & int’Zand, J. J. M. Discovery of GRS 1915+105 variability patterns in the Rapid Burster. Mon. Not. R. Astron. Soc. 450, L52–L56 (2015).

7. Altamirano, D. et al. The faint “heartbeats” of IGR J17091–3624: an exceptional black hole candidate. Astrophys. J. 742, L17 (2011).

8. Janiuk, A. & Czerny, B. On different types of instabilities in black hole accretion discs: implications for X-ray binaries and active galactic nuclei. Mon. Not. R. Astron. Soc. 414, 2186–2194 (2011).

9. Janiuk, A., Grzedzielski, M., Capitanio, F. & Bianchi, S. Interplay between heartbeat oscillations and wind outflow in microquasar IGR J17091-3624. Astron. Astrophys. 574, A92 (2015).

10. Kimura, M. et al. Repetitive patterns in rapid optical variations in the nearby black-hole binary V404 Cygni. Nature 529, 54–58 (2016).

11. Done, C., Wardziński, G. & Gierliński, M. GRS 1915+105: the brightest Galactic black hole. Mon. Not. R. Astron. Soc. 349, 393–403 (2004).

12. Krimm, H. A. et al. Swift reports the detection of a new galactic transient source Swift J1858.6-0814. The Astronomer's Telegram No. 12151 (2018).
13. Hare, J. et al. NuSTAR observations of the transient galactic black hole binary candidate Swift J1858.6-0814: a new sibling of V404 Cyg and V4641 Sgr? *Astrophys. J.* **890**, 57 (2020).

14. Koljonen, K. I. I. & Tomsick, J. A. The obscured X-ray binaries V404 Cyg, Cyg X-3, V4641 Sgr, and GRS 1915+105. *Astron. Astrophys.* **639**, A13 (2020).

15. Buisson, D. J. K. et al. Soft X-ray emission lines in the X-ray binary Swift J1858.6-0814 observed with XMM-Newton Reflection Grating Spectrometer: disc atmosphere or wind? *Mon. Not. R. Astron. Soc.* **498**, 68–76 (2020).

16. Muñoz-Darias, T. et al. The changing-look optical wind of the flaring X-ray transient Swift J1858.6-0814. *Astrophys. J.* **893**, L19 (2020).

17. Castro Segura, N. et al. A persistent ultraviolet outflow from an accreting neutron star binary transient. *Nature* **603**, 52–57 (2022).

18. van den Eijnden, J. et al. The variable radio counterpart of Swift J1858.6-0814. *Mon. Not. R. Astron. Soc.* **496**, 4127–4140 (2020).

19. Buisson, D. J. K. et al. Dips and eclipses in the X-ray binary Swift J1858.6-0814 observed with NICER. *Mon. Not. R. Astron. Soc.* **503**, 5600–5610 (2021).

20. Buisson, D. J. K. et al. Discovery of thermonuclear (Type I) X-ray bursts in the X-ray binary Swift J1858.6-0814 observed with NICER and NuSTAR. *Mon. Not. R. Astron. Soc.* **499**, 793–803 (2020).

21. Belloni, T., Klein-Wolt, M., Méndez, M., van der Klis, M. & van Paradijs, J. A model-independent analysis of the variability of GRS 1915+105. *Astron. Astrophys.* **355**, 271–290 (2000).

22. Fender, R. P. & Pooley, G. G. Infrared synchrotron oscillations in GRS 1915+105. *Mon. Not. R. Astron. Soc.* **300**, 573–576 (1998).

23. Eikenberry, S. S., Matthews, K., Morgan, E. H., Remillard, R. A. & Nelson, R. W. Evidence for a disk-jet interaction in the microquasar GRS 1915+105. *Astrophys. J.* **494**, L61–L64 (1998).

24. Mirabel, I. F. et al. Accretion instabilities and jet formation in GRS 1915+105. *Astron. Astrophys.* **330**, L9–L12 (1998).

25. Tetarenko, A. J. et al. Extreme jet ejections from the black hole X-ray binary V404 Cygni. *Mon. Not. R. Astron. Soc.* **469**, 3141–3162 (2017).

26. Alfonso-Garzón, J. et al. Optical/X-ray correlations during the V404 Cygni June 2015 outburst. *Astron. Astrophys.* **620**, A110 (2018).

27. Hynes, R. I. et al. Optical and X-ray correlations during the 2015 outburst of the black hole V404 Cyg. *Mon. Not. R. Astron. Soc.* **487**, 60–78 (2019).

28. Dallilar, Y. et al. A precise measurement of the magnetic field in the corona of the black hole binary V404 Cygni. *Science* **358**, 1299–1302 (2017).

29. Uemura, M. et al. Optical observation of the 2003 outburst of a black hole X-Ray binary, V4641 Sagittarii. *Publ. Astron. Soc. Jpn.* **56**, 823–829 (2004).

30. Fender, R. P., Bell Burnell, S. J., Williams, P. M. & Webster, A. S. Flaring and quiescent infrared behaviour of Cygnus X-3. *Mon. Not. R. Astron. Soc.* **283**, 798–804 (1996).

31. Motta, S. E. et al. The slow heartbeats of an ultraluminous X-ray source in NGC 3621. *Astrophys. J.* **898**, 174 (2020).

32. Rhodes, L. et al. Long-term radio monitoring of the neutron star X-ray binary Swift J1858.6-0814. *Mon. Not. R. Astron. Soc.* **513**, 2708–2718 (2022).

33. Court, J. M. C. et al. The evolution of X-ray bursts in the ‘bursting pulsar’ GRO J1744–28. *Mon. Not. R. Astron. Soc.* **481**, 2273–2298 (2018).

34. Rothstein, D. M., Eikenberry, S. S. & Matthews, K. Observations of rapid disk-jet interaction in the microquasar GRS 1915+105. *Astrophys. J.* **626**, 991–1005 (2005).

**Publisher’s note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

© The Author(s), under exclusive licence to Springer Nature Limited 2023
**Data reduction: Swift J1858**

**NuSTAR.** The NuSTAR (3–79 keV)³⁵ observed the source between 18:46:48 UTC 5 August 2019 and 08:26:45 UTC 6 August 2019 (OBSID: 90501333002). Barycentre-corrected events were extracted using the software nipupeline. Then, using XSELECT, we filtered the events for both the FPMA and FPMB between 3 and 79 keV (that is, between channels 35 and 1,935) and within a region of 120" from the centre of the source. We used the software HENDRICS (http://ascl.net/1805.019) to create the light curve of 7.8125 ms.

**HST.** Swift J1858 was observed in the far-UV with the Cosmic Origin Spectrograph and the G140L grating using the Primary Science Aperture with a total exposure of 4.9 ks (program GO/DD 15984, principal investigator N. Castro Segura). The observations were obtained in TIME-TAG mode, recording events with a time resolution of 32 ms. Data were reduced using the HST CALCOS and CALSTIS pipelines (provided by the Space Telescope Science Institute, https://github.com/spacetelescope). We extracted light curves from the TIME-TAG event files with a time resolution of 32 ms, using the same regions defined by the pipeline. An empirical background correction was applied. Regions affected by geocoronal airglow emission associated with Lyman α (λ = 1,208–1,225 Å) and O II (λ = 1,298–1,312 Å) were masked when extracting the light curves.

**LT.** We obtained 4.5 ks (1.25 h) of high-time-resolution photometry with the RISE fast-readout camera attached to the 2-m Liverpool Telescope at the Observatorio del Roque de los Muchachos during the first section of the multiwavelength campaign. The observations were carried out using the 720-nm long-pass filter (approximately i+z band). The fast readout time of RISE (approximately 0.04 s) resulted in 10-s time resolution. We used the data-reduction pipeline for the optical-imaging component of the Infrared-Optical suite. The flux calibration against a field star, catalogued in PanSTARRS (broad-band i filter), was performed using astropy/photutils.

**HAWK-I on the Very Large Telescope.** We collected 4 h of near-IR time-resolved photometry in the Ks band (effective wavelength, λeff = 2.2 µm) with HAWK-I³⁶ mounted on the Very Large Telescope UT4 Yepun under program 103.201A.001. The instrument consists of four Hawaii 2RG 2,048 × 2,048-pixel detectors. However, by using the FAST-PHOT mode, the instrument reads only a 128 × 64-pixel window per quadrant, allowing us to achieve a time resolution of 0.2 s. The resulting images were stacked in 'data cubes' of 250 frames, separated by a readout gap of ≈20 s. The data were reduced with the ULTRACAM software tools³⁷. Aperture parameters were derived using a bright reference star (12.51 mag, approximately 22" west from the target) average point spread function and position during each cube. To mitigate the seeing effects, the target and a comparison star (13.11 ± 0.03 mag, about 33" west from the target) were normalized to the reference star. Time stamps were put in a barycentred reference frame using the software developed by Eastman et al.³⁸.

**VLA.** Coordinated radiofrequency observations of Swift J1858 were gathered with the VLA on 6 August 2019 for 3.05 h (about 11 ks) under the program 19A-495 (principal investigator J. van den Eijnden). The array was in the A configuration, observing in C band (4–8 GHz) with the 8-bit samplers. The correlator was set up to generate two base bands centred at 4.5 and 7.5 GHz, with a bandwidth of 1,024 MHz per base band. For the observations, 3C286 was used as a flux calibrator and J1832–1035 was used as a phase calibrator. Flagging and calibration of the data were carried out using the Common Astronomy Software Applications package (CASA v4.7.2 (ref. ³⁹)), with further details on the data reduction reported in the first radio study on this source¹⁸.

To obtain a high-time-resolution light curve of the target source, we split the observation into a set of 3-min-timescale chunks and imaged each chunk in both base bands. We used Briggs weighting with the robust parameter set to 0 to maximize sensitivity and reduce the side-lobe effects of a nearby background source during this imaging process. Flux densities were then measured by fitting a point source in the image plane (using the imfit task) for each time chunk.

**Data reduction: GRS 1915+105**

**Chandra.** Chandra observed GRS 1915+105 for 120 ks beginning on 21 June 2011 at 04:36:42 UTC. The data were taken with the High Energy Transmission Grating Spectrometer in graded, continuous-clocking mode, resulting in a one-dimensional image with 8.55-ms time resolution. The data were processed with standard CIAO data-reduction tools, including `acis_process_events`. We used `axbary` to correct the arrival times of photons to the solar system barycentre. Because the source was so bright, telemetry saturation caused about 5% of the data frames to be lost from the detector chip ACIS-S3. We therefore used `dmextract` to create 1-s light curves from each individual detector chip, then coadded them and interpolated over dropped frames to create the light curve in Fig. 2.

**VLA.** Coordinated radiofrequency observations of GRS 1915+105 were gathered with the VLA on 22 June 2011 for 10.8 h (about 39 ks) under the program SC0666 (principal investigator J. Neilsen). The array was in the A configuration, observing in the C band (4–8 GHz). These data were taken in an early science mode in which the correlator was set up to generate two base bands centred at 4.2 and 7.8 GHz, with a bandwidth of 128 MHz per base band. For these observations, 3C286 was used as a flux calibrator and J1856+0610 was used as a phase calibrator. Flagging and calibration of the data were carried out using CASA v5.6, using standard procedures outlined in the CASA Guides (https://casaguides.nrao.edu) for VLA data reduction (that is, a priori flagging, setting the flux density scale, initial phase calibration, solving for antenna-based delays, band-pass calibration, gain calibration, scaling the amplitude gains and final target flagging). To obtain a high-time-resolution light curve of the target source, we split the observation into a set of 2-min-timescale chunks and imaged each chunk in both base bands. We used natural weighting to maximize sensitivity during this imaging process. Flux densities were then measured by fitting a point source in the image plane (using the imfit task) for each chunk.

**Data analysis**

**CCF.** We measured the CCF between the NuSTAR and HAWK-I with a 0.2-s bin and a maximum shift of 100 s for each NuSTAR orbit. The very strong spikes in the IR can introduce spurious peaks in the correlation coefficient ²⁶,²⁷,⁴⁰,⁴¹. Therefore, we computed the CCF excluding and including the short IR flares (that is, removing data points within ±1 s of the flares). We used the discrete CCF method ⁴²,⁴³, rebinning the X-ray events to the HAWK-I time resolution. The lag was evaluated by simulating $N = 10^4$ CCFs, randomizing each point using a Gaussian distribution with a $\sigma$ corresponding to the empirical error. For each run, the centroid was then computed with a weighted average of the correlation coefficient over a range of 80% from the maximum⁴⁰. The lag was then evaluated from the resulting distribution of centroid lags. The error was obtained from the standard deviation of the obtained distributions. A clear evolution as a function of time of the lags can be seen in Extended Data Fig. 2a.

We also checked the IR/UV lags. The HST light curve has an absolute timing accuracy up to about 1 s. We therefore computed the CCF between the HST and HAWK-I to search for any possible lags with a time resolution of 1 s. As shown in Extended Data Fig. 2b, the CCF peaks at 0. Some asymmetry is seen for lags of roughly 50 s, which can be explained by the different shapes of the flares.
**Fast flaring.** To study the properties of the fast flares, we selected spikes brighter than 0.7 mJy, finding 11 events in the light curve. We then selected 400-s segments centred on the peak values of each flare and averaged them (see Fig. 1d). Given the duration of each data cube of 50 s, we filled the readout gaps using a linear interpolation between the points at the edge, with the addition of Gaussian noise. One flare (taking place around 6,300 s from the start of the observations) showed two peaks separated by about 60 s and was not included in the average. However, its inclusion does not affect the overall trend (Extended Data Fig. 3), confirming the connection between short and long timescales.

If these flares are produced by thermal gas with non-relativistic bulk motion, their typical timescale of about 1 s indicates a maximum size of the emitting region of $R \approx 10^{10}$ cm. Furthermore, their average IR flux density of approximately 1 mJy (which corresponds to about $10^{-14}$ erg cm$^{-2}$ s$^{-1}$), together with the source distance of 13 kpc (ref. $^{20}$), yields a brightness temperature of approximately $10^6$ K. Under these conditions, the thermal gas would emit an X-ray flux of about $10^{-8}$ erg cm$^{-2}$ s$^{-1}$, which is roughly one order of magnitude higher than the observed flux. This discrepancy cannot be solved assuming a smaller IR emitting region, as that would lead to a higher brightness temperature and X-ray flux. Even considering effects of obscuration, with an $N_H$ of $10^{23}$ cm$^{-2}$ and a covering fraction of 0.9 (that is, substantially higher than the values reported by previous X-ray studies$^{13}$), would give an observed X-ray flux of at least several factors higher. We can therefore rule out a thermal origin for the flares. Given also their red colour$^{44}$, we interpret them as direct evidence of optically thin synchrotron emission from transient/short-lived relativistic jet ejections expected to take place$^{4,45}$ during these beat oscillations.

**Reprocessing modelling.** We tested whether the evolving lags between X-ray and IR are compatible with reprocessing, with a previously used code simulating the expected delay between an X-ray flare from the neutron star and the corresponding peak in the IR owing to the irradiation of the disk and the companion$^{46,47}$.

For the binary, we used the parameters of Swift J1858 (ref. $^{19}$). In particular, we set the orbital period to 76,841.3 s and tested inclinations between 70° and 90°. We fixed the neutron star mass to $1.4 M_\odot$ and used the same relation between the inclination and the mass ratio defined in previous studies$^{19}$. For the companion temperature, we adopted 6,100 K (Castro Segura et al. in preparation). Because it is better to have a sharp X-ray flare to measure the delay, we used a triangular function of width 20 s, with the peak at 10 s. We modelled the binary at a range of orbital phases from about 20° to about 90° to cover the five epochs at which we have measured lags.

The delay between the X-ray and the IR peaks depends slightly on the shape we use for the X-ray flare, but we do find lags compatible with the observations for an inclination around 80°. In general, we find that, when the companion is not lying too close to the line between the observer and the neutron star, the lag is dominated by the IR reprocessing off the companion. This lag diminishes as the donor becomes more aligned with the observer-neutron star direction, until a minimum is reached. At that point, the lag begins to increase again as the outer disk reprocessing dominates. The minimum lag allowed increases with inclination $i$, because the mass ratio changes with $i$. We fit our theoretical models to two sets of measured lags: the ones obtained including the flares in the light curves and the ones excluding them (see left panels of Extended Data Fig. 4, top and bottom, respectively). The data points follow the theoretical trend remarkably well in both cases. When we include the flares, the posterior distribution is slightly skewed, preferring higher inclinations, with median at $i = 80.3^\circ$ (0.5th, 16th, 84th and 99.5th percentiles at 76.2°, 78.4°, 82.0° and 84.1°, respectively). When we exclude the flares, owing to the lower value of the lag at phase around 45° (which is also the least precise), the distribution prefers lower inclinations. In any case, the result is compatible with the

previous result, having median at $i = 79.3^\circ$ (percentiles at 75.6°, 77.6°, 81.1° and 83.5°). This result confirms that the lag is indeed owing to reprocessing of the X-rays. These values are fully consistent with the inclination obtained by modelling the eclipse profiles during the final stages of the outbursts$^{48}$.

## Fast IR flare modelling

A consistent interpretation of the fast IR flares needs to reconcile the following features: (1) the IR flare duration time in the observer frame ($t_{IR}$); (2) the observed frequency of the emission ($v_{IR}$); (3) the peak flux density of the flare ($F_{IR}$); and (4) that a bright flare is not simultaneously seen in the radio band. Here we introduce a simplified physical model that is capable of explaining these features and deriving extra information on the source from these constraints. In what follows, we denote quantities in the jet reference frame with a prime. The relevant transformations between reference frames are given by $t = t'/\delta$, $v = \delta v'$ and $F = \delta^2 F'$, in which $\delta(\Gamma, \theta)$ is the jet Doppler factor$^{25}$. For this calculation, we will assume that the jet is launched perpendicularly to the accretion disk, such that $\theta = i$, and with semi-relativistic velocities, $\beta = 0.7 \pm 0.2$ (jet Lorentz factor $\Gamma_i = 1.4 \pm 0.2$).

Both the synchrotron cooling time of a particle, $t_{syn} \propto B^2 E_e^{-1}$, and the characteristic frequency of the emitted photons, $v_{syn} \propto BE_e^2$, depend on the energy of the electron, $E_e$, and the magnetic field intensity in the emitter, $B$ (ref. $^{49}$). We assume that the duration of the flare is given by the electron synchrotron cooling time, that is, condition (1) translates into $t'_{syn}(E, B) = t'_{IR}$. Moreover, the observed photon frequency should match the synchrotron frequency, so that condition (2) reads $v_{IR} = v_{syn}$. Thus, from conditions (1) and (2), we can derive the magnetic field intensity in the jet and the energy of the emitting electrons ($B \approx 553 \pm 22$ G and $E_e \approx 207 \pm 8$ MeV, respectively). One caveat of this is that the magnetic field should not vary substantially during the duration of the flare in the jet frame (about 6 s). Also, one can assume that the observed flux is due to optically thin synchrotron emission from a relativistic electron population of the form $N(E)dE = N_0 E^p dE$, with $p \approx 2.2$ (see the 'Radio modelling' section). We further assume that the minimum and maximum energies of the electrons are 1 MeV and 1 TeV, respectively, and normalize the electron distribution by adopting an equipartition energy condition, $U_e = U_\beta$, which ties the energy densities in relativistic electrons and the magnetic field. With that, we can obtain $N_0$ from $U_e = \int_{E_{min}}^{E_{max}} EN(E)dE$ and use it to calculate the optically thin synchrotron spectral energy distribution at any given frequency$^{49}$. Considering that the emitter is homogeneous and spherical, from condition (3), we can estimate its size to be $R_0 \approx 4.7 \pm 0.2 \times 10^{10}$ cm. Thus, the IR-emitting jets are more compact and with much higher magnetic fields than the radio jets. Finally, we verify that condition (4) holds. By having the emitter size, the electron population, the magnetic field and the synchrotron luminosity, we can calculate the synchrotron self-absorption opacity$^{25}$. We obtain $\tau \approx 10^6 \gg 1$ at 10 GHz and $\tau \approx 10^{-6} \ll 1$ at IR frequencies, showing that the derived values are consistent with the proposed scenario. We note that the overall conclusions do not change if we assume a broader range of jet velocities: adopting $\beta \ll 1$ leads to variations in the results below 25%, whereas for $\beta = 0.99$ ($\Gamma_i = 10$), variations are within a factor of two, except for the flare time in the jet frame, which increases up to about 30 s.

## Radio modelling

Swift J1858. Jet synchrotron emission is known to produce strong radio variability, often in the form of flaring activity showing a distinct observational signature in which longer-wavelength signals present as a delayed, smoothed version of shorter-wavelength signals. As our radio data exhibit this behaviour, we attempted to model them and, in turn, quantify any relationship between the IR flares and the variations observed at longer wavelengths$^{50}$. We implement the Van der Laan synchrotron bubble model$^{50}$, using the formalism outlined in Tetarenko et al. $^{25}$. Briefly, this model describes the emission from ballistically
moving, adiabatically expanding, bipolar jet components, folding in projection and relativistic effects. To perform the modelling, we use a Markov chain Monte Carlo algorithm (emcee⁵¹), simultaneously fitting both radiofrequency bands. Here we allow each ejected pair of jet components to have a different ejection time, peak flux and speed, while tying the opening angle, inclination angle, magnetic field and electron energy distribution (assumed to follow a power law of the form $N(E)dE \propto E^{-p}dE$, in which $p$ is related to the optical depth at the peak of the flare by means of $e^{p} - (2p/3 + 1)\tau_{p} + 1 = 0$ across all components. The radii of the ejecta at the peak of the flare are determined by means

of $R_{p} = \left( \frac{F_{p}d^{2}}{\pi S_{v}(p, B_{p})} \frac{1}{1 - e^{-p}} \right)^{0.5}$, in which $S_{v}(p, B_{p})$ is the synchrotron source function dependent on the electron energy distribution and magnetic field at the peak of the flare. We use wide uniform priors for most model parameters, but sample from the known distance distribution (12.8 ± 0.18 kpc, with bounds of 9–18 kpc), and use previous constraints on inclination in the literature to define this parameter's prior (uniform between 70° and 90° (refs. ¹⁹,²⁰)). We ran two modelling cases: (1) ejection times are unconstrained (see Extended Data Fig. 5 and Extended Data Table 1) and (2) ejection times are constrained to only lie within the IR beats, as we may expect if both the IR and radio are synchrotron in origin (see Fig. 3 and Extended Data Table 2).

We find that both of the modelling runs do a reasonable job of reproducing the radio light curves, capturing most of the flaring structure across both bands, with a similar quality of fit between runs. Our best-fit models contain four jet ejection events, in which the ejection times coincide with IR flares occurring during the rapid variability periods. However, the connection between radio and IR is more pronounced in modelling case (2), as we are better able to match up specific jet components to individual IR flares with lower uncertainty. Figure 3 shows the best-fit model for modelling case (2), overlaid on our radio light curves, and compares best-fit ejection times to the IR light curves (best-fit parameters are given in Extended Data Table 2). Overall, we find that, as we only need four ejections to reproduce the radio emission, most of the IR short flares do not actually become observable radio jet ejections. However, we do notice that the jet ejection times seem to cluster near the boundaries of the rapid variability periods.

**GRS 1915+105.** With the GRS 1915+105 data, we follow the same modelling procedure as that of Swift J1858, allowing each ejected pair of jet components to have a different ejection time, peak flux and speed, while tying the opening angle, inclination angle, magnetic field and electron energy distribution across all components. We use wide uniform priors for most model parameters but sample from the known distance (8.6 ± 2.0 kpc (ref. ⁵²)) and inclination angle (60° ± 5° (ref. ⁵²)) distributions. With GRS 1915+105, we do not place any constraints on ejection times. Here we use an incremental approach, building up our model flare by flare, to find the minimum number of components we need to reasonably reproduce the radio light curves. Furthermore, as we find that the data do not allow us to put tight constraints on the electron energy distribution index, we fix this parameter to $p = 2.2$ (equivalent to $\tau_{p} = 1.6$, as expected from a Fermi acceleration mechanism by means of a single shock, commonly assumed for X-ray binary jets³³⁻⁵⁵).

We find that, in GRS 1915+105, we can reproduce the radio light curves reasonably well with 11 jet ejection events. Figure 3 shows the best-fit model overlaid on our radio light curves and compares best-fit ejection times to the X-ray light curves (best-fit parameters are given in Extended Data Table 3). On comparing the best-fit parameters between Swift J1858 and GRS 1915, we find that the magnetic field strength of the ejecta are different between the two sources. This difference seems to be a natural consequence of the model. Given that the ejecta radii in the model are proportional to magnetic field strength and peak flux (among some other dependencies such as distance), to produce similar ejecta radii in both sources (as expected from the similar flare duration timescales), the lower peak fluxes in Swift J1858 with respect

to GRS 1915 naturally result in a larger magnetic field estimate for Swift J1858. Moreover, similar to what has been observed in previous studies²⁴,⁵⁶, the ejection times in GRS 1915 seem to lie in the dips between the rapid variability periods (seen here in the X-ray rather than the IR of Swift J1858). However, we notice that a subset of lower peak flux components have ejection times that lie close to the boundaries of the rapid variability periods. Notably, these fainter components seem to follow brighter ejections that tend to be launched closer to the beginning of a dip period rather than in the middle of it.

To verify our key result about the timing of the jet ejections in GRS 1915+105 and to evaluate how model dependent such a result was, we performed further modelling runs under different conditions. In particular, we performed two other modelling runs using an eight-jet component model (that is, without the extra lower peak flux components) with the electron energy distribution index fixed and free. The fixed model can reproduce the shoulders in the early flares but still has trouble reproducing the later peaks, whereas the free model can reproduce the flare decays well but not the shoulder features. Furthermore, the free model favours a much lower electron energy distribution index of $p \approx 0.8$, which lies outside the typical range ($p = 2–3$) produced with the Fermi acceleration by means of a single shock mechanism commonly assumed in these types of jet. However, regardless of the differing model properties used in these extra runs, we find that all modelling runs still produce the same ejection time pattern, with the jet ejection times coinciding with the dip periods between variable sections of the X-ray light curve.

## Timescale comparison

The viscous timescale defines the time it takes for a mass-accretion-rate fluctuation to travel through the disk at a given radius. This can be expressed as⁵⁷

$$t_{\text{vis}} \propto \alpha^{-1} M^{-\frac{1}{2}} R^2 \dot{M}^{-\frac{7}{2}}, \quad (1)$$

in which $\alpha$ is the viscosity parameter of Shakura and Sunyaev, $M$ is the mass of the compact object, $R$ is the radius of the disk and $\dot{M}$ is the absolute accretion rate.

Because the duration of the limit cycle is the viscous time in the unstable region of the disk, we find on average a similar $t_{\text{vis}}$ for our black hole (BH) and neutron star (NS) sources. Taking the ratio of equation (1) for both sources and assuming similar values for $\alpha$, we obtain

$$\left(\frac{M_{\text{BH}}}{M_{\text{NS}}}\right)^{-1/2} \left(\frac{R_{\text{BH}}}{R_{\text{NS}}}\right)^{7/2} \left(\frac{\dot{M}_{\text{BH}}}{\dot{M}_{\text{NS}}}\right)^{-2} = 1 \quad (2)$$

As both sources exhibit instabilities, we assume that they are at the same accretion rate in Eddington units, therefore $\dot{M}_{\text{BH}}/M_{\text{NS}} = M_{\text{BH}}/M_{\text{NS}}$. From these assumptions, it follows that:

$$R_{\text{NS}} = R_{\text{BH}} \left( \frac{M_{\text{BH}}}{M_{\text{NS}}} \right)^{-5/7} \quad (3)$$

Defining $r = R/R_{\text{s}}$ in which $R_{\text{s}} = 2GM/c^2$, we find that $r_{\text{NS}} = r_{\text{BH}} \left(\frac{M_{\text{BH}}}{M_{\text{NS}}}\right)^{2/7}$. Given masses of 12.6 $M_{\odot}$ for GRS 1915+105 and 1.4 $M_{\odot}$ for Swift J1858, we have $R_{\text{NS}} = 0.2R_{\text{BH}}$ or $r_{\text{NS}} = 1.8r_{\text{BH}}$.

## Swift J1858 O/IR observations and GRS 1915+105 X-ray heartbeats comparison

To quantify the association between the Swift J1858 O/IR beats and GRS 1915+105 X-ray heartbeats, we compared the mean-normalized flux distribution of the Chandra, HAWK-I, LT and optical time series from two further epochs.

Optical data were obtained in 2019 in two epochs from the 200-inch Hale Telescope at Palomar Observatory. In June 2019, Swift
J1858 was observed using CHIMERA⁵⁸, a two-band rapid photometer mounted on the Hale Telescope. The instrument uses frame-transfer, electron-multiplying charge-coupled devices to achieve 15 ms dead time and is ideally used to study sub-second variability. Observations were obtained in g' and r' filters, whereas CHIMERA was operated using the conventional amplifier with 0.2–1.0-s exposures with 2 × 2 binning. Images were bias subtracted and flat-field corrected, whereas aperture photometry was carried out using the standard CHIMERA pipeline⁵⁸. Fluxes reported in this paper were estimated by means of differential photometry with respect to nearby stars.

On 4 October 2019, the field was observed by Palomar's Wafer-Scale camera for Prime (WASP⁵⁹) wide-field prime-focus camera using 10-s exposures, whereas the extra dead time between observations is of the same order. Bias subtraction, flat-field correction and cosmic-ray removal were performed with IRAF⁶⁰. The night conditions were not photometric and differential aperture photometry was also performed in comparison with nearby stars.

We divided the light curves into 20 flux bins and computed their cumulative distribution functions. As shown in Extended Data Fig. 6, the only apparent difference is between the Chandra and IR flux distributions at a mean-normalized flux of 0.7. However, owing to the different amplitude between the X-ray and O/IR flares, we also computed the distribution of the X-ray after applying a linear rescaling to make the X-ray minimum and maximum match the O/IR. A Kolmogorov–Smirnov test indicates that the GRS 1915+105 X-ray light curve and the Swift J1858 O/IR light curve are consistent with a single flux distribution. This analysis demonstrates that the O/IR variability patterns observed in August in Swift J1858 are commonly seen throughout the entire outburst and that they are compatible with the limit cycles in GRS 1915+105.

## Flux–flux diagram

We computed a flux–flux diagram using the time interval containing all bands from X-ray to IR (that is, from 500 s to 2,000 s). We binned the time series to 10 s (that is, the LT time resolution) and computed the correlation diagram between UV and the other bands. For the HAWK-I light curve, gaps were filled using a linear interpolation between the points at the edge, with the addition of Gaussian noise. Although the correlation between X-ray and UV is relatively poor ($\rho \approx 0.4$), the correlation between UV, optical and IR is clearly visible (see Extended Data Fig. 7a). The power-law index of the correlations were $\beta_{HAWK-I} = 0.3 \pm 0.02$ and $\beta_{LT} = 0.43 \pm 0.03$ for HAWK-I and LT, respectively. This indicates stronger variations in bluer bands⁴⁴, suggesting a thermal reprocessing scenario. This is confirmed when looking at constant IR to optical and UV to optical ratios (that is, the magnitude colour) for X-ray count rates >1 count s⁻¹ (see Extended Data Fig. 7b).

## Comparison with other sources

Previous studies indicated the presence of 'GRS 1915-like' patterns in different neutron stars (most notably the Bursting Pulsar³³, the Rapid Burster⁶ and the ULX NGC 3621 (ref. ³¹)). Our multiwavelength analysis would therefore indicate that similar instabilities would take place also in these sources. As shown for the case of the Bursting Pulsar GRO J1744–28 (ref. ⁶¹), the condition to show accretion disk instabilities is when the magnetospheric radius $R_m$ is larger than the star radius $R_*$. This means that the luminosity $L$

$$L > 2.75 \times 10^{38} \xi^{21/22} \alpha^{-1/11} m^{6/11} R_{*,6}^{7/11} B_{12}^{6/11} \text{ erg s}^{-1}, \quad (4)$$

in which $\xi$ is a factor between 0.5 and 1, $\alpha$ is the viscosity parameter, $m$ is the mass of the neutron star in units of $1.4 M_\odot$, $R_{*,6}$ its radius in units of $10^6$ cm and $B_{12}$ its magnetic field in units of $10^{12}$ G. On the other hand, when the magnetospheric radius is equal to or smaller than the stellar radius, the condition becomes:

$$L > 6.46 \times 10^{35} \alpha^{-1/8} m^{9/16} R_{*,6}^{5/16} \text{ erg s}^{-1}. \quad (5)$$

Finally, the magnetospheric radius is smaller than the star radius when

$$L > 2.83 \times 10^{45} \xi^{7/2} m^{1/2} R_{*,6}^{3/2} B_{12}^{2} \text{ erg s}^{-1}. \quad (6)$$

According to these equations, depending on their magnetic field, neutron stars with a luminosity of at least approximately $10^{36}$ erg s⁻¹ may show a radiation-pressure-dominated inner disk (see Extended Data Fig. 1). Therefore, we compiled a luminosity versus magnetic field diagram for a subset of sources to verify the self-consistency of our proposed scenario. As shown in Extended Data Fig. 1, all accreting neutron stars that show 'GRS 1915-like' repetitive patterns lie on the radiation-pressure-dominated disk region of the parameter space⁶,³¹,³³. The plot also shows the typical ranges of standard classes of accreting neutron stars.

At low magnetic fields ($<10^9$ G), we note that Swift J1858 lies in the radiation-pressure-dominated area and, as indicated by recent X-ray/radio studies³², is also between atolls and Z sources. The latter group is known to have clear evidence of outflows and radio ejecta: this strongly suggests that accretion instabilities may play a notable role in explaining the unique phenomenology of Z sources⁶²,⁶³. We also note that the Rapid Burster lies in the radiation-pressure-dominated regime⁶, even though the source does not always show such variability. Therefore, luminosity alone does not seem to be the only element necessary to trigger 1915-like patterns.

As for intermediate magnetic fields between $10^9$ and $10^{11}$ G, the Bursting Pulsar⁶⁴ is known to show strong repetitive burst patterns similar to GRS 1915+105 (ref. ⁶⁴). In this region of the parameter space, Terzan 5 X-2 also showed outflows close to its peak luminosity⁶⁵ (that is, in the radiation-pressure regime), but non-stationary variability patterns were not observed.

For sources with a very high magnetic field ($>10^{12}$ G), the radiation-pressure-dominated disk solution requires a luminosity $L_x > 5 \times 10^{38}$ erg s⁻¹, that is, almost in the ULXs regime. For such high accretion rates, it is argued that outflows are formed within the spherization radius⁶⁶,⁶⁷. A particularly interesting group of these sources is pulsating ULXs, which—owing to their proximity and high brightness—enable detailed studies of accretion in the super-Eddington regime. Spectral and temporal studies in the past 5 years have allowed us to obtain the magnetic field, as well as evidence for the presence of outflows in several sources⁶⁷⁻⁷¹, which define a clear region of the parameter space beyond the radiation-pressure-dominated disk transition. We also note that, as for the notorious Swift J0243.6+6124, some accreting pulsars can also behave as transient ULXs with luminosities well above the typical values of high-mass X-ray binaries⁷¹. In fact, such systems have been seen already in the Magellanic Clouds⁶⁹,⁷². For these nearby transient ULXs and many of the traditional ULXs, there is evidence of ionized emission and absorption that may be associated with outflows and disk winds⁷³.

The brightest source to show 'GRS 1915-like' patterns is the ULX NGC 3621 (ref. ³¹). However, the nature of its central object is still unknown. Therefore, even though the source is clearly in the radiation-pressure-dominated disk region, it is not possible to fully constrain its position in the $L$ versus $B$ diagram.

## Data availability

All raw data from the Swift J1858 August campaign and GRS 1915+105 Chandra/VLA observations are public and can be downloaded from their archives using the reported codes. All the reduced data from this campaign (including spectroscopic observations, which are not presented here) will be made object of a publication and made accessible
(Castro Segura et al. in preparation). Further analysis of the WASP and CHIMERA data is in progress, thus they are available on request to the authors.

35. Harrison, F. A. et al. The Nuclear Spectroscopic Telescope Array (NuSTAR) high-energy X-ray mission. Astrophys. J. 770, 103 (2013).

36. Pirard, J.-F. et al. in Ground-based Instrumentation for Astronomy, Vol. 5492 (eds Moorwood, A. F. M. & Iye, M.) 1763–1772 (SPIE, 2004).

37. Dhillon, V. S. et al. ULTRACAM: an ultrafast, triple-beam CCD camera for high-speed astrophysics. Mon. Not. R. Astron. Soc. 378, 825–840 (2007).

38. Eastman, J., Siverd, R. & Gaudi, B. S. Achieving better than 1 minute accuracy in the heliocentric and barycentric Julian dates. Publ. Astron. Soc. Pac. 122, 935 (2010).

39. McMullin, J. P., Waters, B., Schiebel, D., Young, W. & Golap, K. in Astronomical Data Analysis Software and Systems XVI, Vol. 376 (eds Shaw, R. A., Hill, F. & Bell, D. J.) 127–130 (Astronomical Society of the Pacific, 2007).

40. Gandhi, P. et al. Furiously fast and red: sub-second optical flaring in V404 Cyg during the 2015 outburst peak. Mon. Not. R. Astron. Soc. 459, 554–572 (2016).

41. Baglio, M. C. et al. A wildly flickering jet in the black hole X-ray binary MAXI J1535–571. Astrophys. J. 867, 114 (2018).

42. Edelson, R. A. & Krolik, J. H. The discrete correlation function: a new method for analyzing unevenly sampled variability data. Astrophys. J. 333, 646 (1988).

43. Gandhi, P. et al. Rapid optical and X-ray timing observations of GX 339–4: multicomponent optical variability in the low/hard state. Mon. Not. R. Astron. Soc. 407, 2166–2192 (2010).

44. Paice, J. A. et al. Blue oscillations and rapid red flares in Swift J1858.6-0814 observed with ULTRACAM/NTT. The Astronomer's Telegram No. 12197 (2018).

45. Janiuk, A. & Czerny, B. Time-delays between the soft and hard X-ray bands in GRS 1915 +105. Mon. Not. R. Astron. Soc. 356, 205–216 (2005).

46. O'Brien, K. et al. Echoes in X-ray binaries. Mon. Not. R. Astron. Soc. 334, 426–434 (2002).

47. Vincentelli, F. M. et al. Discovery of a thermonuclear Type I X-ray burst in infrared: new limits on the orbital period of 4U 1728-34. Mon. Not. R. Astron. Soc. 512, L37–L41 (2020).

48. Knight, A., Ingram, A. & Middleton, M. X-ray eclipse mapping constrains the binary inclination and mass ratio of Swift J1858.6–0814. Mon. Not. R. Astron. Soc. 495, 1908–1920 (2022).

49. Blumenthal, G. R. & Gould, R. J. Bremsstrahlung, synchrotron radiation, and Compton scattering of high-energy electrons traversing dilute gases. Rev. Mod. Phys. 42, 237–271 (1970).

50. Van der Laan, H. A model for variable extragalactic radio sources. Nature 211, 1131–1133 (1996).

51. Foreman-Mackey, D., Hogg, D. W., Lang, D. & Goodman, J. emcee: the MCMC Hammer. Publ. Astron. Soc. Pac. 125, 306 (2013).

52. Reid, M. J. et al. A parallax distance to the microquasar GRS 1915+105 and a revised estimate of its black hole mass. Astrophys. J. 796, 2 (2014).

53. Blandford, R. & Eichler, D. Particle acceleration at astrophysical shocks: a theory of cosmic ray origin. Phys. Rep. 154, 1–75 (1987).

54. Bell, A. R. The acceleration of cosmic rays in shock fronts – I. Mon. Not. R. Astron. Soc. 182, 147–156 (1978).

55. Markoff, S., Falcke, H. & Fender, R. A jet model for the broadband spectrum of XTE J1118+480. Synchrotron emission from radio to X-rays in the Low/Hard spectral state. Astron. Astrophys. 372, L25–L28 (2001).

56. Harmon, B. A. et al. Hard X-ray signature of plasma ejection in the galactic jet source GRS 1915+105. Astrophys. J. 477, L85–L89 (1997).

57. Belloni, T., Méndez, M., King, A. R., van der Klis, M. & van Paradis, J. A unified model for the spectral variability in GRS 1915+105. Astrophys. J. 488, L109–L112 (1997).

58. Harding, L. K. et al. CHIMERA: a wide-field, multi-colour, high-speed photometer at the prime focus of the Hale telescope. Mon. Not. R. Astron. Soc. 457, 3036–3049 (2016).

59. Nikzad, S. et al. High-efficiency UV/optical/NIR detectors for large aperture telescopes and UV explorer missions: development of and field observations with delta-doped arrays. J. Astron. Telesc. Instrum. Syst. 3, 036002 (2017).

60. Tody, D. in Instrumentation in Astronomy VI, Vol. 0627 (ed. Crawford, D. L.) 733 (SPIE, 1986).

61. Mönkönen, J. et al. Evidence for the radiation-pressure dominated accretion disk in bursting pulsar GRO J1744–28 using timing analysis. Astron. Astrophys. 626, A106 (2019).

62. Homan, J. et al. Evidence for simultaneous jets and disk winds in luminous low-mass X-ray binaries. Astrophys. J. 830, L5 (2016).

63. Motta, S. E. & Fender, R. P. A connection between accretion states and the formation of ultrarelativistic outflows in a neutron star X-ray binary. Mon. Not. R. Astron. Soc. 483, 3686–3699 (2019).

64. Miller, J. M., Maitra, D., Cackett, E. M., Bhattacharyya, S. & Strøhmayer, T. E. A fast X-ray disk wind in the transient pulsar IGR J17480–2446 in Terzan 5. Astrophys. J. 731, L7 (2011).

65. Poutanen, J., Lipunova, G., Fabrika, S., Butkevich, A. G. & Abolmasov, P. Supercritically accreting stellar mass black holes as ultraluminous X-ray sources. Mon. Not. R. Astron. Soc. 377, 1187–1194 (2007).

66. Middleton, M. J., Higginbottom, N., Knigge, C., Khan, N. & Wiktorowicz, G. Thermally driven winds in ultraluminous X-ray sources. Mon. Not. R. Astron. Soc. 509, 1119–1126 (2022).

67. Walton, D. J. et al. A potential cyclotron resonant scattering feature in the ultraluminous X-ray source pulsar NGC 300 ULX1 seen by NuSTAR and XMM-Newton. Astrophys. J. 857, L3 (2018).

68. Vasilopoulos, G. et al. NGC 300 ULX1: spin evolution, super-Eddington accretion, and outflows. Mon. Not. R. Astron. Soc. 488, 5225–5231 (2019).

69. Vasilopoulos, G. et al. The 2019 super-Eddington outburst of RX J2029.6-7427: detection of pulsations and constraints on the magnetic field strength. Mon. Not. R. Astron. Soc. 494, 5350–5359 (2020).

70. Kosec, P. et al. Ionized emission and absorption in a large sample of ultraluminous X-ray sources. Mon. Not. R. Astron. Soc. 508, 3569–3588 (2021).

71. Bachetti, M. et al. Orbital decay in M82 X-2. Astrophys. J. 937, 125 (2022).

72. Tsygankov, S. S., Doroshenko, V., Lutovinov, A. A., Mushtukov, A. A. & Poutanen, J. SMC X-3: the closest ultraluminous X-ray source powered by a neutron star with non-dipole magnetic field. Astron. Astrophys. 605, A39 (2017).

73. Shakura, N. I. & Sunyaev, R. A. Black holes in binary systems. Observational appearance. Astron. Astrophys. 24, 337–355 (1973).

74. Illarionov, A. F. & Sunyaev, R. A. Why the number of galactic X-ray stars is so small? Astron. Astrophys. 39, 185–195 (1975).

75. Campana, S. et al. The quiescent X-ray emission of three transient X-ray pulsars. Astrophys. J. 580, 389–393 (2002).

76. Homan, J. et al. Rossi X-ray Timing Explorer observations of the first transient Z source XTE J1701–462: shedding new light on mass accretion in luminous neutron star X-ray binaries. Astrophys. J. 656, 420–430 (2007).

77. Muñoz-Darias, T., Fender, R. P., Motta, S. E. & Belloni, T. M. Black hole-like hysteresis and accretion states in neutron star low-mass X-ray binaries. Mon. Not. R. Astron. Soc. 443, 3270–3283 (2014).

78. Gladstone, J., Done, C. & Gierliński, M. Analysing the atolls: X-ray spectral transitions of accreting neutron stars. Mon. Not. R. Astron. Soc. 378, 13–22 (2007).

79. Coburn, W. et al. Magnetic fields of accreting X-ray pulsars with the Rossi X-ray Timing Explorer. Astrophys. J. 580, 394–412 (2002).

80. Reig, P. Be/X-ray binaries. Astrophys Space Sci. 332, 1–29 (2011).

81. Kaaret, P., Feng, H. & Roberts, T. P. Ultraluminous X-ray sources. Annu. Rev. Astron. Astrophys. 55, 303–341 (2017).

**Acknowledgements** F.M.V. thanks R. Arcodia, P. Casella, G. Marcel, G. Mastroserio, N. Scepi and L. Stella for insightful discussions. The interpretation of the results benefited from discussions held during the meeting 'Looking at the disc-jet coupling from different angles' at the International Space Science Institute in Bern, Switzerland. F.M.V. was supported by the NASA awards 80NSSC19K1456, 80NSSC21K0526 and from grant FJC2020-043334-I financed by MCIN/AEI/10.13039/501100011033 and NextGenerationEU/PRTR. J.N. acknowledges support by the SAO award G01-22036X. A.J.T. is a NASA Einstein Fellow and acknowledges support for this work provided by NASA through the NASA Hubble Fellowship grant no. HST-HF2-51494.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc. for NASA, under contract NAS5-26555. D.A. and N.C.S. acknowledge support from the Science and Technology Facilities Council (STFC) grant ST/VO001000/1. F.M.V., M.A.P. and V.A.C. acknowledge support from the Spanish Ministry of Science and Innovation research project PID2020-120323GB-100. M.A.P. acknowledges support from the Consejería de Economía, Conocimiento y Empleo del Gobierno de Canarias and the European Regional Development Fund (ERDF) under grant with reference ProID2021010132 ACCISI/FEDER, UE. T.B. acknowledges financial contribution from the agreement ASI-INAF n.2017-14-H.0 and from PRIN-INAF 2019 N.15. T.M.D. acknowledges support from the Spanish Ministry of Science and Innovation project PID2021-124879NB-100 and the Europa Excelencia grant (EUR2021-122010). T.R. acknowledges the financial contribution from the agreement ASI-INAF n.2017-14-H.0.

**Author contributions** F.M.V. and J.N. drew the new physical scenario for multiwavelength instabilities. F.M.V. carried out the multi-λ timing analysis of the Swift J1858 dataset. A.J.T. modelled the radio variability with the help of S.d.P. and J.v.d.E. Y.C. modelled the phase dependence of the X-ray/IR lag with the help of N.C.S. For Swift J1858, N.C.S. provided the UV and X-ray data; F.J.-I. provided the optical data from LT; G.V., C.D.B., J.M. and M.O.A. provided optical data from Chimera and WASP; J.v.d.E. and T.R. provided the radio data. J.N. is the principal investigator of the Chandra and VLA proposal on GRS 1915; A.J.T. reduced and analysed the radio data. F.M.V., Y.C., G.V., D.A., T.B., N.D., T.M.D. and J.v.d.E. contributed substantially to the development of the luminosity–magnetic field diagram (Extended Data Fig. 1) to compare the different neutron stars. All authors contributed actively to the discussion and to the final version of the manuscript.

**Competing interests** The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41586-022-05648-3.

**Correspondence and requests for materials** should be addressed to F. M. Vincentelli.

**Peer review information** Nature thanks Emrah Kalemci and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permissions information** is available at http://www.nature.com/reprints.
<figure><img src="image_12.png" /></figure>

**Extended Data Fig. 1 | Luminosity versus magnetic field diagram for accreting neutron stars.** The long dashed line represents the limit in the parameter space for which the inner disk is (blue) or is not (red) radiation pressure dominated, that is, where the transition between zones A and B ($R_{AB}$) of a Shakura and Sunyaev$^{73}$ disk is greater than the magnetospheric radius ($R_m$) or the star radius ($R_*$), whichever is larger. Because the disk rotates with Keplerian velocity, for large $R_m$, the inner disk rotates slower than the neutron star magnetosphere and accretion can be halted owing to the propeller effect$^{74}$. The dashed grey lines correspond to the propeller threshold luminosities as a function of the magnetic field for different spins$^{74,75}$. For all these lines, we set $\xi = 1$, $\alpha = 0.1$, $M_{NS} = 1.4 M_\odot$ and $R = 10$ km. We then marked the sources that showed '1915-like' variability patterns (green ticks) in their light curves: the Rapid Burster$^6$, the Bursting Pulsar$^{61,64}$, the ULX NGC 3621 (ref.$^{31}$) and our target, Swift J1858. Along with these, we also showed the different phenomenological

classes of accreting neutron stars X-ray binaries, depending on their magnetic field and luminosities. At low magnetic field ($\le 10^9$ G), we find classical LMXBs, which, depending on their accretion rate, can manifest as atolls, bright atolls or Z sources$^{76,77}$. X-ray timing studies have also shown that accreting millisecond pulsars (AMXPs), with a magnetic field in the $10^8$ G to $10^9$ G range, are also compatible with atolls in their hard state$^{78}$. For higher magnetic fields of $10^9$ G, the observed accreting neutron stars usually show pulsations, but—owing to the lower propeller threshold—also have slower spin periods with respect to AMXPs. Above $10^{11}$ G, the diagram is mainly populated by high-mass X-ray binaries$^{79,80}$ and pulsating ULXs$^{81}$. For all these classes/objects, we also marked which of them show typical phenomenology linked to accretion instabilities, that is, radio ejecta (yellow star) and outflows (cyan waves). We note that these phenomena tend to be present above the radiation pressure disk threshold (see Methods).
<figure><img src="image_13.png" /></figure>

**Extended Data Fig. 2 | Infrared lag analysis.** a, Lag distribution from the five simultaneous NuSTAR-HAWK-I windows after $10^4$ flux randomizations. An evolution of the lag centroid is visible. b, CCF computed between HST and

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

HAWK-I. Excluding the asymmetry at longer lags, owing to the asymmetry of the flares, the CCF peaks at 0.
<figure><img src="image_15.png" /><figcaption>**Extended Data Fig. 3** | Averaged flare profile for events. Owing to the presence of two nearby flares, these are also present not at the centre. However, the overall connection between long and short timescales is still clear.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>c</figcaption></figure>

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

**Extended Data Fig. 4** | X-ray versus infrared lag modelling. **a**, Fit to the lags obtained including the flares in the light curves. The only parameter allowed to change is the inclination of the binary, *i*. Long and short dashed curves represent 68% and 99% confidence levels, respectively. **b**, Histograms of the posterior distributions of the inclinations. In all plots, dotted lines indicate the

<figure><img src="image_18.png" /><figcaption>d</figcaption></figure>

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

0.5th and 99.5th percentiles, dashed lines indicate the 16th and 84th percentiles and solid lines indicate the median. **c**, Fit to the lags measured excluding the flares from the light curves. **d**, Histogram of the posterior distribution excluding the flares.
<figure><img src="image_20.png" /><figcaption>**Extended Data Fig. 5** | Unconstrained radio modelling of Swift J1858. Same as Fig. 3a-c but modelling the Swift J1858 radio light curve using no constraints on the ejection times. We found similar results with larger errors.</figcaption></figure>
a

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

**Extended Data Fig. 6 | Comparison of the beats at different wavelengths.**

a, Mean-normalized light curves used for quantifying the association between Swift J1858 and GRS 1915: Chandra data from GRS 1915+105 (blue, +1 shift applied), HAWK-I data from Swift J1858 (red), WASP data from Swift J1858

b

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

(green, -1 shift applied) and CHIMERA data from Swift J1858 (purple, -2 shift applied). b, Cumulative distribution function (CDF) of the flux distributions of the light curves.
a

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

**Extended Data Fig. 7** | Flux–flux correlation diagram. **a**, Flux–flux diagram with respect to the HST measurements for NuSTAR (open circles), LT (open squares) and HAWK-I (filled circles). Although the O/IR is well correlated, the X-rays show a non-linear trend. All bands have been normalized to their

b

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

average. **b**, The plot shows the ratio of HAWK-I over LT (filled circles) and HST over LT (opened circles) as a function of the X-ray count rate normalized to 6.5 counts s⁻¹.
## Extended Data Table 1 | Swift J1858 jet modelling results without constraining the ejection time

<table><thead><tr><th>Ejection #</th><th>t<sub>eject</sub> (s)</th><th>F<sub>p</sub> (μJy)</th><th>β (%c)</th><th>R<sub>p</sub> (10<sup>12</sup> cm)</th></tr></thead><tbody><tr><td>1</td><td>2139<sup>+102</sup><sub>-167</sub></td><td>84<sup>+7</sup><sub>-8</sub></td><td>24 ± 2</td><td>0.9<sup>+0.2</sup><sub>-0.3</sub></td></tr><tr><td>2</td><td>4053<sup>+68</sup><sub>-40</sub></td><td>287<sup>+21</sup><sub>-30</sub></td><td>66<sup>+1</sup><sub>-2</sub></td><td>41.6 ± 0.5</td></tr><tr><td>3</td><td>6116<sup>+181</sup><sub>-146</sub></td><td>395<sup>+37</sup><sub>-51</sub></td><td>61<sup>+2</sup><sub>-4</sub></td><td>41.9 ± 0.6</td></tr><tr><td>4</td><td>8389<sup>+407</sup><sub>-200</sub></td><td>384<sup>+42</sup><sub>-71</sub></td><td>62<sup>+2</sup><sub>-5</sub></td><td>1.9<sup>+0.6</sup><sub>-0.7</sub></td></tr></tbody></table>

The third ejecta has values that markedly differ from the others (including very large errors; see also Extended Data Fig. 5), suggesting the need for some further kind of constraint. We further report the radii ($R_p$) at 7.5 GHz based on the best-fit model parameters for each ejecta. Here the subscript 'p' denotes values at the peak of the flare component at 7.5 GHz.
Extended Data Table 2 | Swift J1858 jet modelling results constraining the ejection time

<table><thead><tr><th>Ejection #</th><th>t<sub>eject</sub> (s)</th><th>F<sub>p</sub> (μJy)</th><th>β (%c)</th><th>R<sub>p</sub> (10<sup>12</sup> cm)</th></tr></thead><tbody><tr><td>1</td><td>2471<sup>+153</sup><sub>-203</sub></td><td>100<sup>+15</sup><sub>-14</sub></td><td>36<sup>+6</sup><sub>-7</sub></td><td>0.9<sup>+0.6</sup><sub>-0.4</sub></td></tr><tr><td>2</td><td>3772<sup>+218</sup><sub>-224</sub></td><td>318<sup>+45</sup><sub>-37</sub></td><td>70<sup>+4</sup><sub>-3</sub></td><td>1.6<sup>+1.1</sup><sub>-0.7</sub></td></tr><tr><td>3</td><td>6143<sup>+248</sup><sub>-337</sub></td><td>350<sup>+50</sup><sub>-35</sub></td><td>59<sup>+5</sup><sub>-3</sub></td><td>1.7<sup>+1.1</sup><sub>-0.7</sub></td></tr><tr><td>4</td><td>8511<sup>+440</sup><sub>-756</sub></td><td>331<sup>+98</sup><sub>-47</sub></td><td>60<sup>+8</sup><sub>-5</sub></td><td>1.6<sup>+1.3</sup><sub>-0.7</sub></td></tr></tbody></table>

The four events were constrained to take place during the IR beats. Here the subscript 'p' denotes values at the peak of the flare component at 7.5 GHz. The model also constrains the quiescent 7.5 GHz radio level (50<sup>+0.1</sup><sub>-0.3</sub> μJy), the quiescent radio spectral index (α=-0.92±0.05), the opening angle (φ = 3.3<sup>+0.1</sup><sub>-0.1</sub>), the jet inclination (i = 83.9<sup>+0.1</sup><sub>-0.4</sub>), the optical depth at 7.5 GHz (τ<sub>p</sub> = 1.55<sup>+0.04</sup><sub>-0.06</sub>) and the magnetic field strength (B<sub>p</sub> = 0.04±0.01G). The global values of the jet are consistent with the unconstrained fit.
Extended Data Table 3 | GRS 1915+105 jet modelling results

<table><thead><tr><th>Ejection #</th><th>t<sub>eject</sub> (s)</th><th>F<sub>p</sub> (μJy)</th><th>β (%c)</th><th>R<sub>p</sub> (10<sup>12</sup> cm)</th></tr></thead><tbody><tr><td>1</td><td>400<sup>+8</sup><sub>-11</sub></td><td>8.3 ± 0.1</td><td>14.7<sup>+0.1</sup><sub>-0.2</sub></td><td>1.2±0.03</td></tr><tr><td>2</td><td>1334<sup>+34</sup><sub>-36</sub></td><td>5.5 ± 0.1</td><td>11<sup>+0.1</sup><sub>-0.3</sub></td><td>1.0±0.03</td></tr><tr><td>3</td><td>4610 ± 10</td><td>14 ± 0.1</td><td>18 ± 0.1</td><td>15±0.04</td></tr><tr><td>4</td><td>7995<sup>+11</sup><sub>-10</sub></td><td>14 ± 0.1</td><td>15 ± 0.1</td><td>1.5±0.04</td></tr><tr><td>5</td><td>14074<sup>+12</sup><sub>-11</sub></td><td>15 ± 0.1</td><td>22 ± 0.1</td><td>1.6±0.04</td></tr><tr><td>6</td><td>15426 ± 16</td><td>9 ± 0.1</td><td>18 ± 0.1</td><td>1.2±0.04</td></tr><tr><td>7</td><td>16864<sup>+20</sup><sub>-15</sub></td><td>10 ± 0.1</td><td>16 ± 0.1</td><td>1.3±0.04</td></tr><tr><td>8</td><td>20579<sup>+39</sup><sub>-60</sub></td><td>10<sup>+0.2</sup><sub>-0.1</sub></td><td>18<sup>+0.2</sup><sub>-0.3</sub></td><td>1.3±0.04</td></tr><tr><td>9</td><td>21394<sup>+19</sup><sub>-18</sub></td><td>9 ± 0.1</td><td>15 ± 0.1</td><td>1.3±0.04</td></tr><tr><td>10</td><td>23740<sup>+35</sup><sub>-39</sub></td><td>12<sup>+0.2</sup><sub>-0.5</sub></td><td>16<sup>+0.3</sup><sub>-0.6</sub></td><td>1.4±0.06</td></tr><tr><td>11</td><td>27144<sup>+230</sup><sub>-241</sub></td><td>19<sup>+0.3</sup><sub>-0.7</sub></td><td>20<sup>+0.2</sup><sub>-1.2</sub></td><td>1.8±0.08</td></tr></tbody></table>

We also obtained the best-fit parameters for the quiescent 7.8 GHz radio level (0.63±0.03 mJy), the quiescent radio spectral index (α=-0.6±0.003), the opening angle (φ=7.1°±0.04), the jet inclination (i=62.7°±0.03) and the magnetic field strength (B<sub>p</sub>=0.05±0.002 mG). The optical depth at 7.8 GHz (τ<sub>p</sub>) was fixed to 1.6. We further report the component radii (R<sub>p</sub>) at 7.8 GHz based on the best-fit model parameters for each ejecta. Here the subscript 'p' denotes values at the peak of the flare component at 7.8 GHz.