Constraining the physical properties of large-scale jets from black hole X-ray binaries and their impact on the local environment with blast-wave dynamical models

Relativistic discrete ejecta launched by black hole X-ray binaries (BH XRBs) can be observed to propagate up to parsec-scales from the central object. Observing the final deceleration phase of these jets is crucial to estimate their physical parameters and to reconstruct their full trajectory, with implications for the jet powering mechanism, composition and formation. In this paper we present the results of the modelling of the motion of the ejecta from three BH XRBs: MAXI J1820+070, MAXI J1535$-$571 and XTE J1752$-$223, for which high-resolution radio and X-ray observations of jets propagating up to $\sim$15 arcsec ($\sim$0.6 pc at 3 kpc) from the core have been published in the recent years. For each jet, we modeled its entire motion with a dynamical blast-wave model, inferring robust values for the jet Lorentz factor, inclination angle and ejection time. Under several assumptions associated to the ejection duration, the jet opening angle and the available accretion power, we are able to derive stringent constraints on the maximum jet kinetic energy for each source (between $10^{43}$ and $10^{44}$ erg, including also H1743$-$322), as well as placing interesting upper limits on the density of the ISM through which the jets are propagating (from $n_{\rm ISM} \lesssim 0.4$ cm$^{-3}$ down to $n_{\rm ISM} \lesssim 10^{-4}$ cm$^{-3}$). Overall, our results highlight the potential of applying models derived from gamma-ray bursts to the physics of jets from BH XRBs and support the emerging picture of these sources as preferentially embedded in low-density environments.


INTRODUCTION
Relativistic jets appear as an ubiquitous feature among accreting black holes (BH) in the Universe, from supermassive BHs in active galactic nuclei (AGN) to stellar-mass BHs in galactic X-ray binaries (XRBs).Highly relativistic jets are also produced in energetic transients events, such as gamma-ray bursts (GRBs) and tidal disruption events (TDEs), a large fraction of which are believed to be powered by accreting BHs.The short timescales of evolution (days to weeks) and the relative proximity of BH XRBs make them ideal targets on which to study the properties of relativistic jets (Fender 2006;Romero et al. 2017), some of which appear to also be scale-invariant, and thus valid for all accreting BHs (Körd-⋆ E-mail: francesco.carotenuto@physics.ox.ac.uk † Former NASA Einstein Fellow ing & Falcke 2005).In BH XRBs, different jets are produced in different phases of the outburst (Corbel et al. 2004;Fender et al. 2004).Compact jets, causally connected to the accretion flow, are observed during the hard spectral state (see Remillard &McClintock 2006 andHoman &Belloni 2005 for a review on spectral states), as they emit self-absorbed synchrotron radiation, which dominates in the radio through near-infrared (Corbel et al. 2000;Fender 2001;Markoff et al. 2001;Corbel & Fender 2002;Russell et al. 2013).On the other hand, discrete ejecta are observed to be launched during transitions between the hard and the soft states, producing strong multi-wavelength flares during which the synchrotron emission is initially self-absorbed and then optically thin at radio wavelengths (e.g.Tetarenko et al. 2017).These components consist of bipolar blobs of plasma that travel away from the core, often at apparently superluminal speeds, and might be considered as less-relativistic analogs of what is observed in AGN (e.g.Marscher et al. 2002;Gómez et al. 2008).As of today, spatially resolved discrete ejecta have been observed with radio-interferometric observations in 15 sources: GRS 1915+105 (Mirabel & Rodríguez 1994), GRO J1655-40 (Hjellming & Rupen 1995;Tingay et al. 1995), Cyg X-3 (Mioduszewski et al. 2001), GX 339-4 (Gallo et al. 2004), XTE J1550-564 (Hannikainen et al. 2001;Corbel et al. 2002), XTE J1752-223 (Yang et al. 2010;Miller-Jones et al. 2011), H1743-322 (Corbel et al. 2005;Miller-Jones et al. 2012), XTE J1859+226 (Rushton et al. 2017), MAXI J1535-571 (Russell et al. 2019), V404 Cyg (Miller-Jones et al. 2019), MAXI J1820+070 (Bright et al. 2020;Espinasse et al. 2020), MAXI J1348-630 (Carotenuto et al. 2021), EXO J1846-031 (Williams et al. 2022), MAXI J1803-298 (Wood et al. 2023) and MAXI J1848-105 (Bahramian et al. 2023).This sample represents around 20% of the current population of confirmed and candidate BH XRBs, which are, however, all believed to produce jets (Tetarenko et al. 2016;Corral-Santana et al. 2016).In case of non-detection, this is likely due to source being too far, or with an unfavorable inclination angle (due to the effect of Doppler boosting) or to a failed transition from the hard to the soft state, which was found to happen, approximately, for a third of the observed outbursts (Alabarta et al. 2021).For sources that do display hard-to-soft state transitions, optically thin radio flares can be always detected with an adequate radio monitoring.In some cases, discrete ejecta can propagate up to parsec scales far from the core, displaying re-brightenings and deceleration phases likely due to the interaction with the interstellar medium (ISM), which also result in the production of broadband synchrotron radiation (radio to X-rays) from in-situ particle acceleration, up to TeV energies (Corbel et al. 2002(Corbel et al. , 2005;;Migliori et al. 2017;Espinasse et al. 2020;Carotenuto et al. 2021).
Despite the wealth of multi-wavelength observations collected over the recent years, multiple aspects related to the formation, evolution and overall physics of these jets remain unclear.For instance, the powering mechanism of the jets is still an open problem, as jets could be powered from the extraction of energy from a spinning BH (Blandford & Znajek 1977) or from its accretion disk (Blandford & Payne 1982), or from a combination of the two, as suggested by general relativistic magnetohydrodynamic (GRMHD) simulations (McKinney 2006) and recent Event Horizon Telescope observations probing a possible light spine vs. massive sheath jet structure (Janssen et al. 2021).The plasma composition, either baryonic or purely leptonic, is also unknown, and it is notoriously difficult to constrain as most jets display only a simple featureless synchrotron spectrum (Fender 2006).Moreover, while radio/infrared timing techniques are opening a new window on the physical parameters of compact jets (Casella et al. 2010;Tetarenko et al. 2019Tetarenko et al. , 2021;;Zdziarski et al. 2022), and recent results found evidence for a luminosity dependence of their properties (Prabu et al. 2023), we still lack precise constraints on the physical parameters of discrete ejecta, such as their mass, speed, energy and volume.In particular, a key open problem is the quantification of their total energy content.Measuring the jet's energy is of prime importance not only to estimate the balance between inflows and outflows in BH XRBs, but also because of the implications on the jet composition, powering mechanism and impact on the surrounding environment (e.g.Fender & Muñoz-Darias 2016).
The total energy (internal plus kinetic) of discrete ejecta can be estimated with different approaches.First, given the jet synchrotron emission, it is possible to infer the internal energy of the plasma that is required to produce the observed radiation by relying on the knowledge of the size of the emitting region and of the source distance, while assuming equipartition conditions (Longair 2011).The size of the emitting region can be most easily estimated by directly resolving the plasmon with radio or X-ray observations (e.g.Rushton et al. 2017).When this is not possible, the synchrotronemitting region size can be estimated through the detection of the radio spectral peak due to synchrotron-self absorption (e.g.Fender & Bright 2019 for BH XRBs), or it can be computed assuming a jet expansion speed and an ejection timescale (usually the duration of the rise of the radio flare at the jet's launch), although this may largely underestimate the jet's internal energy (Bright et al. 2020;Carotenuto et al. 2022).An additional way of measuring the jet size relies on simultaneous radio-interferometric observations of the ejecta at the same frequency, but with different angular resolutions probing different spatial scales.By measuring the percentage of flux resolved out between the observations, it is possible to infer the size of the emitting region and subsequently the jet internal energy (Bright et al. 2020).Alternatively, it is possible to identify jet-produced structures in ISM and then use them as calorimeters to measure the mechanical power, and consequently the kinetic energy that the jets need to deposit in those structures in order to create and sustain them (Gallo et al. 2005;Russell et al. 2007;Tetarenko et al. 2018Tetarenko et al. , 2020)).
Independently, focusing on the kinematics of these jets and covering their full trajectory allows us to obtain a complete dataset (angular separation vs. time) of their evolution, which can later be used to test physical models for the jet propagation in the ISM.The application of these models, which are mostly derived from the physics of GRBs (Wang et al. 2003), can yield important constraints on multiple physical parameters of the ejecta, such as their Lorentz factor, mass, inclination angle, ejection time and kinetic energy (Wang et al. 2003;Hao & Zhang 2009;Steiner & McClintock 2012).In particular, covering the jet deceleration phase with dense monitoring campaigns can significantly improve the constraints from these models (Carotenuto et al. 2022).Furthermore, modelling the jet motion is also of prime importance to precisely constrain the time of ejection, which is fundamental to put the jet launch in context with other multi-wavelength observational signatures, such as radio flares, X-ray spectral changes and X-ray quasi-periodic oscillations (QPOs, e.g.Ingram & Motta 2019).This allows us to ultimately obtain a comprehensive view of the source evolution during the state transition, with a special focus on the hot corona of electrons that surrounds the BH, which is responsible for the hard X-ray emission and it is thought to be intimately connected to the jet (e.g.Rodriguez et al. 2003;Markoff et al. 2005;Kara et al. 2019;Méndez et al. 2022;Ingram et al. 2023).
Tracing the jet motion has also turned out to be to especially useful to make use of the jets as probes of the environment surrounding BH XRBs.In fact, different works considering the propagation of ejecta in the ISM have provided strong evidence that BH XRBs are generally located in environments that appear 2-4 orders of magnitude less dense than the canonical Galactic ISM density of 1 particle per cm −3 (at least in the direction of the jet propagation), unless these jets are very narrowly collimated, with opening angles ≪ 1°, or extremely energetic, with kinetic energies above 10 46 erg (Heinz 2002;Hao & Zhang 2009;Carotenuto et al. 2022;Zdziarski et al. 2023).
The wealth of information that can be extracted from this type of modelling was first shown with the application on the large-scale decelerating jets from the BH XRBs XTE J1550-564 (Hao & Zhang 2009;Steiner & McClintock 2012), H1743-322 (Steiner et al. 2012), and, more recently, for MAXI J1348-630 (Carotenuto et al. 2022;Zdziarski et al. 2023).In this paper, as a continuation of the work started in Carotenuto et al. (2022), we expand the sample of sources that displayed unambiguously decelerating discrete ejecta, and for which such modelling has been applied, to include the jets from the BH XRBs MAXI J1820+070, MAXI J1535-571 and XTE J1752-223.These sources displayed resolved, large-scale decelerating jets observed between 2010 and 2018 that were observed to propagate up to ∼15 arcsec far from the core (Yang et al. 2010;Miller-Jones et al. 2011;Yang et al. 2011;Russell et al. 2019;Bright et al. 2020;Espinasse et al. 2020).However, the jet motion in these sources has only been described with basic phenomenological models, mostly applied in order to constrain the ejection date in relation with the simultaneous X-ray activity of the core.Since the quality of the jet angular separation data justifies the application of a physical model to describe the entire jet evolution, we performed such modelling and we present the detailed results in this paper.In particular, we present and discuss new constraints on the jet Lorentz factor, inclination angle, ejection time, as well as upper limits on the maximum energy available to the jets, on the density of the ISM that surrounds the systems and on the mass of the ejecta.For the last part, we also consider the ejecta launched in 2003 by H1743-322, where similar modelling has been already published by Steiner et al. (2012).
In Section 2 we present the sources and the observational data considered for the modelling work, while in Section 3 we discuss in detail the dynamical model that we adopted.Then, in Section 4 we present the results of the application of such model to our data and we discuss our findings in relation to the current understanding of jets from XRBs in Section 5. Finally, we summarize our conclusions in Section 6.

DATA
The data on the ejecta launched by the three BH XRBs considered in this paper have been already published and are therefore available in the literature.In the following sections, we present the sources and discuss the data used for this work.

MAXI J1820+070
The BH XRB MAXI J1820+070 was discovered by the Monitor of All-sky X-ray Image on board the International Space Station (Matsuoka et al. 2009) in March 2018 (Kawamuro et al. 2018), and it was subsequently identified with the optical transient ASASSN-18ey (Denisenko 2018).MAXI J1820+070 is one of the most wellobserved and well-studied BH XRBs in recent years.It harbors a 6.75 +0.64  −0.46 M⊙ BH accreting from a 0.5 ± 0.1 M⊙ companion star (Torres et al. 2019(Torres et al. , 2020;;Mikołajewska et al. 2022).Due to its impressive brightness, primarily in the X-rays (Fabian et al. 2020), it has been the subject of numerous multi-wavelength observing campaigns across the entire electromagnetic spectrum (e.g.Shidatsu et al. 2018;Hoang et al. 2019;Bright et al. 2020;Tetarenko et al. 2021;Abe et al. 2022;Cangemi et al. 2023;Echiburú-Trujillo et al. 2024) yielding a dataset of extremely high quality for a BH XRB in outburst.A model-independent measurement of the distance of 2.96 ± 0.33 kpc is available thanks to Very Long Baseline Interferometry (VLBI) radio parallax observations (Atri et al. 2020).
Bipolar relativistic discrete ejecta from MAXI J1820+070 have been detected and monitored at radio wavelengths for almost one year, with observations at different angular resolutions with the Multi-Element Radio Linked Interferometer Network (eMERLIN), the Very Long Baseline Array (VLBA), the Arcminute Microkelvin Imager Large Array (AMI-LA), the Karl G. Jansky Very Large Array (VLA) and the MeerKAT radio interferometer, showing the jets to propagate out to ∼10 arcsec from the core of the system with a high proper motion (Bright et al. 2020;Wood et al. 2021).Notably, these jets have also been detected at large scales, up to 12 arcsec, in the X-rays with five Chandra X-ray telescope exposures between the end of 2018 and 2019 (Espinasse et al. 2020).These X-ray detections are particularly important because they cover the deceleration phase, not immediately evident from the radio data alone.In this work, we use both the radio and X-ray coordinates of the jets to model their motion, and we also take into account the updated jet coordinates from Wood et al. (2021), obtained with the application of the new dynamic phase centre tracking technique to the VLBA data (see also Wood et al. 2023).

MAXI J1535-571
The BH XRB MAXI J1535-571 was discovered by MAXI in September 2017 (Negoro et al. 2017) when it entered into outburst, and it was subsequently monitored at all wavelengths between radio and the hard X-rays during its 1-year long outburst (e.g.Tao et al. 2018;Huang et al. 2018;Russell et al. 2019;Parikh et al. 2019;Bhargava et al. 2019;Baglio et al. 2018).In particular, the full outburst evolution with the associated state transitions is discussed in Tao et al. (2018) and Nakahira et al. (2018).The source is located at a distance D = 4.1 +0.6  −0.5 kpc (Chauhan et al. 2019), determined from observations of Hi absorption carried out with the Australian Square Kilometre Array Pathfinder (ASKAP).
The radio monitoring campaign presented in Russell et al. ( 2019) covered the evolution of the jets from MAXI J1535-571 throughout the whole outburst, with Australia Telescope Compact Array (ATCA) and MeerKAT observations.Compact jets were detected during an initial brightening in the first hard state and they were subsequently observed to quench as the source transitioned to the intermediate state, displaying intense flaring activity (see also Russell et al. 2020).During the hard-to-soft state transition, MAXI J1535-571 launched a fast single-sided discrete jet that was detected and monitored with MeerKAT and ATCA for almost one year.The relativistic components was observed to propagate and decelerate up to an angular distance of ∼15 arcsec (Russell et al. 2019).Its monitoring allowed the authors to place model-independent constraints on the jet speed, inclination angle and ejection date, which we take into account in the modelling presented in this work.

XTE J1752-223
XTE J1752-223 is a BH XRB discovered by the Rossi X-ray Timing Explorer in 2009 (Markwardt et al. 2009) that remained active for almost one year in outburst and that has been the subject of dense multi-wavelength observing campaigns, mostly focused in the radio and X-ray bands (Shaposhnikov et al. 2010;Ratti et al. 2012;Brocksopp et al. 2013).A recent estimation based on the Bayesian analysis of the soft spectral state and the hard-to-soft state transition yields the following constraint on the source distance: D = 7.11 +0.27  −0.25 kpc (Abdulghani et al. 2024), which is notably more than twice the first distance estimation of 3.5 kpc from Shaposhnikov et al. (2010), and it is consistent with another recent estimation of D = 6 ± 2 kpc based on Gaia DR3 (Fortin et al. 2024).We note that Abdulghani et al. (2024) also provide a first BH mass estimation of 12 ± 1M⊙, based on the same method.In this paper, we adopt 7.1 kpc as the source distance, noting that adopting the 6 kpc value would not substantially change the main conclusions of this work.Also, we do not use the 3.5 kpc distance, as it is obtained with an uncertain X-ray spectral-timing correlation scaling technique based on the evolution of the photon index and the QPO frequency during the outburst (Shaposhnikov & Titarchuk 2009).
Notably, during the outburst and after the first hard state, the source performed a standard transition to the intermediate and then soft state, but then displayed multiple short-lived returns to the intermediate state accompanied by strong radio flaring activity observed with ATCA, implying the production and launch of multiple ejecta (Brocksopp et al. 2013).These ejecta, at least three separated approaching components, were eventually imaged with the European VLBI Network (EVN) and with VLBA (Yang et al. 2010;Yang et al. 2011;Miller-Jones et al. 2011), appearing to propagate only at small scales, i.e. less than one arcsec.At least one component (labeled "A" in Yang et al. 2010) displayed evidence of deceleration (Miller-Jones et al. 2011), while no receding component was detected and no ejecta were detected at larger distances from the core with interferometers probing larger angular scales (such as ATCA).In this paper we will only focus on this decelerating component, as not enough data are available for the other ejecta.

THE DYNAMICAL MODEL
We adopt a numerical blast-wave dynamical model to describe propagation of jets in the ISM.The model was originally developed to describe GRB afterglows (Piran 1999;Huang et al. 1999), and has been applied to describe the evolution of mildly relativistic ejecta in BH XRBs for the first time in Wang et al. (2003), including the transition from relativistic to non-relativistic motion.In particular, we consider the same implementation of Carotenuto et al. (2022).
The model considers a symmetric pair of confined conical ejecta launched simultaneously in opposite directions, at an inclination angle θ with respect to the line of sight.The ejecta start to move away from the core with an initial Lorentz factor Γ0 and kinetic energy E0, and expand with a constant half-opening angle ϕ inside an ambient medium with constant density nISM.Matter in the same ambient medium is entrained by the jets, which are therefore continuously decelerating, and during this process their kinetic energy is converted into internal energy of the swept-up ISM through external shocks, in a similar fashion as GRB afterglows (Wang et al. 2003).In particular, a forward shock develops at the contact discontinuity between the jet and the ISM.The shock continuously heats the encountered ISM and randomly accelerates particles.In this context, the radiated energy is assumed to be negligible and the jet expansion is considered adiabatic throughout its whole evolution (Chiang & Dermer 1999;Huang et al. 1999;Wang et al. 2003;Hao & Zhang 2009), an assumption that has been proven to be robust in recent works (Steiner & McClintock 2012;Bright et al. 2020;Carotenuto et al. 2022;Zdziarski et al. 2023).Considering one of the two ejected components, it is possible to write the energy conservation equation where the two terms on the right-hand side are, respectively, the instantaneous kinetic energy of the ejecta and the internal energy of the swept-up ISM.More in detail, Γ is instantaneous jet bulk Lorentz factor, M0 is the mass of the ejecta and Γ sh is the Lorentz factor of the shocked ISM.Here, σ is a numerical factor equal to 6/17 for ultrarelativistic shocks and ∼0.73 for non-relativistic shocks (Blandford & McKee 1976), and it is possible to interpolate between the two regimes with the following numerical scaling (Huang et al. 1999;Wang et al. 2003;Steiner & McClintock 2012): where β = (1 − Γ −2 ) 1/2 is the intrinsic jet speed in units of c.The mass of the entrained material msw can be written as where R is the distance from the compact object, n is the ambient density and mp is the proton mass.The Lorentz factor of the shocked ISM can be expressed as a function of the jet bulk Lorentz factor by using the jump conditions for an arbitrary shock (Blandford & McKee 1976;Steiner & McClintock 2012): where γ is the adiabatic index that varies between 5/3 for ultrarelativistic shocks and 4/3 for non-relativistic shocks.As the jet decelerates, we interpolate between the two limits with The relativistic kinematic equations for the approaching and receding components are (Rees 1966;Blandford et al. 1977;Mirabel & Rodríguez 1994): where the ∓ refers, respectively, to the approaching and receding jet, and t is the arrival time of the photons at the observer.
The measurable projected angular separation from the core is where D is the source distance.
In total, the model depends on 7 parameters: the jet initial kinetic energy E0 and Lorentz factor Γ0, the inclination angle of the jet axis θ and the jet half-opening angle ϕ, the source distance D, the external ISM density nISM and the ejection time tej.It is crucial to note that a degeneracy exists in this model between the three parameters E0, ϕ and nISM, as they appear as a single term in Equation 1 (taking into account the expression for msw in Equation 3).Hence, only the factor E0/nISMϕ 2 can be independently constrained by the application of this model.Similarly to Steiner & McClintock (2012), we refer to such factor as "effective energy", which here we define as Therefore, in order to obtain a reliable estimate of the jet kinetic energy, one needs to independently measure the two parameters ϕ and nISM, or to assume reasonable values for them (see Section 5.4).
For every set of parameters that compose the model, it is possible to obtain the proper motion curve of the jet on the plane of the sky and then to compare it with the data.This can be done by integrating Equation 6 starting at a time tej from an assumed distance R0 = 10 8 cm (from which there is only a weak dependence), and then numerically solving Equation 1 at every time step for the instantaneous jet Lorentz factor.The information on the instantaneous speed is then used to update the distance traveled by the jet, which is converted to the angular separation α (Equation 7), that can be directly compared to the observational data.

Fit setup
We fit the data for the three BH XRBs considered in this work with the dynamical model presented in Section 3. We adopt a Bayesian approach, applying a Monte Carlo Markov Chain (MCMC) code implemented with the emcee package (Foreman-Mackey et al. 2013).
For every point of the parameter space, Equation 6 was integrated using odeint from the SciPy package (Virtanen et al. 2020).
We include the maximum amount of available information in the choice of our priors, which are physically motivated from our knowledge of the source in question and of BH XRBs in general.We discuss the specific choices in the following sections dedicated to each source.Every MCMC run was conducted using 110 walkers.For each run, after manual inspection, we consider that convergence is reached when the positions of the walkers in the parameter space are no longer significantly evolving.Once the chains have converged, the best fit result for each parameter is taken as the median of the one-dimensional posterior distribution obtained from the converged chains, while the 1σ uncertainties are reported as the difference between the median and the 15th percentile of the posterior (lower error bar), and the difference between the 85th percentile and the median (upper error bar).

MAXI J1820+070
We first consider the bipolar ejecta from MAXI J1820+070.The angular separation for the two components is shown in Figure 1, including the measurements presented in Bright et al. (2020), Espinasse et al. (2020) and Wood et al. (2021).The approaching and receding components are marked, respectively, by red and blue points.We performed a joint fit of the dynamical model presented in Section 3 to the approaching and receding components.We adopted a flat prior for Γ0 (between 1 and 100) and a log-flat prior for Ẽ0 (between 10 35 and 10 55 erg).We further assumed a normal prior for the source distance centered on D = 2.96 kpc and with a width of 0.3 kpc (Atri et al. 2020), and we assumed a flat prior for tej centered on the ejection time of component C (MJD 58305.95)presented in Wood et al. (2021) and ranging between MJD 58300 and 58310.For the inclination angle, again relying on Wood et al. (2021), we used a normal distribution centered on 64°and with a width of 5°, while truncated outside the interval 0°-90°.
The best fit is shown in Figure 1, along with the proper motion of the two jet components, and the results are reported in Table 1.The statistical uncertainty range on the plot is represented as the ensemble of trajectories corresponding to the final positions of the walkers in the parameter space.From Figure 1, it is possible to see that the model fits exceptionally well to the data, and the agreement with observations can be seen from the residuals on the bottom panel of the same figure.The deceleration of both jets can be adequately described by a single Sedov phase in a homogeneous environment.This type of deceleration has also been modeled using a simple polynomial fit in Espinasse et al. (2020) and Wood et al. (2021), but we note that in our case the whole jet motion can be described by a single physical model.The statistical uncertainty on the fit is remarkably small thanks to the fact that we detected both components and that we had VLBI observations taken at the beginning of the jet motion (Bright et al. 2020;Wood et al. 2021), which allowed us to constrain with great accuracy the ejection time.The high-resolution Chandra observations at the end of the monitoring are also important to cover the deceleration phase (Espinasse et al. 2020).According to this model, the jet is launched at tej = MJD 58305.96+0.02 −0.02 with a bulk Lorentz factor Γ0 = 2.61 +0.54  −0.39 , an effective energy of Ẽ0 = 2.6 +0.4  −0.4 × 10 46 erg, and a medium-to-high inclination angle θ = 59.6°+ 1.2°− 1.0°.The source distance is D = 2.96 +0.11 −0.13 kpc, which tracks the prior choice based on the the radio parallax measurement by Atri et al. (2020).The posterior distributions for the parameters of the model are shown the Appendix A (Figure A1), where we present the corner plot displaying the one-dimensional posterior distribution for all the parameters and the two-parameters correlations.

MAXI J1535-571
We fit the external shock model to large scale jet data from MAXI J1535-571, using the measurements reported in Russell et al. (2019).Similarly to the fit already presented in Carotenuto et al. (2022), we fit the data for the approaching ejection only, given the non-detection of the receding counterpart.The associated proper motion is shown in Figure 2. As for MAXI J1820+070, we adopted a flat prior for Γ0 (between 1 and 100) and log-flat prior for Ẽ0 (between 10 35 and 10 55 erg).We assumed a flat prior for tej between MJD 58005 and 58025 and a normal prior for the source distance centered on D = 4.1 kpc, with a width of 0.5 kpc (Chauhan et al. 2019).For the inclination angle, we used a uniform distribution in cos θ truncated outside the interval 0°-45°, following the constraints reported in Russell et al. (2019).
The best fit results are reported in Table 1 and are shown in Figure 2, from which can be seen that the model fits the data remarkably well.The bottom panel of the figure displays the residuals, revealing a good agreement with the observations.Also in this case, the jet deceleration can be accurately described by a single Sedov phase in a homogeneous environment.From the fit, we can place constraints on the jet ejection date tej = MJD 58017.4+4.0 −3.8 , on its bulk Lorentz factor Γ0 = 1.6 +0.2 −0.2 and on its medium-to-low inclination angle θ = 30.3°+ 6.3°− 6.3°.The source distance is D = 4.2 +0.8 −0.9 kpc, which, also in this case, tracks the prior choice based on the the Hi absorption measurement by Chauhan et al. (2019).We also constrain the effective energy of the jets to be E0 = 5.8 +16.6  −4.0 × 10 48 erg.The posterior distributions for the parameters of the model are shown in Figure A2.

XTE J1752-223
Finally, we fit the dynamical model to the VLBI data of the approaching ejecta launched by XTE J1752-223, as reported in Miller-Jones et al. (2011).As for the previous cases, we adopted a flat prior for Γ0 (between 1 and 100) and a log-flat prior for Ẽ0 (between 10 35 and 10 55 erg).We assumed a uniform distribution in cos θ, with a truncation outside the interval 0°-45°, consistent with the constraints reported in Miller-Jones et al. (2011).Moreover, we assume a normal prior for the source distance centered on D = 7.1 kpc and with a width of 0.3 kpc, from Abdulghani et al. (2024).Given that only four data points are available, we fixed the ejection date tej to MJD 55217, just before the 20 mJy peak of the radio flare observed with ATCA (Brocksopp et al. 2013), and one day before the transition from the hard-intermediate state (HIMS) to the soft-intermediate state (SIMS) occurred around MJD 55218 (Shaposhnikov et al. 2010).
The best fit results are reported in Table 1 and  −0.6 × 10 45 erg, at a medium-to-low inclination angle of θ = 18.4°+ 2.5°− 2.3°.Given the scarcity of data on this source, we are unable to provide a robust value for the jet initial Lorentz factor Γ0, but from the posterior we can constrain Γ0 > 3.4 (99.7% of confidence), according to the choice of ejection date discussed above.As shown in Figure A3, the data provide a median value of ≃5.4,but this constraint depends directly on the chosen tej.In the same way, the lower limit on Γ0 can be relaxed if we assume an earlier  2021).The un-shaded, gray and seashell regions mark periods in which the source was, respectively, in the hard, intermediate and soft state (Shidatsu et al. 2018).The black horizontal dashed line represents the zero separation from the core, while the black continuous line represents the best fit obtained with the external shock model.The orange shaded area represents the total uncertainty on the fit and it is obtained by plotting the jet trajectories corresponding to the final positions of the MCMC walkers in the model parameter space.Residuals ([data -model]/uncertainties) are reported in the bottom panel.The model appears to provide an excellent description of the motion of both the approaching and receding ejecta, with a low statistical uncertainty.ejection date.For illustration, performing the same fit, but assuming tej = MJD 55210 or 55215 leads to Γ0 > 1.8 (median value ≃ 2.1) and Γ0 > 2.9 (median value ≃ 4.2), respectively.On the other hand, for ejections at later times, fixing tej = MJD 55218 at the transition from the HIMS to the SIMS (Shaposhnikov et al. 2010) results in an equally acceptable fit, yielding similar constraints for the initial jet Lorentz factor: Γ0 > 3.4 (99.7% of confidence, but with a median value ≃ 7.4).Similarly, fixing tej = MJD 55221 at the transition from the SIMS to the soft state (an ejection time which we deem more unlikely for BH XRBs) yields Γ0 > 3.5 (99.7% of confidence), with an extremely high median value of ≃ 13.4.We note that, for our initial Lorentz factor posterior distributions, fixing tej in the range MJD 55217-55821 results in very similar lower limits on Γ0, while the median value is much more sensitive to the choice of tej.

DISCUSSION
We successfully modeled the motion of the large-scale jets from three BH XRBs, MAXI J1820+070, MAXI J1535-571 and XTE J1752-223 with a dynamical blast-wave model based on external shocks.This physical model is found to provide an excellent description of the propagation of the ejecta in the ISM, regardless of whether we detect or not both ejected components, and it also allows us to place meaningful constraints on various parameters of the ejecta.After the application of this model to the jets of XTE J1550-564, H1743-322 and MAXI J1348-630 (Steiner & McClintock 2012;Steiner et al. 2012;Carotenuto et al. 2022), and considering the results shown in the previous section, it appears that all the large-scale jets display a deceleration consistent with a Sedov phase.The goodness of the fits obtained for all the sources confirm the validity of the application of physical models derived from our knowledge of GRBs to XRBs (e.g.Wang et al. 2003), and highlights the potential of using the largely developed set of theoretical GRB models (including their entire multi-wavelength emission) to explain even more features observed in the less-relativistic jets from XRBs, as for instance the presence of forward and reverse shocks within the jet.We discuss in the following sections the constraints on the jet physical parameters and on the source environment that we obtained in this work, comparing them with our current knowledge of jets from BH XRBs.

Lorentz factor
We first discuss the constraints on the initial Lorentz factor Γ0 for the ejecta in our sample.It is generally difficult to constrain this parameter from the simple observation of the jet propagation, especially if the ejecta are significantly superluminal.The reason is that a source of significantly relativistic jets (with proper motions µapp and µrec) will usually be observed close to its maximum allowed distance Dmax = c(µappµrec) −1/2 , where Γ0 diverges (Fender et al. 1999;Fender 2003).Therefore, only lower limits on Γ0 are available for most of the sources displaying discrete ejecta (e.g.Fender 2003;Miller-Jones et al. 2006;Bright et al. 2020).
For MAXI J1820+070, we obtain an interesting estimate of the initial Lorentz factor of the jets: Γ0 = 2.6 +0.5 −0.4 , which implies a mildly relativistic ejecta.Such constraint is consistent from the previous lower limit Γ0 > 2.1 (Bright et 2011) and information on the spectral states from Shaposhnikov et al. (2010) and Brocksopp et al. (2013).Once the ejection date is fixed (here at MJD 55217), the model with a single Sedov phase appears to fit reasonably well the data.
Table 1.Parameters of the blast-wave model applied in this work inferred from the Bayesian fit described in Section 4. The values quoted are the median parameter and the 1σ confidence intervals derived from the MCMC run.The effective energy is defined as Ẽ0 = E 0 /n ISM ϕ 2 (see Section 3, Equation 8).
† : The constraint on Γ 0 for XTE J1752-223 strongly depends on the preferred ejection date (see Section 4.4).and, interestingly, this determination has been possible despite the fact that MAXI J1820+070 is located close to its Dmax = 3.1 kpc (Wood et al. 2021).

Parameter
In the case of MAXI J1535-571, we are also able to place an important constraint on the initial Lorentz factor of the approaching component Γ0 = 1.6 +0.2 −0.2 , implying a relatively slow ejecta, traveling initially at ∼0.77c.Jets from MAXI J1535-571 appear to be among the less relativistic in the observed sample of discrete ejecta (Miller-Jones et al. 2006;Steiner & McClintock 2012;Steiner et al. 2012;Carotenuto et al. 2022).Interestingly, the ejecta from MAXI J1535-571 is the one that propagates to the largest distance from the core (up to 0.8 pc).Launched with a high Ẽ0 (see Table 1), it is likely among the most energetic and the least relativistic jets observed so far, similar in nature to MAXI J1348-630, which also displayed large scale ejecta propagating up to 0.6 pc from the core, with Γ0 = 1.85 +0.15  −0.12 (Carotenuto et al. 2021(Carotenuto et al. , 2022)).Since it appears that jets from MAXI J1348-630 and MAXI J1535-571 are among the most energetic observed so far, this likely implies that the mass content of the ejecta is probably the driving factor that determines the large distance at which the ejecta propagate, with M0 being the dominant factor in the (Γ0−1)M0c 2 kinetic energy term in Equation 1 (see also the discussion in Zdziarski & Heinz 2024).At the same time, again similarly to MAXI J1348-630, such determination of Γ0 is one of the most precise and robust (although model-dependent) constraints on the Lorentz factor of a jet from a BH XRB to-date, with this being due to the ideal combination of low Γ0 and low inclination angle, which is less affected by the degeneracy between θ and the source distance (Fender 2003, Fender et al., in prep.).
We can only provide a lower limit Γ0 > 3.4 for the ejecta of XTE J1752-223, which, however, directly depends on the chosen tej, and it can be relaxed if we assume an earlier ejection date, as mentioned in Section 2.3.On the other hand, if tej is fixed to closer to the state transition, our robust lower limit on Γ0 does not vary (see Section 4.4).However, if the ejecta were launched on MJD 55218, at the peak of the radio flare on the day of the HIMS-to-SIMS transition, they would have a most likely Γ0 ≃ 7.4, which would be the highest Lorentz factor ever observed for these objects and it would challenge the common assumption that jets from BH XRBs are only mildly relativistic, unlike what is observed in AGN and GRBs (e.g.Jorstad et al. 2005;Ghirlanda et al. 2018).Observations of the ejecta closer to the core would have greatly helped to determine the initial Lorentz factor of this source, which is likely to be higher than the average among the available sample of ejecta.Overall, the data available so far appear to suggest that BH XRBs are able to accelerate relativistic jets in the mildly relativistic range 1 ≲ Γ0 ≲ 2, but a growing number of possible exceptions that suggest even faster jets (as for instance MAXI J1820+070 and XTE J1752-223), and this appears to be similar to the range of generally inferred bulk Lorentz factor for compact jets (e.g.Casella et al. 2010;Saikia et al. 2019;Tetarenko et al. 2019Tetarenko et al. , 2021;;Zdziarski et al. 2022).

Inclination angle
The inclination angle that we obtain for MAXI J1820+070 is θ = 59.6°+ 1.0°− 1.2°.We note that our posterior is still consistent with the chosen normal prior distribution θ = 64°± 5°from Wood et al. (2021).However, the peak of the posterior is slightly lower than the peak of the prior, which is based on the first part of the jet motion.With the values obtained for Γ0 and θ, we can compute the Doppler factor for the two ejecta.The Doppler factor for the approaching component at launch is δapp = Γ −1 0 (1 − β0 cos θ) −1 ≃ 0.7, while for the receding is δrec = Γ −1 0 (1+β0 cos θ) −1 ≃ 0.25.This implies that at the beginning of their motion both components are Doppler de-boosted, and their intrinsic luminosity is decreased, respectively by a factor δ 3−α app ≃ 0.3 and δ 3−α rec ≃ 0.008, using the formalism for discrete jet components, where the flux density follows Sν ∼ ν α (e.g.Blandford et al. 1977) and a spectral index α = −0.6 (Espinasse et al. 2020).
We constrain the inclination angle of the jet axis in MAXI J1535-571 to be θ = 30.3°+ 6.3°− 6.3°.Already, from the simple measurement of the jet proper motion and the source distance, Russell et al. (2019) strongly constrained the maximum jet inclination to θ < 45°(which is included in our choice of the prior on θ).Notably, such a value does not appear to be consistent with the inclination angle of the inner edge of the accretion disk obtained from NICER X-ray observations Miller et al. (2018).The authors report an angle i = 67.4°±0.8°f rom the spectral fitting of the relativistic reflection component in the high-SNR X-ray spectrum obtained during the intermediate state (they also report a near maximal BH spin of a = 0.994 ± 0.002, Miller et al. 2018).At the same time, our value is broadly consistent with the inclination angle (i = 37°+ 22°− 13°) of the region emitting a narrow and asymmetric iron line (Miller et al. 2018).Interestingly, the authors explain such difference in i with potential disk warping.The results of this comparison appear counter-intuitive, as we would generally expect jets to be launched along the direction of the BH spin, which, in turn, should be aligned with the inner edge of the accretion disk (Bardeen & Petterson 1975).However, at least a case in which jets were launched along the axis of a rapidly precessing inner disk has been observed (Miller-Jones et al. 2019).As in the case of MAXI J1348-630 (Carotenuto et al. 2022), a measurement of the orbital plane of the system would be useful to test the potential 3D alignment between the disk and the jet axis, taking into account that cases of disk/jet misalignment have been observed (Miller-Jones et al. 2019;Poutanen et al. 2022) and this is also supported by Table 2. Time delay ∆t ej,X between the inferred ejection times t ej and the possibly associated observed X-ray signatures.Here a positive ∆t ej,X means that the the jet is launched before the first appearance of the X-ray signature.The uncertainties on ∆t ej,X are the same as the ones obtained for t ej in this work and in Carotenuto et al. (2022).For MAXI J1820+070, we include both Components A (slow jet) and C (fast jet) as labeled in Wood et al. (2021), reminding the reader that in this work we only focus on Component C. The uncertainty on MAXI J1535-571 is due to the large uncertainty on t ej and to the fact that Type-B QPOs might have been present between MJD 58016.8 and 58025 (Stevens et al. 2018).GRMHD simulations (e.g.Liska et al. 2018).The Doppler factor for the approaching component at launch is δapp ≃ 1.9, implying a Doppler boosting with a factor δ 3−α app ≃ 10, while for the receding component we can compute δrec =≃ 0.37, implying a de-boosting of a factor δ 3−α rec ≃ 0.03, again adopting α = −0.6.Such a high Doppler de-boosting easily explains the non-detection of the receding component, as the radio flux density was pushed below the ATCA and MeerKAT sensitivity limit for the available exposure times (Russell et al. 2019).

Ejection time
The ejection time is a crucial piece of the puzzle in the current effort to reconstruct the precise sequence of events that lead to the formation and launch of discrete ejecta, which is still unclear, and modelling the jet motion is a reliable way of obtaining such information.
In the case of MAXI J1820+070, we infer an ejection date tej = MJD 58305.96+0.02 −0.02 , which is completely consistent with the most up-to-date estimation reported in Wood et al. (2021).We note that the quality of the data from Bright et al. (2020), especially the dense VLBI monitoring at the hard-to-soft state transition, and the additional significant improvements obtained with the new dynamic phase-center tracking technique adopted in Wood et al. (2021), already allowed the authors to constrain the ejection time with great accuracy, and it is worth noting that a precision of roughly 30 minutes has rarely been achieved for this type of events.It is also worth nothing that the obtained tej places the jet launch close to the peak of the radio flare observed at 15.5 GHz with AMI-LA (Bright et al. 2020;Homan et al. 2020).Moreover, this jet appeared to be launched approximately 6 hours after the first detection of Type-B QPOs in this outburst, defining the beginning of the soft-intermediate state (Homan et al. 2020).However, it is important to remark that Wood et al. ( 2021) associate the radio flare and the detection of Type-B QPOs with the ejection of a different pair of ejecta (with the approaching jet labeled as Component A), which had an intrinsic speed β ∼ 0.3, were not detected beyond the milli-arcsec scale and were ejected ∼9 h before the fast large-scale jets that we consider in this work (Wood et al. 2021).Specifically, Component A displayed an elongated structure and its ejection was inferred to last ∼6 h, hence partially overlapping with the first appearance of Type-B QPOs in this source.Interestingly, Wood et al. (2021) do not identify any Xray or radio flare counterpart to Component C. It has been suggested that Type-B QPOs could correspond to the time of jet launching (e.g.Fender et al. 2009;Miller-Jones et al. 2012), implying a strong causal relation between the two phenomena that has been particularly highlighted in Homan et al. (2020).However, for MAXI J1820+070 such link appears to be stronger with Component A than with Component C. Therefore, if this timing signature is indeed linked to jet ejections, it is unclear at the moment whether there is any connection with the fast ejecta observed to propagate at large scales.Such connection is also not confirmed for other sources (Miller-Jones et al. 2012;Russell et al. 2019;Carotenuto et al. 2021), for which in general the ejections are inferred to happen from hours to days before the detection of Type-B QPOs, as can also be seen in Table 2 (for which we include data from MAXI J1348-630, Carotenuto et al. 2022).
Regarding MAXI J1535-571, we infer tej = MJD 58017.4+4.0 −3.8 , with a much larger uncertainty with respect to MAXI J1820+070, mostly due to the lack of early-time VLBI observations of the ejecta.Interestingly, our new tej places the ejection in the soft intermediate state (Tao et al. 2018), updating the previous estimation in the hard-intermediate state (Russell et al. 2019).Given the large radio flare reported on MJD 58017.4 by Russell et al. (2019), it is more likely that the ejecta was launched before this date than later, despite our statistical uncertainty being almost symmetrical around the same MJD 58017.4.Interestingly, our preferred ejection date is now roughly 4 days after the well-monitored quenching of the compact jets, which, from the tracking of the evolution of the break frequency from the infrared to the radio bands, appeared to be switched off over a timescale of 1 d on MJD 58013 (Russell et al. 2020).If such a result is confirmed, it would imply that discrete ejecta do not result immediately from the destruction of the compact jets, but that instead they are formed and launched sometime afterwards (see also Echiburú-Trujillo et al. 2024 for MAXI J1820+070).
For MAXI J1535-571, a tentative detection of Type-B QPOs with NICER was reported in Stevens et al. (2018).Specifically, possible Type-B QPOs were detected when stacking the NICER data in the range between MJD 58016.8 and 58025, but it is worth remarking that the authors could not clearly differentiate between Type-A and Type-B QPOs in the data (Stevens et al. 2018).The reported QPO interval overlaps with our inferred ejection date tej = MJD 58017.4+4.0 −3.8 , but, given the 4-d uncertainty, we cannot precisely conclude on the exact sequence of events.However, if Type-B (or Type-A) QPOs were present during the whole stacked interval and lasted until MJD 58025, then we could at least conclude that they persist after the ejecta are produced.
Unfortunately, we are unable to provide constraints on the ejection date for XTE J1752-223, given that tej is a fixed parameter in our modelling.This is also due to the fact that, considering the four available detections, the jets were already strongly decelerating (Yang et al. 2010;Miller-Jones et al. 2011).The radio flare peaking peaking on MJD 55218 and reported in Brocksopp et al. (2013) suggests that the ejection might have happened on that day or before, while it is unlikely to have happened at later times.Our choice of tej on MJD 55217 (one of the possible options) places the ejection in the hard-intermediate state.Type-B QPOs were instead reported on MJD 55218 and 55220 (Shaposhnikov et al. 2010), possibly again suggesting that these variability features could persist after the launch of discrete ejecta.However, the data on XTE J1752-223 do not allow us to draw any firm conclusion on this particular aspect of the jet production.
A summary of the time delay ∆tej,X between the jet ejection date and the first appearance of Type-B QPOs for the sources considered in this work is reported in Table 2.We mention that, rather than the appearance/disappearance of QPOs, it has been recently proposed that jet ejections in GRS 1915+105 could be linked to a change in the coronal geometry observed through changes in the Type-C QPO frequency and a change of sign in the phase lags at the QPO frequency, along with the simultaneous radio emission (Méndez et al. 2022).This result (see also García et al. 2022) appears to be also consistent with previous findings on the non-trivial geometry of the corona in MAXI J1820+070 (Kara et al. 2019) and MAXI J1348-630 (García et al. 2021).

Jet kinetic energy and external ISM density
We discuss in this section our results for the jet kinetic energy.Due to the degeneracy between E0, ϕ and nISM (see Section 3), only the effective energy Ẽ0 = E0/nISMϕ 2 can be independently constrained by our fit.This parameter can be easily interpreted as the kinetic energy required for a ϕ = 1°jet to propagate to the observed distance in a 1 cm −3 ISM.Considering our sample, we constrain an effective energy of Ẽ0 = 2.6 +0.4 −0.4 × 10 46 erg for MAXI J1820+070, a higher value Ẽ0 = 5.8 +16.6  −4.0 × 10 48 erg for MAXI J1535-571 and a lower value Ẽ0 = 1.1 +1.2 −0.6 × 10 45 erg for XTE J1752-223.The highest value obtained for MAXI J1535-571 can be explained by the fact that, assuming the same nISM, the jets from MAXI J1535-571 propagate up to a larger angular distance with less evident deceleration when compared to the other two sources.
Through the constraints on Ẽ0, we plot in the first three panels of Figure 4 the explicit dependence of E0 on nISM for different values of the jet opening angle.Given that no independent information is available regarding the exact values of nISM and ϕ for any of our sources, it is not possible to provide a preferred estimate of the jet kinetic energy.Regarding the jet half-opening angle, the value ϕ = 1°i s generally adopted in the literature (Steiner & McClintock 2012;Steiner et al. 2012;Carotenuto et al. 2022).Such value is consistent with the large number of observational upper limits available for these jets (e.g.Kaaret et al. 2003;Miller-Jones et al. 2006;Russell et al. 2019;Carotenuto et al. 2021;Wood et al. 2021;Williams et al. 2022).Moreover, for cases in which these jets have been resolved in radio or X-rays, the inferred opening angles were of the same order of magnitude (e.g.Bright et al. 2020;Espinasse et al. 2020;Chauhan et al. 2021).Therefore, we assume ϕ = 1°as a reasonable choice, but we also show our solutions for narrower or less-likely wider jets in Figure 4.
On the other hand, much less information is available on the density of the environment of BH XRBs.Two other BH XRBs that displayed large-scale decelerating jets, XTE J1550-564 and MAXI J1348-630, were inferred to be located in a low-density ISM cavity, for which an reasonable external ISM density of 1 cm −3 was assumed (Steiner & McClintock 2012;Carotenuto et al. 2022).Applying the same model to those jets allows us to place independent constraints on the density jump between the interior and exterior of the cavity, but in these cases the value of the internal density still relies on the aforementioned assumption on the external ISM density.In the following, we discuss how our Ẽ0 solutions, coupled to independent constraints on the jet energy, allow us to place some very valuable and informative constraints on nISM for jets decelerating in an uniform ISM.Before that, we note that the jets launched from H1743-322 during its 2003 outburst are considered to be decelerating in a uniform ISM, and Steiner et al. (2012) report an effective energy Ẽ0 = 1.0 +2.2 −0.7 × 10 47 erg.Therefore, we choose to include H1743-322 in our sample for the following considerations, and the dependence of E0 on nISM is shown on the fourth panel of Figure 4.
While the kinetic energy of these jets cannot be directly measured, it is possible to estimate an independent upper limit on the total energy that the system can provide to the jets during the ejection.In fact, under several assumptions, it is fairly straightforward to estimate an upper limit on the jet energy by considering the power available from accretion during the timescale associated to the ejection, and assuming that no energy is transferred to the jets after the launch.If the jet is accelerated during a timescale ∆t, we can write for a single component: where Pjet is the jet power in the rest frame of the source, and can be expressed as the fraction ηjet of the accretion power Ṁ c 2 , which can in turn be traced by the simultaneous X-ray bolometric luminosity LX corrected by the radiative efficiency of the accretion flow η rad (as LX = η rad Ṁ c 2 ).As generally done in the literature, we assume that the duration of the ejection phase ∆t should be roughly equivalent to the duration ∆t obs of the rising phase of the radio flares observed at the moment of launch, which is observable through a radio monitoring with adequate cadence.We note that ∆t can be shorter than ∆t obs if the rise is due to synchrotron self-absorption.Furthermore, in Equation 9we adopt the standard assumption η rad = 0.1 (e.g.Frank et al. 2002;Coriat et al. 2012) and, since we are interested in an upper limit on Ejet, we consider the (roughly) maximum possible jet power by simply assuming ηjet = 1.We note that ηjet could also exceed 1 in presence of a magnetically arrested accretion disk (MAD) and with a contribution from the BH spin (e.g.Bisnovatyi-Kogan & Ruzmaikin 1974;Narayan et al. 2003;Tchekhovskoy et al. 2011;McKinney et al. 2012;Davis & Tchekhovskoy 2020).
After obtaining an upper limit on the jet energy Emax, it is possible to combine such information with the constraints on Ẽ0 to place an upper limit on nISM by: Following this approach, we now discuss the constraints on the external ISM density for each of the sources considered in our sample.

MAXI J1820+070
Starting with MAXI J1820+070, we consider the radio flare produced by the ejection and observed on MJD 58306 with AMI-LA, which was characterized by a rising timescale of 6.7 h and a peak flux density of ≃ 50 mJy at 15.5 GHz (Bright et al. 2020;Homan et al. 2020).We also consider the simultaneous X-ray luminosity from Fabian et al. (2020), extrapolated using the measured spectral parameters to the 0.5-200 keV energy range using the multicomponent feature of the webpimms1 tool.Using Equation 9, we The dependence is obtained through the constraints on the effective energy Ẽ0 and it is here shown for the four sources considered in Section 5.4.The horizontal dot-dashed orange line represents the maximum energy Emax available to the jet from the simultaneous accretion power, which sets a strong upper limit on n ISM for all sources except XTE J1752-223.The dotted black line shows instead the minimum energy E flare in the jet frame inferred from the radio flare associated to the ejection, here implying that such value is likely a large underestimation of the jet kinetic energy, given that it would require extremely low value of n ISM for the jet to propagate up to the observed distances.Regions excluded by our constraints on Emax and E flare are shaded in grey.
Table 3. Source and jet parameters used for the calculations discussed in Section 5.4.From left to right, we list the name of the source, the peak bolometric X-ray luminosity L X used to compute the available accretion power, the derived maximum kinetic energy Emax available to the jet, the duration ∆t obs and the peak flux density S ν,peak at the frequency ν of the radio flare used to compute the jet-frame minimum energy E flare from equipartition, and the inferred upper limits on the external n ISM and jet mass M 0 .Given the multiple sources of uncertainty in the computations of these numbers, we associate a conservative 30% uncertainty to L X and Emax, and a 50% uncertainty to E flare .The references from which these data are obtained are reported in Section 5.4.obtain Emax ≃ 2 × 10 43 erg, which is shown as a horizontal dashdotted orange line in Figure 4, and reported in Table 3. Notably, such upper limit is broadly consistent with the jet internal energy Eint measured by Bright et al. (2020), which obtained values between 10 41 and 10 43 erg, thanks to a reliable estimation of the size of the jet emitting region 90 days after the jet launch.In addition, Espinasse et al. (2020) obtained Eint ≃ 5 × 10 41 erg by resolving the ejecta in the X-ray band with Chandra and measuring the broadband radioto-X-ray spectrum.Both these results were obtained with minimum energy calculations, assuming equipartition between electrons and magnetic fields in the jet plasma (Longair 2011).

Source
Using Equation 10, we can place the following constraint on the ISM density surrounding MAXI J1820+070, assuming ϕ = 1°: nISM,J1820 ≲ 10 −3 cm −3 , (11) with the upper limit being relaxed in case the jet is significantly narrower than 1°(but unlikely given that the jet appeared to be resolved in one of the X-ray detections reported in Espinasse et al. 2020).We also note that Tetarenko et al. ( 2021) measured ϕ = 0.45°+ 0.13°− 0.11°f or the compact jets in the same source.Interestingly, such result on the nISM appears to place MAXI J1820+070 in a low-density region filled with hot/coronal-phase ISM (Cox 2005), similarly to MAXI J1348-630 and XTE J1550-564 (Steiner & McClintock 2012;Carotenuto et al. 2022), but in this case the proper motion data can be adequately described with the propagation in a uniform, low-density ISM.If MAXI J1820+070 is in a low-density cavity, it is possible that these jets were not sufficiently energetic to travel up to the cavity "wall", if present, or, alternatively, the ISM around MAXI J1820+070 might have a much smoother distribution compared to the two sources mentioned before.Finally, we can also obtain a lower limit on the jet kinetic energy, which is represented by the internal energy inferred from the radio flare observed at the moment of ejection, when applying minimum energy calculations (e.g.Fender 2006).The internal energy E flare,obs in the observer frame is computed considering that the flare spectrum evolves from optically thick to optically thin, under the assumption that such evolution is due to decreasing optical depth to synchrotron self-absorption (Fender & Bright 2019).As in most of our cases, if the flare has been monitored at a single frequency ν, we can write: where S ν,peak is the peak flux in mJy and D is the source distance in kpc.We note that this approach has the advantage of being independent of the flare rising timescale, and it can be applied to any flare for which there is evidence for self-absorption and a measurement of the peak flux is available.Most of the flares from BH XRBs detected so far have shown this particular evolution (e.g.Tetarenko et al. 2017;Russell et al. 2019;Carotenuto et al. 2021;Fender et al. 2023).Adopting the formalism from (Fender & Bright 2019), we compute the same energy in the jet frame and we account for the relativistic bulk motion by where δ is the Doppler factor for the approaching jet and we rely on the assumption that α = 0 at the time of peak flux (Fender & Bright 2019).We further note that this lower limit on E0 gives the lowest jet mass M0 = E0/(Γ0 − 1)c 2 and corresponds to the jet composition of pure e ± pairs.For MAXI J1820+070, using Γ0, θ from the fit and Snu, ν, ∆t obs from the AMI-LA flare, we obtain, through Equation 13, E flare,RF ≃ 9 × 10 37 erg, which is also shown with a black dotted line in Figure 4 and reported in Table 3.Given the significant uncertainties associated to this method, we associate a conservative 50% error to these estimations.We note that, from Figure 4, a jet with such kinetic energy would necessarily be propagating in an extremely low-density environment, with nISM ≪ 10 −4 cm −3 .It is reasonable to suggest that the E flare and Emax lines enclose the physical parameter space for the ejecta, and the true (E0, nISM) value should lie between the two horizontal lines.

MAXI J1535-571
Similarly to MAXI J1820+070, we consider the radio flare associated to the jet ejection to estimate the maximum and minimum energies available to the jet.An upper limit on the jet kinetic energy can be obtained with Equation 9, by considering the radio flare produced by the ejection of the S2 component (Russell et al. 2019).Such flare, observed with MeerKAT and ATCA on MJD 58017, had an approximate rising timescale of ≃ 24 h and a peak flux density of ≃ 600 mJy at 1.3 GHz.At the same time, we consider the simultaneous X-ray luminosity from Tao et al. (2018), extrapolated using the measured count rates and reported spectral parameters to the 0.5-200 keV energy range using webpimms.We infer a maximum jet energy of Emax ≃ 5 × 10 44 erg, which, assuming ϕ = 1°, leads to a more stringent constraint on the ISM density surrounding MAXI J1535-571: nISM,J1535 ≲ 10 −4 cm −3 .( Such upper limit is roughly one order of magnitude lower than what obtained for MAXI J1820+070, and also lower than what inferred for MAXI J1348-630 and XTE J1550-564 (Steiner & McClintock 2012;Carotenuto et al. 2022;Zdziarski et al. 2023).This is consistent with the picture of the jets from from MAXI J1535-571 propagating up to a larger distance without the abrupt deceleration observed in the two latter sources.Considering again the radio flare, we can also use the rising timescale, peak flux density and observing frequency to compute with a minimum energy of E flare,RF ≃ 10 39 erg in the jet frame, with the same procedure outlined in the previous section and using Equation 13.Therefore, the true kinetic energy of the jets from MAXI J1535-571 likely lies in the interval between 10 40 and 10 44 erg, which, for instance, is in line with what estimated for the ejecta from by GRS 1915+105 (Mirabel & Rodríguez 1994;Fender et al. 1999;Zdziarski 2014).All the constraints are shown in Figure 4 and reported in Table 3.

XTE J1752-223
We apply the same method used above to the jets launched by XTE J1752-223.First, we consider the radio flare observed with ATCA on MJD 55217, characterized by a rising timescale of ≃ 24 h and a peak flux density of 20 mJy at 5.5 GHz (Brocksopp et al. 2013).
Including the simultaneous X-ray bolometric luminosity reported in the same paper and extrapolated to the 0.5-200 keV range, we use Equation 9 to compute again the maximum energy available to the jet Emax ≃ 5 × 10 44 erg, represented in the third panel of Figure 4.The application of Equation 10 yields the following upper limit on the external ISM density, when assuming ϕ = 1°: which is below the standard galactic nISM ≃ 1 cm −3 , but does not stringently constrain the environmental density as in the cases of MAXI J1820+070 and MAXI J1535-571.On the other hand, given the short distance up to which the jets from XTE J1752-223 propagate (see Figure 3), this result would be consistent with a scenario in which the ejecta have energies similar to the other sources, but travel in a much denser environment.
Considering again the radio flare, we combine the rising timescale, peak flux density and observing frequency to estimate a minimum energy of E flare,RF ≃ 10 37 erg with the same method as before (and Equation 13).Since we do not have a preferred value for the initial Lorentz factor, we assume Γ0 = 3.5, consistent with the lower limit from the fit.Again, we associate a conservative 30% uncertainty to this estimation, and we report all values in Table 3.In particular, we mention that the main source of uncertainty is the duration of the ejection ∆t, which could be overestimated in case of a sparse radio monitoring.Nevertheless, the minimum and maximum energies together with the low effective energy appear to point to a denser environment with respect to the other sources.We note that XTE J1752-223 displayed a complex flaring activity during the 2010 outburst, with rapid oscillations between the intermediate and the soft state, and with the likely production of additional, undetected ejecta  9and from the simultaneous bolometric X-ray luminosity, and it is here shown for the four sources considered in Section 5.4.For ISM densities above the upper limits presented in Section 5.4, marked with dotted vertical black lines, the higher energies required would imply the full accretion power to be supplied to the jets over timescales not compatible with the observed flares and state transitions.Regions excluded by our constraints on ∆t and n ISM are shaded in grey.(Brocksopp et al. 2013).Hence, it might also be possible that jets from this source are intrinsically less energetic than the single pair of ejecta displayed by MAXI J1348-630 or XTE J1550-564 (with the counter-example of GRS 1915+105, Fender & Pooley 2000).Such behavior has been also observed in a significant number of other sources, which are inferred to produce multiple subsequent ejecta, which are, however, rarely detected and spatially resolved (e.g.Homan et al. 2001;Brocksopp et al. 2001;Brocksopp et al. 2002;Fender et al. 2009;Tetarenko et al. 2017;Carotenuto et al. 2021;Fender et al. 2023).As of today, it is unclear if these multiple ejecta that should result from the complex flaring activity are not detected due their lower energy content, or for different reasons, possibly linked to the source environment that might be affected by the previous ejecta.

H1743-322
Lastly, we consider H1743-322 and the radio flare associated to the ejecta that was observed to peak on MJD 52768 with the VLA (McClintock et al. 2009).The flare had a rising timescale of ≃ 24 h, with a peak flux density of ≃ 35 mJy at 4.8 GHz.Combining this information with the extrapolated X-ray bolometric luminosity yields a maximum energy available to the jet Emax ≃ 1 × 10 44 erg (Equation 9), represented in the fourth panel of Figure 4, which then translates through Equation 10 to the following upper limit on the ISM density: nISM,H1743 ≲ 10 −3 cm −3 . (16) Assuming a half-opening angle of 1°, we confirm the requirement for a low density environment, as already pointed out by Hao & Zhang (2009) and Steiner et al. (2012).Considering the radio flare, and again assuming Γ0 = 3 and an inclination angle θ = 75° ( Steiner et al. 2012), the jet in its rest frame has a minimum energy of E flare,RF ≃ 7 × 10 39 erg from Equation 13.

Ejection duration
As an alternative way to represent the constraints obtained on the density of the ISM surrounding our sources, we can use Equation 9to compute the ejection duration associated to the measured LX for any possible value of E0, which is then linked to a specific value of nISM through the definition of Ẽ0.The results are shown in Figure 5 for different values of ϕ.Here, the orange dash-dotted line represents the ejection duration used in the previous subsections to estimate the maximum and minimum energies available to the ejecta.Considering our constraints on Ẽ0, we show that, for MAXI J1820+070, MAXI J1535-571 and H1743-322, if the jet were to propagate in a denser environment with respect to the upper limits reported in Table 3, the ejection would have required a sustained supply of accretion power Ṁ c 2 over timescales much larger than the ones associated with the radio flares, by orders of magnitude.For instance, if the ejecta from MAXI J1535-571 were propagating in a 10 −2 cm −3 ISM with ϕ = 1°, to reach the same distance the system should have supplied energy to the jet at the measured rate for more than 1000 h, which is longer than the time interval between the beginning of the outburst and the inferred jet ejection.
It is crucial to remark that these arguments strongly rely on the assumption that the rising timescale of the radio flares is a good proxy for the true ejection duration, to which we have no direct observational access.We note that this is not true for models that consider continuous jets instead of discrete ejections.In a continuous jet model, the resolved radio knots that we observe are believed to be caused by internal shocks between plasma shells accelerated at different speeds (Kaiser et al. 2000;Jamil et al. 2010;Malzac 2013).While this does not decrease the total amount of energy required for the jet, it relaxes the requirement on the ejection timescale, since the majority of energy is stored in the material of the unseen continuous jet (Kaiser et al. 2000).

BH XRBs in low-density environments
The propagation of the ejecta considered in this work can be adequately described with a single deceleration phase in an homogeneous environment, with an assumed constant ISM density.For three out of the four sources considered in the previous section, combining the upper limits on the maximum available energy with the strong constraints on the effective energy provides us with robust upper limits on the external ISM density, implying that MAXI J1820+070, MAXI J1535-571 and H1743-322 are all harbored in a low-density region in the ISM.For XTE J1752-223, the results are instead not conclusive.At the same time, the motion and the light curve evolution of the ejecta observed from MAXI J1348-630 and XTE J1550-564 also suggest the presence of low-density ISM cavities in which these sources might be embedded, with internal nISM = 10 −3 ∼ 10 −2 cm −3 (Hao & Zhang 2009;Steiner & McClintock 2012), and with either a sharp border (Carotenuto et al. 2022), or a more physical, smooth transition layer (Zdziarski et al. 2023).Despite this difference, it is remarkable to note that for all of the sources displaying large-scale decelerating jets it has been necessary to invoke an environment with a density up to 4 orders of magnitude lower than the canonical ISM density of 1 cm −3 .As already argued in Heinz (2002) for GRS 1915+105 and GRO J1655-40, a low-density environment seems to be a necessary requirement for the jet to propagate up to such large distances (fractions of pc) far from the central compact object.This highlights the importance and the great potential of using the current and future observations of large-scale jets as probes of the environment surrounding BH XRBs.
Despite the emerging scenario, it is currently unclear how such low-density environments might be produced.BH XRBs might be preferentially located in regions occupied by the hot ISM phase (e.g.Ferrière 2001), or, more likely, the low-density region/cavity could formed from the feedback of the system itself.The BH surroundings might have been evacuated by the supernova explosion that created the compact object or by a different type of outflow.Several possibilities (with different degrees of plausibility, see discussion in Hao & Zhang 2009) include the winds from the progenitor star (e.g.Gaensler et al. 2005), winds from the companion star (e.g.Sell et al. 2015) and winds from accretion disk (e.g.Miller et al. 2006;Fuchs et al. 2006;Muñoz-Darias et al. 2016), in addition to the previous activity of the jet itself, whether collimated (e.g.Gallo et al. 2005;Russell et al. 2007;Heinz et al. 2007Heinz et al. , 2008;;Yoon et al. 2011;Coriat et al. 2019) or, as more recently proposed, uncollimated (Sikora & Zdziarski 2023).Interestingly, such scenario appears to be sup-ported by laboratory experiments testing multiple supersonic plasma ejections in a drift chamber (Kalashnikov et al. 2021).After the passage of each ejection, a low-density region can be observed, called "vacuum trace", which causes the subsequent ejections to encounter much less environmental resistance in their propagation.
It is worth mentioning that the properties, distribution and chemistry of the ISM around these systems can be investigated through independent observation at different wavelengths, of which a prime example is the mapping of the molecular line emission from the material shocked by the jets, as already done for GRS 1915+105, GRS 1758-258 and 1E 1740.7-2942(Mirabel et al. 1998;Chaty et al. 2001;Tetarenko et al. 2018Tetarenko et al. , 2020)).Alternatively, it is possible to search for optical Hα emission from the same shocked ISM and independently infer its density from the measurement of the integrated luminosity in the diagnostic line and from realistic assumptions on the shock velocity (Dopita & Sutherland 1996;Russell et al. 2007).Therefore, the importance of these approaches resides also in their potential for partially solving the E0, ϕ, nISM degeneracy that is currently present in our model.

Jet mass
The mass of the ejecta is a parameter of great importance in the study of these systems, and it has not yet been constrained with sufficient accuracy for any source.A mass measurement can directly give us information on the long sought-after composition of the jets from BH XRBs, which is still an open problem.While we have evidence for baryons in the jets from SS 433 from the detection of Dopplershifted iron emission lines in the X-rays (Kotani et al. 1996;Migliari et al. 2002), no information is available from the synchrotron spectra of the other jets, both compact and discrete (Fender 2006), although there have been attempts to model compact jets SEDs with hadronicleptonic models (e.g.Romero et al. 2005;Pepe et al. 2015;Romero et al. 2017;Kantzas et al. 2021).In this context, finding evidence for massive ejecta could strongly suggest the presence of cold protons, balanced by a long tail of non-relativistic electrons (Carotenuto et al. 2022).This is also supported by Zdziarski et al. (2023), which argue that a massive jet is unlikely to be pair dominated, and by more recent results reported in Zdziarski & Heinz (2024), which suggest the existence of a fundamental difference in composition between compact jets, pair-dominated, and discrete ejecta, which should instead have a baryonic composition.Evidence of protons in these jets bolster arguments that BH XRBs could represent a class of PeV cosmic ray sources (e.g.Fender et al. 2005;Cooper et al. 2020).However, the mechanism for which protons are loaded in the jets remains unclear, as they could be already present at the jet formation or could be entrained during the early phases of the jet motion (see for instance O' Riordan et al. 2018 andKantzas et al. 2023).
In this work, we cannot directly provide estimations of the masses of the ejecta for the four sources considered, but we can still place interesting upper limits on it by writing and using the best value of Emax and Γ0 from the sections above.
With the values reported in Table 1, we find the ejecta launched by MAXI J1820+070 to have a mass which is equivalent to an upper limit of ≲ 8 × 10 −12 M⊙.This estimation is consistent with the ≃ 10 20 g obtained with minimum energy calculations, assuming one proton per electron (Espinasse  et al. 2020), and it is also consistent with jet masses estimated in other sources with the same method (e.g.Fender et al. 1999;Gallo et al. 2004).
Regarding MAXI J1535-571, we infer a higher upper limit: or ≲ 5 × 10 −10 M⊙.It is likely that the ejecta from MAXI J1535-571 have a larger density contrast than MAXI J1820+070 with the surrounding ISM, as this is probably required in order to propagate a larger distances.Notably, the ejecta from MAXI J1535-571 were unresolved in all the ATCA and MeerKAT detections (Russell et al. 2019), hence it is not expected to have a volume significantly larger than the jet from MAXI J1820+070.Assuming a most likely value Γ0 = 3.5, the mass of the jets from XTE J1752-223 is inferred to be equivalent to an upper limit of ≲ 1×10 −10 M⊙.In this case, while the mass upper limit lies in between what obtained for the two previous sources, the early jet deceleration observed in XTE J1752-223 likely points towards a denser environment (see Section 5.4).Similarly to XTE J1752-223, by assuming Γ0 = 3 in the case H1743-322, we can constrain equivalent to an upper limit of ≲ 3 × 10 −11 M⊙.All the inferred upper limits on M0 are reported in Table 3.Assuming a high accretion rate of ∼10 18 g s −1 (roughly 0.1L Edd ), typical of the hard-to-soft state transition (Maccarone 2003), and assuming that during the ejection the majority of the accreted mass is channeled into the jets, it would take from minutes to hours to accumulate a jet mass in the range 10 20 ∼ 10 22 g.On the other hand, it is generally known that the mass outflow rate for thermal winds can be up to ten times higher than the simultaneous accretion rate (e.g.Higginbottom & Proga 2015;Dubus et al. 2019).This appears to be consistent with the general picture in which, for BH XRBs, most of the outflow mass is carried by winds, while most of the kinetic feedback is carried by jets (Fender & Muñoz-Darias 2016).

The sample of large-scale jets
After discussing the results of the modeling work presented in this paper, we can consider for the first time the entire sample of large-scale jets detected so far and for which we have information on the source parameters, with the aim of looking for possible interesting and informative trends/correlations.The current sample includes the three sources considered in this work, namely MAXI J1820+070, MAXI J1535-571 and XTE J1752-223, with the addition of MAXI J1348-630 (Carotenuto et al. 2021(Carotenuto et al. , 2022)), XTE J1550-564 (Sobczak et al. 2000;Hannikainen et al. 2001;Wu et al. 2002;Corbel et al. 2002;Steiner & McClintock 2012) and H1743-322 (McClintock et al. 2009;Steiner et al. 2012).We further add GX 339-4, that displayed large-scale jets in 2003 and for which, albeit not detecting deceleration, we have a lower limit on the initial Lorentz factor Γ0 > 2.3 (Gallo et al. 2004).
We first compare the initial jet speed (Γ0) with the measured dimensionless BH spin a * , as can be seen from panel (a) of Figure 6.At visual inspection, it is not clear whether there is any evidence of correlation between the two parameters, even if we note that we only have three estimations of Γ0 in our sample, while the rest are lower limits.Furthermore, we note that these spin values are obtained with different methods that often yield different results (e.g.Reynolds 2021; Draghis et al. 2023).Specifically, the spin of H1743-322 is obtained through continuum fitting (Steiner et al. 2012), while the spins of MAXI J1348-630, XTE J1752-223, MAXI J1535-571 and GX 339-4 have been obtained with relativistic reflection modeling (Parker et al. 2016;García et al. 2018;Miller et al. 2018;Jia et al. 2022).Lastly, the spins of XTE J1550-564 and MAXI J1820+070 were obtained with the application of the RPM (Relativistic Precession Model, e.g.Stella & Vietri 1998, 1999;Stella et al. 1999), from Motta et al. (2014) and Bhargava et al. (2021), respectively.From our sample, the BH spin does not seem to have a significant effect on the initial jet speed.
With the aim of comparing the jet speed with the simultaneous accretion rate, we show in panel (b) Γ0 as a function of the simultaneous bolometric X-ray luminosity LX obtained from the literature and converted to the 0.5 − 200 keV energy range with webpimms, as in Section 4.2.We plot LX in units of L Edd , which is the Eddington luminosity L Edd ≃ 1.3 × 10 39 (MBH/10 M⊙) erg s −1 that represents the limit for spherically stable hydrogen accretion (Frank et al. 2002).To estimate L Edd we need a measurement of MBH, but we note that a dynamically confirmed BH mass is only available for XTE J1550-564 and MAXI J1820+070 (Orosz et al. 2011;Torres et al. 2020).In this plot we updated the mass of MAXI J1348-630, which was previously estimated from the normalization of the disk blackbody component in the X-ray spectral fitting reported in Tominaga et al. ( 2020), for a non-spinning BH.Adopting the recent spin a * = 0.78±0.02measurement by Jia et al. (2022), we computed the spin-dependent radius of Innermost Stable Circular Orbit (ISCO) and then use this updated parameter to obtain MBH = 15.2 ± 2.3 M⊙ (see Tominaga et al. 2020 for the explicit dependence of the mass on the ISCO radius).Moreover, we included the new BH mass estimation MBH = 12 ± 2 M⊙ obtained for H1743-322 through X-ray reflection spectroscopy (Nathan et al., submitted) and the new mass estimation of MBH = 12 ± 1 M⊙ obtained for XTE J1752-223 through the analysis of the the soft state and the soft-to-hard spectral state transition (Abdulghani et al. 2024).As for the previous panel, we do not have enough robust estimations of Γ0 and MBH to draw a conclusion, but we note that this plot, after Fender et al. (2004), is starting be populated and will be highly relevant once more measurements become available.
Another interesting comparison, that we show in panel (c) of Figure 6, can be done between the jet speed and the internal energy inferred from the radio flare observed at the moment of ejection, when applying minimum energy calculations (Fender & Bright 2019, Equation 12) and converting the internal energy to the jet frame, accounting for the bulk relativistic motion with Equation 13 (with 50% uncertainties associated).We notice that most flares in our sample had a minimum energy ranging between 10 37 and 10 40 erg.For H1743-322 and XTE J1550-564, we assumed Γ0 = 3 and inclination angles of, respectively, θ = 75° (Steiner et al. 2012) and θ = 70° (Steiner & McClintock 2012), and for which we obtain E flare,RF ∼ 10 41 erg.We caution that the energies reported in this plot are affected by large uncertainties, both in the peak flux (the peak could be missed in the monitoring or it could be optically thin) of the emitting region and in the conversion from the observer frame to the jet frame.Again, a visual inspection seems to show that there is no clear correlation between these two parameters.In the current understanding, Γ0 is the bulk Lorentz factor of the whole jet, while E flare,RF should only result from the relativistic electrons present in the jet plasma.If discrete ejecta have a predominant baryonic composition (e.g.Zdziarski & Heinz 2024), the mass of the proton will be more important in determining Γ0 than the energy contained in the relativistic electrons.
Lastly, we compare the de-projected distance traveled by the jet with Γ0, as shown in panel (d) of of Figure 6.Due to the uncertainties on the source distance and jet inclination angle, we assume a conservative 30% uncertainty on the physical distance traveled by the jet.The initial jet speed does not appear to be a driving factor in determining the distance at which the jet propagates to.We might in fact expect that the distance could be more correlated with the jet mass (or with the jet/ISM density contrast, Savard et al. in prep.)than with its speed.In this context, a massive jet with a low Lorentz factor will propagate further in a given ISM density than a lighter jet with a higher Γ0 and the same kinetic energy.

CONCLUSIONS
In this paper, we have presented a physical modelling of the motion of the decelerating jets launched by MAXI J1820+070, MAXI J1535-571 and XTE J1752-223.Adopting a Bayesian approach, we fitted the jet angular distance data with the dynamical blast-wave model developed by Wang et al. (2003), and we found that the model provides an excellent description of the jet motion, from the first phase of ballistic motion to the final deceleration phase.In particular, a single Sedov phase in a homogeneous ISM appears to adequately capture the dynamics of the decelerating jets for the entire sample considered.The results obtained from a simple model derived from GRBs demonstrate the high potential of applying some of the welldeveloped theoretical advancements on GRBs to the jets from BH XRBs.These discrete ejecta can be considered as less-relativistic analogues of GRB jets, with the advantage of providing better access to their physics due to their location in the Galaxy.
From the fits, we are able to place constraints on multiple physical parameters of the jets, including a first estimation of the initial Lorentz factor of the ejecta from MAXI J1820+070 (Γ0 = 2.6 +0.5 −0.4 ) and estimates for the Lorentz factor (Γ0 = 1.6 +0.2 −0.2 ), ejection date (MJD 58017.4+4.0 −3.8 , soft-intermediate state) and inclination angle θ = 30.3°+ 6.3°− 6.3°f or the ejecta produced by MAXI J1535-571.By considering the constraints on the effective energies and on the maximum energy available to the jets from the accretion power, we are also able to provide new upper limits on the jet mass and on the ISM density surrounding our sources.Overall, our results support the emerging scenario for which BH XRBs displaying large-scale jets appear to be mostly located in low-density environments.
Considering the current sample of large-scale jets, we find no clear correlations between the initial jet speed and the BH spin, the simultaneous accretion rate, the jet minimum energy inferred from the flare at the moment of ejection, or the distance traveled by the jet.It is worth mentioning that the lack of evident correlation between most of the parameters considered is nevertheless informative for our current understanding of jets from BH XRBs.While our sample is currently limited to a small number of sources, more observations of decelerating ejecta from BH XRBs are needed in order to test our results, an this will lead us to a significant progress in our understanding of the jet production, propagation, and feedback on the surrounding environment.In this context, two new large-scale decelerating ejecta have been recently discovered in MAXI J1848-015 (Tremou et al. 2021;Bahramian et al. 2023) and4U 1543-47 (Zhang et al. in prep.), and these jets represents ideal targets for the continuation of this work.In the near future, this approach will be greatly improved by the joint modelling of kinematics and radiation from these ejecta (Cooper et al. in prep.) and will also benefit from the comparison with the first relativistic hydrodynamic simulations of these objects (Savard et al. in prep.).
It is worth to mention that the data currently available for MAXI J1820+070 are of outstanding quality, with a monitoring campaign that covered the first phases of jet motion (VLBI) as well as the final deceleration phase, allowing us to obtain much smaller statistical uncertainties on the model parameters with respect to the other sources.In the coming years, similar monitoring campaigns can be planned and performed in order to cover the entirety of the jet evolution.Milliarcsec-resolution VLBI observations are crucial to observe the jets much closer to the compact object and hence to obtain more stringent constraints on their physical parameters, especially the ejection date (Miller-Jones et al. 2019;Wood et al. 2021Wood et al. , 2023)).In the future, the new generation Event Horizon Telescope (ngEHT) will also be able to perform extremely high-resolution (order of ∼10 µas) observations of the jets in the mm range (e.g.Johnson et al. 2023).At the same time, dense and long-term monitoring campaigns with sensitive interferometers such as MeerKAT, which already detected a large number of discrete ejecta, but also the future SKA-MID (Braun et al. 2015) and ngVLA (Selina et al. 2018) will be fundamental to follow the deceleration phase and to probe the final phases of the jet evolution.

APPENDIX A: POSTERIOR DISTRIBUTIONS
We show in this section the corner plots with the posterior distributions of the model parameters for the three sources considered in this work.
This paper has been typeset from a T E X/L A T E X file prepared by the author.
are shown in Figure 3, along with the proper motion of the jet, while the posterior distributions for the model parameters are shown in Figure A3.According to the model, the jet is launched with an effective energy of Ẽ0 = 1.1 +1.2

Figure 1 .
Figure 1.Angular separation in arcsec between the discrete ejecta and the position of MAXI J1820+070, with data from Bright et al. (2020); Espinasse et al.(2020)  andWood et al. (2021).The un-shaded, gray and seashell regions mark periods in which the source was, respectively, in the hard, intermediate and soft state(Shidatsu et al. 2018).The black horizontal dashed line represents the zero separation from the core, while the black continuous line represents the best fit obtained with the external shock model.The orange shaded area represents the total uncertainty on the fit and it is obtained by plotting the jet trajectories corresponding to the final positions of the MCMC walkers in the model parameter space.Residuals ([data -model]/uncertainties) are reported in the bottom panel.The model appears to provide an excellent description of the motion of both the approaching and receding ejecta, with a low statistical uncertainty.

Figure 2 .Figure 3 .
Figure 2. Same as Figure 1, but for MAXI J1535-571, with data from Russell et al. (2019) and information on the spectral states fromTao et al. (2018).Dark and light gray regions differentiate, respectively, between the HIMS and the SIMS.The model appears to fit the data remarkably well, implying that a Sedov phase is an adequate physical scenario for the jet deceleration.

Figure 4 .
Figure 4. Explicit dependence of the jet kinetic energy E 0 on the external ISM density n ISM for different values of the half-opening angle ϕ.The dependence is obtained through the constraints on the effective energy Ẽ0 and it is here shown for the four sources considered in Section 5.4.The horizontal dot-dashed orange line represents the maximum energy Emax available to the jet from the simultaneous accretion power, which sets a strong upper limit on n ISM for all sources except XTE J1752-223.The dotted black line shows instead the minimum energy E flare in the jet frame inferred from the radio flare associated to the ejection, here implying that such value is likely a large underestimation of the jet kinetic energy, given that it would require extremely low value of n ISM for the jet to propagate up to the observed distances.Regions excluded by our constraints on Emax and E flare are shaded in grey.

Figure 5 .
Figure 5. Explicit dependence of the ejection duration ∆t on the external ISM density n ISM for different values of the half-opening angle ϕ.The dependence is obtained through the constraints on the effective energy Ẽ0 , from Equation9and from the simultaneous bolometric X-ray luminosity, and it is here shown for the four sources considered in Section 5.4.For ISM densities above the upper limits presented in Section 5.4, marked with dotted vertical black lines, the higher energies required would imply the full accretion power to be supplied to the jets over timescales not compatible with the observed flares and state transitions.Regions excluded by our constraints on ∆t and n ISM are shaded in grey.

Figure 6 .
Figure 6.Comparison between the inferred initial Lorentz factor Γ 0 of the ejecta in our sample of large-scale jets and: (a) -the dimensionless spin parameter a * , (b) -the bolometric X-ray luminosity L X simultaneous to the ejection in Eddington units, (c) -the jet frame internal energy E flare inferred from the radio flare associated to each ejection (see text for details), (d) the de-projected distance traveled by the jet.We find no clear correlation between Γ 0 and the parameters shown in the four panels, and more sources are needed to increase the sample of large-scale jets.

Figure A1 .
Figure A1.Corner plots showing the constraints on the physical parameters of the ejecta from MAXI J1820+070.The panels on the diagonal show histograms of the one dimensional posterior distributions for the model parameters, including the jet initial Lorentz factor, effective energy, inclination angle and ejection time (here represented as MJD −58300), as well as the source distance.The median value and the equivalent 1σ uncertainty are marked with vertical dashed black lines.The other panels show the 2-parameter correlations, with the best-fit values of the model parameters indicated by green lines/squares.The plot was made with the corner plotting package (Foreman-Mackey 2016).

Figure A3 .
Figure A3.Corner the constraints on the physical parameters of the ejecta from XTE J1752-223, same as Figure A1.Note the log scale for Γ 0 .