© The Authors 2025

1. Introduction

Accretion of matter onto celestial bodies is a universal phenomenon giving rise to the formation of jets and other types of outflows observed across various scales and scenarios, from proto-planetary disks to gamma-ray bursts and from stellar-mass to super-massive black holes (SMBHs). In the case of accreting stellar-mass black holes (BHs), this process reaches its extremes: These objects contrive to feed back to the environment via jets and winds a large fraction of the energy and matter they could potentially have swallowed, thereby actively impacting their environment rather than behaving only as sinks. The significance of this feedback extends across diverse scales, from SMBHs fuelling active galactic nuclei (AGNs) to stellar-mass BHs in X-ray binaries, which are considered small-scale analogues of SMBHs (see e.g. Merloni et al. 2003; Falcke et al. 2004; Motta et al. 2017; Fernández-Ontiveros & Muñoz-Darias 2021). The influence of SMBH activity is believed to play a crucial role in regulating galactic-scale activities, for example, by triggering/quenching star formation and affecting the chemical enrichment of the intergalactic medium as well as the overall evolution of galaxies and large structures (Magorrian et al. 1998; Croton et al. 2006). Similarly, stellar-mass BHs impact their surroundings by reintroducing a fraction of the infalling gas and liberating gravitational energy, thus energising the interstellar medium (ISM) and increasing/decreasing its density and even injecting high energy particles into the ISM (LHAASO Collaboration 2024; Martí-Devesa & Olivera-Nietov 2025). These phenomena induce interstellar turbulence that may stimulate local star formation and may even introduce seed magnetic fields into the ISM, potentially contributing to the average magnetic field observed in the Galaxy (see e.g. Heinz et al. 2008; Mirabel et al. 2015.

Theoretical models describing the interaction between relativistic jets from BHs and their environment have postulated the emergence of a shock. Ejected jet particles are expected to induce the formation of a radio lobe that expands to give rise to a shock-compressed gas bubble populated by relativistic electrons, which should generate non-thermal emissions (Scheuer 1974; Begelman & Cioffi 1989; Kaiser & Alexander 1997; Kaiser et al. 2004). Large extended structures around powerful extragalactic radio sources of type FR II (Fanaroff & Riley 1974) ascribed to this phenomenon have been observed for decades. Although, similar albeit smaller structures are expected around stellar-mass BHs (Kaiser et al. 2004), only a handful of examples exist. A possible explanation is that the formation of such features in the ISM may not be universally observed at each interaction site, this being subject to the properties of both the local environmental conditions (e.g. density and chemical composition) and the intrinsic properties of the jets (e.g. velocity, launch direction, duty cycle, lifetime). Yet, the attributes of jet–ISM interaction regions encapsulate significant insights into poorly understood characteristics of jets, specifically the total jet power, radiative efficiency, speed, and matter content.

Analogous to the studies of the radio lobes associated with AGN jets, the jet–ISM interaction regions can be used as calorimeters of the ‘power × lifetime’ product of the jet. This approach has already been adopted in the case of Cyg X-1. Deep radio observations of this source revealed a semi-circular shock region with a diameter of ∼5 pc aligned with the resolved AU-scale radio jet (Stirling et al. 2001) that develops at the location where the collimated jet impacts the ISM, where a radio lobe is inflated (Gallo et al. 2005). To date, this structure and similar ones previously discovered near the peculiar system SS433 (Begelman et al. 1980; Dubner et al. 1998) are the clearest examples of jet-inflated radio lobes near a BH X-ray binary (or microquasar) in the Galaxy (but see Sell et al. 2015 for an alternative explanation for the case of Cyg X-1). A few more examples have also been found, including around the neutron star X-ray binary Cir X-1 (Sell et al. 2010; Coriat et al. 2019), in the vicinity of the peculiar BH X-ray binary in the nearby galaxy NGC 7793 (Pakull et al. 2010; Soria et al. 2010; Hyde et al. 2017), in the Large Magellanic Cloud (Hyde et al. 2017), and in NGC 300 (Urquhart et al. 2019).

GRS 1915+105 is one of the best-studied Galactic BH X-ray binaries. Located at a radio parallax distance of 9.4 ± 0.7 kpc (Reid & Miller-Jones 2023), GRS 1915+105 hosts a stellar-mass BH with mass ∼11.4 M_⊙ (Reid & Miller-Jones 2023), which – when active – is believed to accrete erratically at a mass accretion rate close to its Eddington limit (Done et al. 2004). This system was the first Galactic BH observed to display relativistic, apparently super-luminal, radio ejections (Mirabel & Rodríguez 1994), and it is still regarded as the archetypal Galactic source of relativistic jets.

GRS 1915+105 appeared as a bright transient in August 1992 (Castro-Tirado et al. 1992) and remained very bright in X-rays and radio until recently (Motta et al. 2021). Over the past 30 years, it has alternated between extraordinary activity phases with strong and variable radio and X-ray emission, and long periods of reduced activity, with relatively low radio and X-ray fluxes, with the latter being associated with a hard state X-ray energy spectrum (Belloni et al. 1997). In such a hard state, the source is thought to be ‘jet-dominated’; that is, the majority of the liberated accretion power is in the form of a radiatively inefficient jet and is not dissipated as X-rays in the accretion flow (Fender et al. 2003). In radio, an extended milli-arcsecond scale core radio jet is commonly observed (Dhawan et al. 2000), which is variable on relatively short timescales (< 1 h). Transitions to softer accretion states are usually accompanied by strong radio flares, which are associated with relativistic ejections that produce extended jets. Such jets have been repeatedly resolved, propagating out at approximately 30° from the north-south direction on the plane of the sky to distances between a few (Dhawan et al. 2000; Rushton et al. 2007) to tens or hundreds of milli-arcseconds (Fender et al. 1999; Miller-Jones et al. 2007).

Searches for extended jet–ISM interaction structures around GRS 1915+105 at radio wavelengths have been unsuccessful (Rodriguez & Mirabel 1998). Only two relatively compact regions with flat radio spectra symmetric about GRS 1915+105 – IRAS 19124+1106 and IRAS 19132+1035 – were identified as candidate interaction sites as early as the late 1990s (Rodriguez & Mirabel 1998). Both of these sources are located approximately 17 arcmin away from GRS 1915+105, at about 40 pc if located at the same distance as GRS 1915+105, and their flat radio spectra are incompatible with synchrotron emission. The H I distance to IRAS 19132+1035 places the hot spot at a distance of 6.6 ± 1.4 kpc from the observer (Chaty et al. 2001), which is marginally consistent with the most recent parallax distance of GRS 1915+105 (9.4 ± 0.7 kpc, Reid & Miller-Jones 2023). Kaiser et al. (2004) developed a self-similar model for the large-scale jets from GRS 1915+105, applying a fluid dynamical model originally developed for extragalactic jets to this microquasar to infer the properties of the jet presumably responsible for the interaction (Kaiser & Alexander 1997). These authors concluded that while a bow shock structure should exist, it had not been observed around GRS 1915+105 (nor in the vicinity of other X-ray binaries) due to the limited sensitivity of the available observations. More recently, using ALMA line observations, Tetarenko et al. (2018) showed that the southern jet from GRS 1915+105 collides with a molecular cloud coincident with IRAS 19132+1035, which is likely heated by a young medium-mass star cluster. These results further supported the connection of IRAS 19132+1035 with the jets from GRS 1915+105, although they indicated that the interaction between the BHXB jet and the ISM in this region occurs on smaller scales than previously hypothesised.

In this paper, we report on the discovery of a large-scale structure near GRS 1915+105 in sensitive MeerKAT observations with properties consistent with those expected if the jet developed a bow shock, as predicted by Kaiser et al. (2004). The paper is structured as follows: Section 2 describes the observations we carried out and the data reduction. Section 3 describes the properties of the structure we identified in the images. Section 4 reports the detailed calorimetry of the jets, which we employed to estimate the energy transported in the jets and transferred to the ISM. Section 5 discusses our results and their implications. Finally, Sect. 6 contains our concluding remarks.

2. Observations and data reduction

We observed GRS 1915+105 for a total of 65 times with the MeerKAT interferometer as part of the ThunderKAT Large Survey Project (Fender et al. 2016), which between 2018 and 2023 routinely observed active X-ray binaries, cataclysmic variables, supernovae, and gamma-ray bursts. We started observing GRS 1915+105 (J2000 19^h15^m11.56^s +10°56′44.9″) with MeerKAT soon after the start of operations, on 2018-12-08, and subsequently every several weeks until March 2020 (Motta et al. 2021). On March 18th, 2020 the AMI-LA interferometer, which has been carrying out a long-term monitoring of GRS 1915+105, interrupted operations due to the Covid-19 outbreak (Motta et al. 2021). Thus, we started weekly monitoring of the target with MeerKAT, which ended in April 2021. After this date, sparse MeerKAT observations of GRS 1915+105 were taken until June 2023. The observations used in this work are listed in Appendix A.1 and are all publicly available on the SARAO archive¹.

We used the telescope’s L-band receivers and obtained data at a central frequency of 1.28 GHz across a 0.86 GHz bandwidth (856–1712 MHz). The correlator was configured to deliver either 4096 or 32 768 channels, with an 8-second integration time per visibility point. In either case, data were binned down to 1024 channels for consistency before any further analysis. Between 58 and 64 of the 64 available dishes were used in the observations, with a maximum baseline of 7.698 km. The first MeerKAT observation (observing block 1544268663) consisted of a 90 min run, of which 60 min was on-source, 20 min on the primary (flux and bandpass) calibrator (J1939−6342), and 3 minutes on the nearby secondary (phase) calibrator (J2011−0644). All the other observations consisted of 15 minutes of on-source time, bookended by two 2-minute scans of the secondary calibrator and a 10-minute observation of a primary calibrator.

We conducted the subsequent analysis via a set of PYTHON scripts specifically tailored for the semi-automatic processing of MeerKAT data (O XKAT², Heywood 2020). The CASA package (CASA Team 2022) was used to flag the first and final 100 channels from the observing band, autocorrelations and zero amplitude visibilities. Further flagging of the data was performed to remove RFI in the time and frequency domain. Flux density scale, bandpass and delay corrections were derived from the scans of the primary calibrator, while the complex gain and delay corrections were obtained from the scans of the secondary calibrator. A spectral model for the phase calibrator was derived starting from the flux and bandpass calibrator, by temporarily binning the data into 8 equal spectral windows. The above corrections were applied within CASA to the target field, which was split to a different measurement set and further flagged using the TRICOLOUR package³, after averaging the data in time (8 s) and frequency (8 channels) for imaging purposes.

We used WSCLEAN (Offringa et al. 2012) to obtain a first image of the entire square-degree MeerKAT field by combining a subset of 10 observations where GRS 1915+105 showed a relatively low flux density (below ∼10 mJy; see Motta et al. 2021). By combining several observations we obtained a better uv-coverage compared to that of individual observations, which is crucial for a good image deconvolution, while by selecting only observations where the target was relatively faint we mitigated the artefacts that may arise in imaging bright and/or variable sources with a sparse uv-coverage. We used the resulting image to generate a deconvolution mask, which we adopted to image each observation individually. Each measurement set was then self-calibrated using CUBICAL (Kenyon et al. 2018) to solve for phase and delay corrections for every 32 seconds of data.

Since GRS 1915+105 showed dramatic flux density variations, spanning a dynamic range of over 2 orders of magnitude (i.e. from ≈1 mJy to over 600 mJy), we uv-subtracted it from each observation before combining several observations to avoid the effects of a strong and variable source near the phase centre⁴ when building a deep image of the field. We first used WSCLEAN to generate a model of the field with higher frequency resolution (16 channels instead of 8) for each observation. Then we partitioned the model into a second FITS cube containing only GRS 1915+105, which we subsequently subtracted from the model visibilities of the full sky. Finally, we used WSCLEAN to combine all the observations, to obtain a deep image of the field without the point source corresponding to GRS 1915+105. We exclude from the final stacked data-set those measurement sets that individually yielded images affected by severe artefacts.

The resulting map has been primary beam corrected using the KATBEAM package⁵, blanking the map beyond the nominal 30% level⁶. The resulting wide-field image is presented in Fig. 1, which constitutes the deepest radio image of this field in L-band at the time of writing, with an overall exposure of 14 hours.

Fig. 1. Field of GRS 1915+105 as seen by the MeerKAT interferometer at 1.28 GHz in a total on-source time of 15.5 hr. The restored beam for this map is circular and has a radius of 6.6″. The circle marks the position of GRS 1915+105, and the two arrows mark the direction of the jets identified in the ’90s (Fender et al. 1999). In the image, north is up and east is left.

3. Results

The MeerKAT image of the GRS 1915+105 field is presented in Fig. 2. The field is complicated with multiple large-scale emission regions and compact sources. GRS 1915+105 is in the middle of the field, represented by the white dot. The direction of the radio jets, which have been repeatedly resolved as extended radio jets on a range of scales from ∼1 milli-arcsecond to hundreds of arcseconds (Dhawan et al. 2000; Miller-Jones et al. 2005, 2007; Rushton et al. 2007; Fender et al. 1999), is denoted in Fig. 1 with white arrows. We discover an arched structure that appears approximately 17 arcmin to the south-east of the GRS 1915+105 position, with an apparent diameter of 10 arcmin, and an average flux density between 0.1 and 0.2 mJy (see Fig. 2). The arch appears to be connected with a large HII region (CHIMPS 4894⁷; see Fig. 2) that extends on either side of the bright extended source IRAS 19132+1035 and aligns with the direction of the extended radio jets from GRS 1915+105. These properties suggest that the structure we observe is compressed material behind a bow shock, which may have formed at the head of a large-scale dark jet (i.e. not directly imaged in the radio band) that rams into dense ISM material, carving out a cavity in the ISM.

Fig. 2. Zoom-in of the jet–ISM interaction region. The bow shock structure is marked by the cyan region, which has been defined by eye. The position of GRS 1915+105 and IRAS 19132+1035, as well as the jet direction, are marked.

The observed configuration is reminiscent of what was observed around Cyg X-1, where a bow shock formed as a large-scale jet connected with the tail of a nearby HII nebula (Sh2-101, or the Tulip nebula, Gallo et al. 2006). The large-scale dark jet in this case was able to blow a cavity into the ISM with a diameter of approximately 5 pc, at a distance of 6 pc from the BH⁸. For the case of Cyg X-1, Sell et al. (2015) stressed that the possibility of a massive O-star wind playing a role in inflating the bubble marked by the bow shock cannot be ruled out. In principle, a similar situation could exist for GRS 1915+105, for which we have no final proof that the structure we observe is really connected with a bow shock due to the jets, and not a fore/background structure. However, we note that the companion star in GRS 1915+105 is not massive, but an evolved low-mass companion (Reid et al. 2014), which is not likely to generate strong winds.

Based on the above, in this work, we postulate that the structure we discover in the MeerKAT data is induced by the jet interacting with the surrounding medium. Under this assumption, we note that the cavity blown by the jet into the ISM is significantly larger than in the case of Cyg X-1. At a distance of ≈9.4 kpc and assuming a source inclination angle i of 60° (Reid et al. 2014), the bubble near GRS 1915+105 has a physical diameter of approximately 30 pc and formed ≈42 pc away from the binary, and the associated arched structure is 16 times more luminous than that near Cyg X-1.

The postulated jet–ISM interaction site near GRS 1915+105 shows a complex structure formed by various areas. In Fig. 3 (panel A), we show a sketch of the region marking its various components, which we also describe below:

• IRAS 19132+1035 (henceforth the IRAS region), which emits flat-spectrum thermal radiation (see Rodriguez & Mirabel 1998). This feature is also characterised by molecular line emission, which indicates that the gas in the region is a mix of hot and cold gas, the latter possibly hosting young massive stars, which could contribute significantly to the emission from the region (Tetarenko et al. 2018).

• A non-thermal feature (henceforth the northern feature) to the north-east edge of IRAS 19132+1035 pointing towards GRS 1915+105. This region is the site of particle acceleration that originates at the end of the large-scale jet. This structure was already tentatively associated with the jet from GRS 1915+105 in the late ’90s (Rodriguez & Mirabel 1998) and Tetarenko et al. (2018) confirmed the association by showing the presence of molecular lines near this region, consistent with the effect of an active jet pushing into cold gas. The northern feature emits steep-spectrum synchrotron radiation (see Rodriguez & Mirabel 1998), and in a AGN would be identified with a hot-spot (Kaiser et al. 2004).

• An arched structure (henceforth the bow shock structure) with an average brightness of 0.15 mJy/beam that extends to either side of the IRAS region. This feature, which was not known before and has been discovered by MeerKAT, likely marks the presence of shock-compressed material and of a bow shock. In the old VLA data (Rodriguez & Mirabel 1998) the structure was not detected due to limited sensitivity, although Kaiser et al. (2004) predicted it, and estimated that the brightness of the feature should have been of the order 0.1 mJy/beam. This brightness level is well below the sensitivity of the VLA observations, but remarkably consistent with what we observed with MeerKAT. The superior surface brightness sensitivity of MeerKAT is also likely responsible for the detection of the dim diffuse emission that other instruments would resolve out.

Fig. 3. Panel A: schematic picturing the jet–ISM interaction region as observed in the MeerKAT image. Panel B: sketch of the model employed to interpret the observed structure.

4. Analysis

To interpret the interactions of the large-scale jet of GRS 1915+105 with its environment, we employed the model originally developed by Kaiser & Alexander (1997) for the jets of radio galaxies and later adapted to Galactic sources by Kaiser et al. (2004). We assumed that the supersonic jets emerging from a central source end in strong shock fronts at the location where they impact the ambient medium. The plasma from the jets will inflate two lobes that are overpressured with respect to the environment, and the pressure in the lobes should be high enough that the jets are confined while propagating through this region so that the energy transported by the jet is deposited at the jet end shock without significant losses. The lobes, therefore, expand forwards and sideways, and thus a bow shock forms, while the jets – expanding in a low-density but high-pressure medium – are protected against disruption due to the turbulence in the gas in the environment. The jet direction is assumed to remain constant over time, as the model requires that the jets propagate in an environment that has been cleared out by previous ejections.

The sketch in Fig. 3 (panel B) shows the model we employ to interpret the structure we observed in the MeerKAT image. The site of particle acceleration at the end of the jet, the excited multi-phase gas which the jet interacts with, and the shocked compressed material are visible in the radio band as the northern feature, the IRAS region, and the bow shock structure, respectively. The large-scale jet from GRS 1915+105 and the lobe fed by the jet are not visible, possibly due to sensitivity reasons. This morphology could be seen as a small-scale analogue of what is seen around, for example, Cygnus A (Carilli et al. 1991).

With this model in mind, we could determine the properties of the jet. We proceeded as follows:

1. We first derived the density of the shocked material at the end of the jet by assuming that the jet hits an overdense region, corresponding to the IRAS region, and by assuming that the continuum emission is predominantly Bremsstrahlung. Then, we derived the pre-shock ISM density.

2. We modelled the non-thermal northern feature assuming that the continuum emission is predominantly due to synchrotron and derived the minimum energy magnetic field and pressure.

3. We modelled the bow shock structure assuming that the emission is due to synchrotron in analogy to what is done for the northern feature, and we derived the minimum energy magnetic field and pressure within the lobes. For completeness, we also considered the case where the emission is due to Bremsstrahlung.

4. We derived the jet characteristic scale to verify the applicability of the jet model by Kaiser et al. (2004).

5. We estimated the jet advancing velocity assuming that it equals the velocity of the bow shock in the ISM and derived the jet age.

6. We estimated the jet energy transport rate.

7. Finally, we derived the jet energy transfer rate via two more simplified calorimetry methods based on the estimate of the enthalpy of the jet–ISM system and on a different use of the minimum energy argument previously applied to AGNs.

For all the quantities we derive below, we provide either ranges or lower limits, depending on the assumptions made. Ranges were obtained using a Markov chain Monte Carlo approach in order to derive a posterior distribution for each given quantity, from which we obtained an upper and lower boundary as the values corresponding to the 68% confident range (16th and 84th percentile) unless otherwise specified. This choice was dictated by the large uncertainties on some of the parameters we used or by the need to make conservative assumptions. The calculations behind the treatment below are presented in the accompanying Jupyter notebook.

4.1. Assumptions

In the interest of clarity, in this section we state the assumptions that are made in our analysis. This work hinges on two major assumptions, which have been already introduced above. The first is that the IRAS region emits radiation at least partly due to interaction with a jet from GRS 1915+105. The second assumption is that the arched structure we discovered in the MeerKAT images is the signature of a bow shock induced by the jet from GRS 1915+105, marking the edge of a lobe. Other assumptions we made are listed below and discussed in the following sections:

• The IRAS region is filled by multi-phase gas, i.e. a mix of hot, fully ionised pure hydrogen gas at a temperature between 10⁴ K and 3 × 10⁶ K, and cold, denser molecular gas with temperatures around 10 − 100 K.

• The emission from the hot gas in the IRAS region is predominantly Bremsstrahlung, as indicated by the flat continuum spectrum, (Rodriguez & Mirabel 1998).

• The IRAS region can be modelled as a sphere with a certain filling factor to account for the inhomogeneity of the gas (see Sec. 4.2.1).

• The emission from the northern feature is predominantly synchrotron radiation, as indicated by the steep continuum spectrum from this region, (Rodriguez & Mirabel 1998), with a spectral index α = −0.75.

• The northern feature can be modelled as a cylindrical shape.

• The gas within the northern region is in equipartition (see Sect. 4.2.2).

• The lobe blown by the jet, partly delimited by the bow shock structure, can be modelled as an ellipsoidal shape with a major axis equal to the distance between GRS 1915+105 and the IRAS region and minor axes equal to the bow shock width (see Sect. 4.2.3).

• The jet opening angle is constrained between 0 and 10 degrees.

• The adiabatic indices of the material in the lobe cavity, in the external medium, and within the jet (Γ_j) at the jet–ISM interaction site are assumed to equal 5/3 (see Sect. 4.3.2).

4.2. The structure of the jet–interstellar medium interaction site

4.2.1. The IRAS region

The continuum emission of IRAS 19132+1035 is characterised by a flat radio spectrum, which indicates a Bremsstrahlung origin from a gas emitting at a temperature of T ∼ 10⁴ K. Below this temperature, hydrogen is not significantly ionised and hence cannot produce Bremsstrahlung radiation efficiently. This is further supported by the detection of a H92α recombination line with width ∼25 km s⁻¹ (Rodriguez & Mirabel 1998), which implies a temperature of at least T ∼ 10⁴ K. We assumed a conservative upper limit of T = 3 × 10⁶ K to the temperature of the Bremsstrahlung emitting gas, assuming that the shock forming at the jet end is radiative (and thus the electrons need to have sufficient energy to emit photons effectively when decelerated) and noting that Bremsstrahlung cooling becomes significant at temperatures typically above 10⁶ K for particle densities typical for the ISM and HII regions (0.1–100 cm⁻³; see e.g. Gould 1980).

From the monochromatic emissivity ϵ_ν of ionised hydrogen due to thermal Bremsstrahlung (Longair 1994), we can derive the electron density n_e:

$$\begin{array}{c}\hfill n_{\mathrm{e}}=\sqrt{\frac{{ϵ}_{\nu }}{g(\nu ,T)\sqrt{T}{C}_{\mathrm{radio}}exp(\frac{h\nu }{k_{b}T})}},\end{array}$$(1)

where C_radio = 6.8 × 10⁻³⁸ erg s⁻¹ cm⁻³ Hz⁻¹, and g(ν) is the Gaunt factor (see Appendix A.2 for a detailed treatment). We assume that the gas in the region is pure hydrogen, and thus the number density of electrons n_e equals the number density of protons n, and that the density and temperature are uniform within the volume V considered, which in the case of the IRAS region for simplicity we modelled as a sphere with apparent diameter 36 arcsec. Since the IRAS region contains also cold and dense molecular gas, which does not contribute to the continuum emission, we assume a sphere filling factor of 0.5 to account for the gas inhomogeneity.

Considering a temperature range of T ∼ 10⁴ − 3 × 10⁶ K, we derived an electron density of the shock-compressed gas between 395 and 630 cm⁻³. At the above temperatures, such particle densities imply that hydrogen must be fully ionised, and thus the electron density equals the proton density. We note that our approach differs from that taken in Gallo et al. (2005) for the case of Cyg X-1, who assumed an ionisation fraction of ≈2%, which resulted in a significantly higher total particle density (see also Atri et al. 2025). Assuming that the shock-compressed gas is approximately 4 times denser than the surrounding medium, we deduce that the pre-shock gas density (equal to the un-shocked ISM electron density) is between ≈100 and 160 particles per cm³. This range is consistent with the typical density of an HII region (and over 10 times larger than the typical ISM density). Since the bow shock structure seems to be indeed connecting with an HII region, these density values are reasonable.

4.2.2. The northern feature

Both the MeerKAT radio map and the older VLA radio maps (see Rodriguez & Mirabel 1998) of IRAS 19132+1035 show an elongated feature extending from the main structure towards the direction of GRS 1915+105 for approximately 17.4 arcsec, with a steep radio spectrum and an integrated flux density of 1.57 mJy. In the VLA images the spectrum of the feature is consistent with a power law F ∝ ν^α with α = −0.75, typical of synchrotron emission, as opposed to the flat spectrum of the IRAS region. The synchrotron spectrum implies the presence of a relativistic plasma consisting of a magnetic field and a population of relativistic electrons, for which we can calculate the total energy, and from there infer the magnetic field corresponding to minimum energy, and the associated minimum energy pressure.

We assumed a radio spectrum with some spectral index α between upper and lower frequencies ν₁, ν₂, and α = (1 − p)/2 = −0.75 (Mirabel & Rodríguez 1994). The energy density in magnetic fields is B²/8π (in cgs units); therefore, for a volume V, the total energy in the magnetic field is

$$\begin{array}{c}\hfill E_B=\frac{B^2}{8\pi }fV,\end{array}$$(2)

where f is a filling factor, which represents the fraction of the physical source’s volume that is occupied by the magnetic field and relativistic particles. In our case, we assume f = 1.

Given a population of electrons with a power law distribution with index −p between energy limits E₁ and E₂, the total energy in these electrons can be written as

$$\begin{array}{c}\hfill E_{\mathrm{e}}=c_{12}(p,{\nu }_1,{\nu }_2)B^{-3/2}L,\end{array}$$(3)

where we characterised the synchrotron emission as coming from a single frequency ν = c₁₂BE² for a given electron energy and magnetic field. The pseudo-constant c₁₂ encapsulates all the information about the frequency range and the slope of the energy spectrum between those limits. In Appendix A.2 a more detailed treatment is provided, which includes the explicit form of the c₁₂ constant.

Since the energy stored in the magnetic field is proportional to the magnetic field strength squared, there must be a minimum of the total energy as a function of B, if we assume that there are no other contributions to the energy budget. Such a minimum occurs close to the equipartition of the energy, which corresponds to an equipartition magnetic field B_eq given by

$$\begin{array}{c}\hfill B_{\mathrm{eq}}={\left(6\pi \frac{\eta }{f}{c}_{12}\frac{L}{V}\right)}^{2/7},\end{array}$$(4)

where η indicates how much energy is stored in the protons that accompany the electrons (in the case of a normal baryonic plasma), and a total pressure is given by

$$\begin{array}{c}\hfill p_{\mathrm{eq}}=\frac{7}{9}\frac{B_{\min }^2}{8\pi }(k+1),\end{array}$$(5)

with k equal to the ratio of the internal energy stored in any particles that do not contribute to the synchrotron emission to the sum of the energy in the magnetic field and in the relativistic electrons.

For simplicity, we assume a cylindrical geometry for the northern feature, with a radius of 3.9 arcsec and a length of 17.4 × sin(i)⁻¹ arcsec (Mirabel & Rodríguez 1994), where i is the angle of the long axis of the cylinder to our line of sight (which we assumed to be 60°). Thus, the volume of the northern feature is V_cyl = 2.6 × 10⁵⁴ cm³, and its luminosity L_cyl = 2 × 10²⁹ erg s⁻¹. The predicted minimum energy magnetic field is B_min > 2.2 × 10⁻⁵ G, and the minimum energy pressure is p_Bmin = 1.5 × 10⁻¹¹ erg/cm³.

4.2.3. The lobe

In the case of jets from AGNs, a substantial fraction of the energy transported by the jet is transferred to a population of relativistic electrons and a magnetic field in the lobes, which become bright due to synchrotron emission. The lobes are thus generally easily detected in radio in the case of AGNs, and similar lobes may be expected to emerge around microquasars due to a similar mechanism if the environment density is large enough (see also Heinz 2002). The search for such lobes around GRS 1915+105 has not been successful so far, likely due to the insufficient sensitivity of such searches. In this regard, Kaiser et al. (2004) argued that, under a reasonable set of assumptions, a synchrotron lobe would have had a specific intensity of approximately 0.1 mJy/beam at a resolution of 4 arcsec (i.e. comparable with the MeerKAT resolution), which was below the 1-σ rms noise of the VLA data taken in the ’90s. Such a brightness is remarkably consistent with that we measured in the bow shock structure in the MeerKAT data, i.e. 0.15 mJy/beam.

Let us assume that the bow shock structure corresponds to the brightened limb of the lobe, and let us estimate the magnetic field implied by the measured average flux density and the associated pressure inside the lobe, which we can compare with the pressure within adjacent jet regions. We assume that the conditions of minimum energy hold in the lobe, and so we can infer the amount of energy necessary to generate the detected luminosity, and from there the magnetic field corresponding to the minimum energy criterion B_min using Eq. (4) (see also e.g. Longair 1994). We neglected any energy losses of the relativistic electrons due to radiation processes. The minimum energy magnetic field B_min can be compared with the magnetic field density we obtained for the non-thermal northern feature (see Sect. 4.2.2). Due to the low luminosity of the bow shock we cannot directly measure its spectral slope α. Hence, following Kaiser et al. (2004), based on the fact that the relativistic plasma powering the northern region inflates the cavity, we assume that also the bow shock structure emits synchrotron radiation with a spectral slope α = −0.75.

We model the radio lobe as an ellipsoidal shape, with a major axis equal to the distance between GRS 1915+105 and the outermost edge of the IRAS region, and a semi-minor axis equal to the width of the bow shock in the direction orthogonal to the jet (see Table 1). We obtain a magnetic field density of $B_{\min }^{\mathrm{lobe}}\sim 6.4\times {10}^{-6}$ G, and a minimum energy lobe pressure is p_lobe > 1.3 × 10⁻¹² erg cm⁻³, which are consistent with the values we obtained for the northern feature. The pressure is also significantly larger than the pressure in the unperturbed ISM, as expected based on our assumptions.

If, instead, we assume that the bow shock structure emits Bremsstrahlung radiation similarly to the IRAS region, we can infer the implied electron density using Eq. (1) following the procedure described in Sect. 4.2.1. In this case, we assume that the depth of the ring is of the order of the projected width of the bow shock on the plane of the sky (ΔR) because of limb brightening, and we determine the unit volume associated with the average specific intensity of the region (0.15 mJy/beam) as the beam surface area (6.6 arcsec²) times the ring thickness (60 arcsec, i.e. 2.7 pc). The associated pre-shock electron density ranges between 13 and 26 particles per cm³, i.e. approximately a factor of 6 lower than what we inferred in Sect. 4.2.1.

4.3. Jet calorimetry

The jet model assumes a constant energy transport rate, Q₀, for each side of the jet, and is averaged over the lifetime of the jets. We further assume that the ambient medium has a constant density ρ₀. Under these circumstances the evolution of the large-scale structure created by the jet (lobe and bow shock), is self-similar once the jet extends beyond its characteristic length scale, which is given by Eq. (5) in Kaiser et al. (2004) (see also Falle 1991).

The length scale depends on the density of the gas through which the jet is propagating, and the jet Lorentz factor γ_j. We derived the former in Sect. 4.2.1. The latter is more difficult to estimate (see e.g. Fender et al. 2003), as many assumptions that can scarcely be verified are generally involved. For instance, the compact and steady jets (which are rarely resolved) may be fundamentally different from the discrete relativistic ejections more often observed around X-ray binaries, and hence the motion of radio ‘knots’ in the jets may simply reflect the motion of shocks along the jets rather than the bulk motion of the jet material. Any estimate of γ_j will also depend on the distance to the source and on the orientation of the jet to our line of sight, which, when known, is generally affected by large uncertainties. In the interests of generality, here we use a conservative lower limit γ_j = 1.02 (corresponding to β = 0.2), which maximises L₀.

For a reasonable value of Q₀ ∼ 10³⁶ erg s⁻¹ (see below) and for the minimum density derived in Sect. 4.2.1, L₀ is of the order 10⁻⁵ pc, that is, orders of magnitude smaller than the jet length L_j (> 40 pc), which confirms the applicability of the model by Kaiser & Alexander (1997), as already showed by Kaiser et al. (2004).

4.3.1. Jet velocity and age

Kaiser & Alexander (1997) showed that the velocity of the jet advancing into the ISM is roughly equal to the expansion velocity of the bow shock driven by the jet itself, which implies that the age of the bow shock equals the age of the jet that produced the (most recent) shock. For a strong shock in a monatomic gas, the expansion velocity $\dot{L}$ is set by the temperature of the shocked gas, and is given by

$$\begin{array}{c}\hfill \dot{L}=\sqrt{\frac{16k_{\mathrm{b}}}{3m_{\mathrm{p}}}T},\end{array}$$(6)

where m_p is the mass of the proton. If the shock emits radiation, the initial post-shock temperature can be higher than that of the thermalised, bremsstrahlung-emitting gas. We assume that the temperature of the gas is constrained between 10⁴ and 10⁶ K, based on the similarity with the bow shock structure observed near Cyg X-1 (see Gallo et al. 2005, and accompanying paper by Atri et al. 2025). We estimate a shock velocity between 20 and 360 km s⁻¹, which implies that the bow shock will be highly supersonic with respect to the unshocked ISM.

For an ISM with constant density Kaiser et al. (2004) showed that the growth of the jet length within the lobe as a function of time t is given by

$$\begin{array}{c}\hfill L_j=C_1{\left(\frac{Q_0}{{\rho }_0}\right)}^{1/5}{t}^{3/5},\end{array}$$(7)

where ρ₀ is the density of the un-shocked surrounding medium, and Q₀ is the time-averaged jet power being transferred to the ISM. ρ₀ can be inferred from Eq. (1). C₁ is a dimensionless constant between 2 and 9 and median ≈3 that depends on the thermodynamical properties of the jet material, and on the aspect ratio, R, of the lobe inflated by the jet. The explicit form of C₁ is given in Appendix A.3.1. Taking the time derivative of Eq. (7), one obtains

$$\begin{array}{c}\hfill {\dot{L}}_j=\frac{3L_j}{5t},\end{array}$$(8)

which implies

$$\begin{array}{c}\hfill t=\frac{3L_j}{5\dot{L_j}}.\end{array}$$(9)

From Eq. (8) we infer the jet age by assuming that L_j corresponds to the distance between GRS 1915+105 and the outermost edge of the IRAS region, corrected by the jet inclination angle i with respect to the line of sight. Based on our assumptions, the jet age is between 0.09 and 0.2 Myr, with an absolute upper limit of 1.3 Myr. This value essentially only depends on the shock advancing velocity, which in turn is calculated from the gas temperature T. We note that these estimates are robust as they are set by conservative assumptions, although better constraints could be obtained by refining the temperature range considered in the calculations.

4.3.2. Jet power

Assuming that the direction of the jet remains constant over time, and again that the jet is colliding with a medium of density ρ₀, the power Q_jet being transported by the jets into the surrounding averaged over its lifetime depends on the properties of the jet itself. Combining Eqs. (7) and (8), we obtained the jet power:

$$\begin{array}{c}\hfill Q_{\mathrm{jet}}={\left(\frac{5}{3}\right)}^3\frac{{\rho }_0}{{\overline{C_1}}^5}{L}_j^2{\dot{L_j}}^3.\end{array}$$(10)

This expression depends strongly on the constant $\overline{C_1}$, which encapsulates several properties of the jet, including the jet opening angle θ, which is the main driving parameter, the aspect ratio R of the lobe inflated by the jet, and the adiabatic indices of the jet, lobe cavity, and un-shocked ISM (see Appendices A, Appendix A.3.1). The indices can be assumed to be all equal to 5/3 (see Kaiser & Alexander 1997; Kaiser et al. 2004, for a discussion). Since changes in the indices only mildly affect the value of the constant, we choose to ignore possible variations.

For extragalactic jet sources, R is in general lower than 2 (Kaiser et al. 2004). In the case of the lobe in GRS 1915+105, assuming an ellipsoidal jet lobe, the aspect ratio R of the lobe equals the axial ratio (given by the length of the lobe divided by its width), which gives R = R_ax ≈ 49 pc/25 pc ∼ 2 (see also Appendix A and Appendix A.3). This in principle provides a constraint on the jet opening angle since the variable $\overline{C_1}$ is a function of the opening angle (see Appendix A). In particular, $\overline{C1}\approx 2$ implies a jet opening angle of approximately 30 degrees. This value is large and at odds with the opening angles typically inferred for the transient jets seen in X-ray binaries, which are estimated or measured to be at most 4 degrees (see e.g. Miller-Jones et al. 2006). However, it is unclear whether the large-scale jet interacting with the ISM at large distances from the launching site has the same properties as either the transient intermediate-state jets often observed as a resolved radio-bright outflow expanding outwards on milli-arcsecond to arcsecond scales from the binary position or the continuous hard-state jets resolved at a few milli-arcsecond scales (Fender et al. 2009). Thus, for the sake of generality and to avoid extreme scenarios (i.e. very high or low collimation), we assumed a jet opening angle ranging between 1 and 10 degrees⁹. We caution the reader that the choice of opening angle heavily affects the final jet energy estimate as larger opening angles correspond to increasingly more powerful jets, and thus it constitutes an important caveat in our treatment.

For an opening angle between 1 and 10 degrees, $\overline{C1}$ ranges between 3.5 and 7.5, where the latter value comes from smaller opening angles. The one-sided (only one side of the bi-polar jet is considered) time-averaged jet energy transferred to the environment Q_jet from Eq. (10) is Q_jet ∼ 3.3 × 10³⁷ − 1.5 × 10³⁹ erg s⁻¹, where the upper limit is largely driven by the upper limit imposed to the jet opening angle.

The pressure inside the lobe is given by

$$\begin{array}{c}\hfill p=0.0675\frac{C_1^{10/3}}{R^2}{\left(\frac{{\rho }_0{Q}_0^2}{L_j^4}\right)}^{1/3},\end{array}$$(11)

which ranges between 2.8 × 10⁻¹¹ and 7.8 × 10⁻¹⁰ erg cm⁻³. This means that the lobe pressure is consistent with the pressure we obtain for the cylindrical feature stemming from IRAS 19132+1035 (see Eq. (5)), and significantly overpressured with respect to the un-shocked ISM, as required by the jet model. For a temperature of 50 K, a typical ISM pressure would be of the order 9 × 10⁻¹³ erg cm⁻³ assuming ideal gas conditions.

4.3.3. The ‘Enthalpy’ method

The estimate of the time-averaged energy rate transported by the jet described above hinges on the assumption that the emission from the IRAS region is Bremsstrahlung. In order to confirm our results we derived the energy deposited into the ISM via two more simplified though very generic calorimetry methods.

In the case of AGNs, the jet power transport rate is often estimated by calculating the enthalpy of the regions powered by the jets, and dividing by an appropriate timescale (e.g. Bîrzan et al. 2008; Cavagnolo et al. 2010). A change in enthalpy at constant pressure is given by dH = dU + pdV, and it includes the internal energy stored in a volume and the work done in excavating it. Hence, it measures the total energy/heat ratio input by the jet. The enthalpy of an ideal gas with adiabatic index γ is given by

$$\begin{array}{c}\hfill H=U+PV=\frac{\gamma PV}{(\gamma -1)},\end{array}$$(12)

which for a relativistic ideal gas (i.e. γ = 4/3), as is commonly assumed in the AGN case, gives H = 4PV (for a non-relativistic gas, γ = 5/3 and H = 5PV/2). The jet power is thus given by

$$\begin{array}{c}\hfill Q=\frac{4PV}{t},\end{array}$$(13)

where P is the pressure, which we assumed is equal to the minimum energy condition pressure we calculated above for the lobes (see Sect. 4.2.3), namely, 1.3 × 10⁻¹² erg cm⁻³. The term V is the volume excavated by the jet from GRS 1915+105, and t is the jet age. We assumed that the cavity excavated by the jet is an ellipsoidal shape with a length equal to the source-to-bow shock distance, and a width equal to the bow shock diameter, which has a volume of 6 × 10⁵⁹ cm³. The upper limit to the jet age t is around 0.4 Myr (95th percentile) from the self-similar estimate (see Sect. 4.3.1), which is consistent with the time scale obtained just by dividing the jet length by the shock velocity. The resulting one-sided jet transferred power ranges between 1.2 × 10³⁷ erg s⁻¹ and 1.5 × 10³⁹ erg s⁻¹. Despite the large uncertainty on this estimate, the value we obtained (especially its lower limit) is consistent with our estimate from the self-similar model in Sect. 4.3.2. We note that these values could easily increase due to the presence of a non-radiating pressure component (and/or a non-volume-filling synchrotron plasma) in the lobes.

4.3.4. The ‘Hot-spot’ method

For the case of AGNs, under the assumption of equipartition, Godfrey & Shabala (2013) derived the jet power being supplied to the hot spot, which is given by

$$\begin{array}{c}\hfill Q_{\mathrm{GS}}=Ac\frac{B_{\mathrm{eq}}^2}{8\pi }g=2\pi R^{2}c\frac{B_{\mathrm{eq}}^2}{8\pi }g,\end{array}$$(14)

where g is a function of various parameters of the flow causing the emission, but for simplicity’s sake, we assumed g ≈ 2 based on the AGN case¹⁰. In Eq. (14), A is the area of the hot spot surface, which we assumed is a circular surface with radius R ≈ 7.6 × 10¹⁸ cm based on the area of the northern feature, and B_eq is the equipartition magnetic field, which from Sect. 4.2.2 is approximately 22 μG. We thus have Q_GS ∼ 1.7 × 10³⁶ erg s⁻¹. We note that this value is a lower limit since departures from equipartition will increase it. Still, it is consistent with the lower limit to the power estimate given in Sect. 4.3.2.

5. Discussion

In the MeerKAT observations of the BH X-ray binary GRS 1915+105 we found an extended arch-like structure located to the south-east of the system, extending at either side of the bright source IRAS 19132+1035. IRAS 19132+1035 is located at a distance consistent with that of GRS 1915+105, and was identified as a potential jet–ISM interaction zone attributed to GRS 1915+105 (Rodriguez & Mirabel 1998; Chaty et al. 2001). The association of IRAS 19132+1035 with GRS 1915+105 was recently confirmed by Tetarenko et al. (2018), who mapped the molecular line emission across IRAS 19132+1035, and found evidence that supports the presence of a jet–ISM interaction at this site. As already noticed by Rodriguez & Mirabel (1998), to the north of IRAS 19132+1035 a linear feature with non-thermal spectrum is observable (referred to as the northern feature in this paper). Such a feature has been interpreted as associated with the large-scale jet, an association that was confirmed in Tetarenko et al. (2018) based on molecular line features which indicate that the jet is pushing into dense gas.

5.1. Bow shock structure comparison with Cyg X-1

Based on the similarity in morphology to the bow shock structure observed near Cyg X-1 (Gallo et al. 2005), and on the apparent association with IRAS 19132+1035, we postulated that the arched structure near GRS 1915+105 formed as the result of the interaction of the jets with the ISM. In particular, we interpret the structure as the result of the jets from the binary ending in strong shocks at the location where they interact with the surrounding medium and blow a lobe into the ISM. The gas surrounding the cavity is thus compressed, and this results in a bow shock forming at the jet end.

Compared to the case of Cyg X-1, the jet-induced structure near GRS 1915+105 shows greater complexity, featuring a dim bow shock structure, previously unknown and similar to that discovered near Cyg X-1 (Gallo et al. 2005), a bright, flat-spectrum extended emission region (the IRAS region), and a steep-spectrum linear structure adjacent to it, pointing towards the position of GRS 1915+105. In the jet-related scenario that we depicted, the IRAS region corresponds to compressed and shock-heated ISM in front of the jet, and the northern feature corresponds to the strong shock that originates at the end of the postulated large-scale jet, which in AGNs would be the hot spot. The bow shock structure corresponds to the shock-compressed material at the edges of a jet-blown cavity and was not found in previous searches for similar features (Rodriguez & Mirabel 1998; Heinz 2003). Under the assumption that the structure near Cyg X-1 and GRS 1915+105 are both jet-induced, the main difference between the two – i.e. the absence of a hot spot in Cyg X-1 – could be explained in terms of particle acceleration at the jet end, which is efficient in GRS 1915+105, and inefficient or possibly absent in the case in Cyg X-1. It is also worth noting that the jet of GRS 1915+105 terminates at a molecular cloud, and hence it is presumably ramming into a denser medium than in the case of Cyg X-1, where the jet seems to interact with a less dense molecular cloud.

5.2. Properties of the interstellar medium

Beside IRAS 19132+1035, also a candidate jet–ISM interaction spot with the counter-jet was identified to the north-west of GRS 1915+105 (IRAS 19124+1106). However, we do not detect any radio lobe near this region. If IRAS 19124+1106 is indeed an interaction site with the jet, then a possible bow shock structure (and lobe) would overlay a bright extended HII region (which IRAS 19124+1106 may or may not be part of). In either case, it would be hard to identify a possible bow shock structure, assuming its brightness would be similar to the structure at the south-east of GRS 1915+105.

From the flux density of the thermal emission from the shock compressed gas we inferred the density of the ISM in the region of the shock. We find a particle density of the order 100 particles per cm³, which is relatively high, but not excessively so given that this region of the sky is located on the Galactic plane, and is rich in dense star-forming regions. The density we found supports the hypothesis that the bow shock formed because of the interaction of the jets from GRS 1915+105 with a large HII region, similar to what was proposed for the case of Cyg X-1. This would help reconcile the apparent contrast between the gas density we derive, and the low values predicted by, for example, Heinz (2002, 2003), Carotenuto et al. (2024), who suggest that X-ray binaries may be located in low-density cavities in the ISM.

5.3. Jet properties

Based on our conservative estimate of the temperature of the emitting gas within the IRAS region, we inferred the jet age and the jet time-averaged energy transport rate. This was under the assumption that the IRAS region is dominated by Bremsstrahlung emission indicated by its flat spectrum. We estimate that the jet that generated the bow shock must be between 0.09 and 0.22 Myr old, with an upper limit of ≈0.4 (95th percentile; see Sect. 4.3.1). If our assumptions are correct, the time required for the jets from GRS 1915+105 to slice across the IRAS region as the system moves through the Galaxy should be longer than the jet age itself. On the other hand, if the projected motion of GRS 1915+105 implies that its jet crossed the IRAS region faster than the age of the jet calculated using the temperature of the gas, the motion of GRS 1915+105 will give a stronger constraint on the jet age suggesting that this might be a better indicator of jet age than the analysis presented in this work. The projected proper motion of GRS 1915+105 on the plane of the sky is oriented approximately north-east to south-west (top left to bottom right in the image in Fig. 1), i.e. essentially orthogonal to the direction of its jets (Dhawan et al. 2007). Thus, we assume that the time GRS 1915+105 will need to move across the IRAS region is similar to the time the jets from GRS 1915+105 could have interacted with the IRAS region. Estimates of the potential kick velocity of GRS 1915+105 show that the BH in the system was born with a direct collapse, and thus the system’s peculiar velocity is comparable to the local standard of rest (Atri et al. 2019; Reid & Miller-Jones 2023). Considering these two facts, and using the most recently reported non-circular velocity of GRS 1915+105 of 20 km s⁻¹, we estimate that the time GRS 1915+105 (along with its jets) will take to move across the IRAS region is approximately 0.9 Myr, assuming a width for the IRAS region of 7.2 arcmin. As expected, this time period is larger than the jet age, from which we can estimate a much more stringent lower limit on the jet power.

The one-sided time-averaged energy transport rate ranges between 3.3 × 10³⁷ and 1.5 × 10³⁹ erg s⁻¹, which we derived using the self-similar model proposed by Kaiser et al. (2004). This estimate is confirmed by independent though simplified calculations based on arguments usually applied to AGNs. The first method is based on the estimate of the enthalpy of the regions powered by the jets, and the second yields the jet power under the assumption of equipartition applied to the jet hot spot. In both cases, we obtain a strict lower limit to the jet power of 1.7 × 10³⁶ erg s⁻¹. The jet power estimated above exceeds by at least three orders of magnitude the radiative luminosity¹¹ of the nebula, ∼1 × 10³⁴ erg s⁻¹. This substantial difference indicates that the jet supplies far more energy than is required to power the observed radiation from the nebula. The excess power inferred from the jet’s energy budget is likely deposited in ongoing particle acceleration and other processes that drive the physical evolution of the nebula. A similar situation is observed in the case of Cyg X-1 (Atri et al. 2025), where the bow shock radiative luminosity (∼3 × 10³¹ erg s⁻¹ Hz⁻¹) is relatively low if compared to the jet power (between ∼1 × 10³⁵ ergs s⁻¹–1 × 10³⁸ ergs s⁻¹).

Our estimates of the time-averaged jet power differ quite significantly from the estimate derived by Tetarenko et al. (2018) (∼10³² erg s⁻¹), which was also obtained following Kaiser et al. (2004) and arguments similar to those we adopted in this paper. The difference in the estimate can be ascribed to slightly different assumptions on the jet properties and on the composition of the IRAS region, and in particular (i) the jet opening angle (0–4 degrees in Tetarenko et al. (2018) versus 1–10 degrees in this work), (ii) the density within the IRAS region (4000 particles/cm³ versus ∼80–180 particles/cm³), (iii) the jet front expansion velocity (1 km s⁻¹ versus 20–360 km s⁻¹). The assumption on the jet opening angle is rather speculative and does not affect significantly the lower limit to the time-averaged energy transport rate (set by the lower limit on the jet opening angle). The particle density in Tetarenko et al. (2018) is estimated based on the CO density maps from ALMA data, which are used to estimate the HII column density. This value is necessarily associated with the cold gas in the region, which in this work is postulated to fill only 50% of the IRAS region, while the remaining volume is filled with hot gas (see Sect. 4.2.1), which is generating the radio continuum emission considered here. The highest impact (a difference of the order 10⁴), is given by the shock front propagation velocity, which in Tetarenko et al. (2018) is derived based on the equivalent width of the SiO lines, also arising from the cold gas, and thus not necessarily correct for the hot gas in the region.

5.4. The IRAS region

The existence of multi-phase gas with temperatures of ∼10 − 100 K and > 10⁴ K, as already suggested by Kaiser et al. (2004) is required by the simultaneous presence of both the flat radio continuum and the lines emission. The ISM is inherently multi-phase (Tielens 2005) – but how does this cold phase survive being shocked? This apparent conundrum can be resolved by considering the different shock propagation velocities through the two gas phases. If the shock strikes a multi-phase ISM region in pressure equilibrium, then the density contrast χ between the molecular clouds and warm phase is the inverse of their temperature ratio, implying χ ∼ 10² − 10³. The question of how long the molecular clouds survive is critical and corresponds to the well-studied ‘cloud-crushing’ problem, where the time for a shock to pass through an overdense cloud (in our case the cold, denser molecular gas) is obtained from ram pressure balance (Klein et al. 1994) and given by

$$\begin{array}{c}\hfill {\tau }_{\mathrm{cc}}={\chi }^{1/2}{R}_{\mathrm{c}}/v_{\mathrm{s}},\end{array}$$(15)

where v_s is the shock velocity in the underdense, warm region, and R_c is the radius of a dense cloud. Hence, we argue that the molecular line emission is observed because (i) the clouds have not had sufficient time to be shocked and disrupted (perhaps implying that R_c is fairly significant compared to the width of IRAS 19132+1035) and/or (ii) the lower shock velocity within the molecular clouds means the post-shock temperature is much lower than that in the warm phase. The basic scenario is that of a shock passing quickly through the warm phase and leaving behind embedded dense molecular clouds that are disrupted in a few cloud-crushing times. This behaviour is essentially that depicted in McKee & Ostriker (1977) (Fig. 2), and also occurs in simulations of the interaction of AGN jets with a two-phase ISM (Mukherjee et al. 2018).

The morphology of the IRAS region in radio and in the sub-millimetre band and the presence of cold gas and molecular lines that can be ascribed to young star activity prompted Tetarenko et al. (2018) to evaluate what influence star formation had on the region and consider the possibility that star formation was triggered by the jet itself interacting with the environment. While the thermal component of the radio continuum could be explained in terms of the activity of a young star cluster, we argue that the discovery of a bow shock marking the presence of a cavity inflated by the jets from GRS 1915+105 suggests that the jet’s impact on the region must have been significant. Hence, the star activity might not have been the dominant factor in shaping the IRAS region. As discussed in Tetarenko et al. (2018) and other authors (see e.g. Podsiadlowski et al. 2003), GRS 1915+105 is old enough that the jet activity might have indeed triggered star formation. The transient nature of GRS 1915+105 and the age of the jet estimated above imply that the bow shock structure we observe today may only be the last one formed by the jet, and much older activity might have occurred. We note in particular the presence of arched structures to the south-east of the IRAS region, following the same orientation of the bow shock structure. Such structures may be related to older bow shocks formed in previous bursts of jet activity, which expanded further out with respect to the IRAS region, and may be between 0.1 and 0.5 Myr old.

5.5. Implications

In the scenario we have depicted, the energy rate transported by the large-scale jet that induced the formation of the bow shock is comparable with the average energy typical of the small-scale jets observed during outbursts. The average energy rate is estimated to be between 10³⁷ and 10⁴² erg s⁻¹ for both the quasi-continuous jets typical of the harder states, and the transient discrete jet observed at sub-arcseconds scales during the intermediate states (Fender et al. 1999; Fender & Pooley 2000; Heinz & Grimm 2005). Especially the latter is believed to be launched on short time scales, which may imply that the dark jet either transports similar amounts of energy as the transient jets and is active over similarly short scales or it transports less energy but is active for longer times. In this regard, it might be possible that the presence of a hot spot in the bow shock structure near GRS 1915+105 – a feature that is absent in the case of Cyg X-1 – might be connected with the possible turning on and off of the jet activity, as opposed to the quasi-continuous jet from Cyg X-1 (Cyg X-1 also goes through occasional soft states). It is also interesting to notice that the minimum time-averaged energy transferred by GRS 1915+105 is at least an order of magnitude larger than what is estimated for the case of Cyg X-1 (1.4 × 10³⁵ ergs s${}^{-1}<Q_{\mathrm{jet}}^a<1.4\times {10}^{36}$ ergs s⁻¹, (Atri et al. 2025; Gallo et al. 2005).

Since GRS 1915+105 is a transient source and not a persistent one as Cyg X-1, it is hard to correct the time-averaged energy rate transported by the dark jet we estimate above to account for the varying activity period of the system. The estimated age of the bow shock near GRS 1915+105 is at least 80 kyr, which indicates that the structure must have formed much earlier than during the current outburst (started in 1995), that is, either during previous outbursts (undetected, because they were precedent to the birth of X-ray astronomy) or during quiescence. In the former case, the time-averaged transferred energy in GRS 1915+105 may be up to two orders of magnitude larger than what we estimated, as it must take into account the fact that an outburst that may last around 30 years is between 25 and 100 times shorter than the age of the jet that generated the bow shock. If, instead, the bow shock formed during quiescence, then it must be connected with the radio-faint continuous jets observed in these phases, when the accretion energy dissipated as X-ray is negligible (Fender et al. 2004). If this is the case, our estimate of the time-averaged transferred energy would be more accurate.

6. Conclusions

We have reported the discovery of a bow shock structure towards GRS 1915+105 observed in recent MeerKAT observations at 1.28 1GHz. Assuming the structure is associated with GRS 1915+105, we attributed it to the presence of a large-scale dark jet that is interacting with the local environment several parsecs away from the launching site. This finding provides further direct evidence of jet–ISM interactions, a crucial ingredient of BH feedback, the understanding of which has been based so far mainly on AGN studies. Our analysis reveals that the energy dissipated by the jets in this system is comparable to the energy liberated via accretion, indicating that stellar-mass BHs can transfer substantial amounts of energy back into their environments through jet activity.

Our findings have important implications for the study of X-ray binaries as a population. The detection of jet-induced large-scale structures first near Cyg X-1 (a persistent system) and now in the vicinity of GRS 1915+105 (a very variable and transient system, though it has been in outburst for 30 years at the time of writing) demonstrates that the impact of jets from Galactic compact objects may extend far beyond their immediate vicinity, potentially influencing the ISM by enriching it with matter, re-heating it with energy, and even potentially triggering star formation, hence significantly contributing to the broader ecology of the Galaxy.

In the broader context of accreting BHs at all scales, our results indicate that jet-induced feedback mechanisms are not exclusive to SMBHs but are also relevant for their stellar-mass counterparts and not necessarily only those hosted in X-ray binary systems. This expands our understanding of how BHs of different masses impact their cosmic environments, from local star formation to the evolution of galaxies.

Data availability

The un-calibrated MeerKAT visibility data presented in this paper are publicly available in the archive of the South African Radio Astronomy Observatory at https://archive.sarao.ac.za. The capture blocks considered in this work are listed in Appendix A. The Continuum MeerKAT observations were taken as part of the ThunderKAT Large Survey Program, project code SCI-20180516-PW-01. The python source code used to perform the calculations is available online on GitHub https://github.com/saramotta/DarkJets.

Acknowledgments

The authors would like to thank Dave Russell, David Williams-Baldwin, Cristina Baglio, Fraser Cowie, and Andrew Hughes for useful feedback on this work, and Evangelia Tremou and Joe Bright for scheduling the observations that led to the results reported in this paper. SEM acknowledges support from the INAF Fundamental Research Grant (2022) EJECTA. We acknowledge the use of the Ilifu cloud computing facility – www.ilifu.ac.za, a partnership between the University of Cape Town, the University of the Western Cape, Stellenbosch University, Sol Plaatje University, and the Cape Peninsula University of Technology. The Ilifu facility is supported by contributions from the Inter-University Institute for Data Intensive Astronomy (IDIA – a partnership between the University of Cape Town, the University of Pretoria and the University of the Western Cape), the Computational Biology division at UCT and the Data Intensive Research Initiative of South Africa (DIRISA). PA is supported by the WISE fellowship program, which is financed by NWO. JvdE acknowledges a Warwick Astrophysics prize post-doctoral fellowship made possible thanks to a generous philanthropic do-nation. This work made use of the CARTA (Cube Analysis and Rendering Tool for Astronomy) software (DOI 10.5281/zenodo.3377984 – https://cartavis.github.io).

References

Appendix A: Additional information

A.1. MeerKAT observations

A.2. The IRAS region

The monochromatic emissivity ϵ_ν of ionised hydrogen due to thermal bremsstrahlung is given by (Longair 1994):

$$\begin{array}{c}\hfill {ϵ}_{\nu }=\frac{L_{\nu }}{V}=C_{\mathrm{radio}}g(\nu ,T)\frac{n_e^2}{\sqrt{T}}exp(\frac{h\nu }{k_{b}T})erg{s}^{-1}c{m}^{-3}H{z}^{-1}\end{array}$$(A.1)

where C_radio = 6.8 × 10⁻³⁸erg s⁻¹ cm⁻³ Hz⁻¹, and the Gaunt factor g(ν) is given by:

$$\begin{array}{c}\hfill g(\nu ,T)\approx \frac{\sqrt{3}}{2\pi }[ln\left(\frac{128{ϵ}_0^2{k}^3{T}^3}{m_e{e}^4{\nu }^2{Z}^2}\right)-{\gamma }^2]\end{array}$$(A.2)

From Eq. A.2 we can derive the electron density n_e:

$$\begin{array}{c}\hfill n_e=\sqrt{\frac{{ϵ}_{\nu }}{g(\nu ,T)\sqrt{T}{C}_{\mathrm{radio}}exp(\frac{h\nu }{k_{b}T})}}\end{array}$$(A.3)

where we assumed that the number of electrons n_e equals the number of protons n (i.e. pure hydrogen gas), and that the density and temperature are uniform within the volume V considered, i.e. in the case of the IRAS region a sphere with apparent diameter 36 arcsec.

A.3. The northern feature

To model the non-thermal region, we assume a radio spectrum with a spectral index α between upper and lower frequencies ν₁, ν₂. We used the integrated radio luminosity between these limits (and not the specific radio luminosity at a certain frequency¹²) given by

$$\begin{array}{c}\hfill L=4\pi D^2{\displaystyle {\int }_{{\nu }_1}^{{\nu }_2}}{F}_{\nu }d\nu =4\pi D^2{F}_{{\nu }_2}{\nu }_2^{-\alpha }\left(\frac{{\nu }_2^{\alpha +1}-{\nu }_1^{\alpha +1}}{\alpha +1}\right),\end{array}$$(A.4)

where we have assumed α = (1 − p)/2 = -0.75 (Mirabel & Rodríguez 1994).

The energy density in magnetic fields is B²/8π (in cgs units); therefore, for a volume V, the total energy in the magnetic field is

$$\begin{array}{c}\hfill E_B=\frac{B^2}{8\pi }fV,\end{array}$$(A.5)

where f is a filling factor, which represents the fraction of the physical source’s volume which is actually occupied by the magnetic field and relativistic particles. In our case f = 1.

Given a population of electrons with a power law distribution with index −p between energy limits E₁ and E₂, the total energy in these electrons is given by

$$\begin{array}{c}\hfill E_e=c_2^{-1}L{B}^{-2}\frac{(p-3)}{(p-2)}\frac{(E_1^{2-p}-E_2^{2-p})}{(E_1^{3-p}-E_2^{3-p})}.\end{array}$$(A.6)

For a given electron energy and magnetic field, we can characterise the synchrotron emission as coming from a single frequency ν = c₁BE², that is, E = ν^1/2c₁^−1/2B^−1/2. So, we can write

$$\begin{array}{c}\hfill E_e=c_2^{-1}{c}_1^{1/2}\stackrel{\sim }{c}(p,{\nu }_1,{\nu }_2)B^{-3/2}L=c_{12}(p,{\nu }_1,{\nu }_2)B^{-3/2}L,\end{array}$$(A.7)

where

$$\begin{array}{c}\hfill c_{12}=c_2^{-1}{c}_1^{1/2}\stackrel{\sim }{c}(p,{\nu }_1,{\nu }_2)\end{array}$$(A.8)

and

$$\begin{array}{c}\hfill \stackrel{\sim }{c}(p,{\nu }_1,{\nu }_2)=\frac{(p-3)}{(p-2)}\frac{{\nu }_1^{(2-p)/2}-{\nu }_2^{(2-p)/2}}{{\nu }_1^{(3-p)/2}-{\nu }_2^{(3-p)/2}}.\end{array}$$(A.9)

Since E_B ∝ B², the total energy is given by

$$\begin{array}{c}\hfill E_{e+B}=E_e+E_B=const.B^{-3/2}+const.B^2,\end{array}$$(A.10)

and thus there must be a minimum of the total energy as a function of B if we assume that there are no other contributions to the energy budget. Such a minimum occurs at the equipartition of the energy, that is, when E_e and E_B are equal (which can be calculated by differentiating Eq. A.10 with respect to B):

$$\begin{array}{c}\hfill E_B=\frac{3}{4}\eta E_e,\end{array}$$(A.11)

where η indicates how much energy is stored in the protons which accompany the electrons (in the case of a normal baryonic plasma). For the minimum energy condition η = 1 is generally assumed. Using Eqs. A.5, A.7, and A.11, we can express the equipartition magnetic field as

$$\begin{array}{c}\hfill B_{\mathrm{eq}}={\left(6\pi \frac{\eta }{f}{c}_{12}\frac{L}{V}\right)}^{2/7}.\end{array}$$(A.12)

For minimum energy conditions, the total pressure is then given by

$$\begin{array}{c}\hfill p=\frac{7}{9}\frac{B_{\mathit{min}}^2}{8\pi }(k+1),\end{array}$$(A.13)

with k equal the ratio of the internal energy stored in any particles that do not contribute to the synchrotron emission to the sum of the energy in the magnetic field and in the relativistic electrons.

A.3.1. The $\overline{C1}$ dimensionless constant

The term $\overline{C1}$ is a dimensionless constant that depends on the thermo-dynamical properties of the jet material, and on the axial ratio R_ax of the lobe inflated by the jet, which equals the length of the lobe divided by its width if the lobe is assumed to be cylindrical. Such a constant can be calculated by making assumptions on the jet properties and on the environment, which in our case is treated as uniform in density (i.e. the density profile index β equals 0). The constant depends heavily on the jet opening angle θ, as well as on the adiabatic indices of the material in the lobe cavity (Γ_c), in the external medium (Γ_x), and within the jet (Γ_j).

The explicit form of $\overline{C1}$ is the following:

$$\begin{array}{c}\hfill \overline{C_1}={\left(\frac{\overline{C_2}}{\overline{C_3}{\theta }^2}\frac{({\mathrm{\Gamma }}_x+1)({\mathrm{\Gamma }}_c-1){(5-\beta )}^3}{18[9\{{\mathrm{\Gamma }}_c+({\mathrm{\Gamma }}_c-1)\frac{\overline{C_2}}{4{\theta }^2}\}-4-\beta ]}\right)}^{\frac{1}{(5-\beta )}},\end{array}$$(A.14)

with

$$\begin{array}{c}\hfill \overline{C_2}={(\frac{({\mathrm{\Gamma }}_c-1)({\mathrm{\Gamma }}_j-1)}{4{\mathrm{\Gamma }}_c}+1)}^{\frac{{\mathrm{\Gamma }}_c}{({\mathrm{\Gamma }}_c-1)}}\frac{({\mathrm{\Gamma }}_j+1)}{({\mathrm{\Gamma }}_j-1)}\end{array}$$(A.15)

and

$$\begin{array}{c}\hfill \overline{C_3}=\frac{\pi }{4R_{\mathrm{ax}}^2},\end{array}$$(A.16)

where R${}_{\mathit{ax}}=\sqrt{\frac{1}{4}\frac{\overline{C_2}}{{\theta }^2}}$.

All Tables

All Figures

	Fig. 1. Field of GRS 1915+105 as seen by the MeerKAT interferometer at 1.28 GHz in a total on-source time of 15.5 hr. The restored beam for this map is circular and has a radius of 6.6″. The circle marks the position of GRS 1915+105, and the two arrows mark the direction of the jets identified in the ’90s (Fender et al. 1999). In the image, north is up and east is left.
In the text

	Fig. 2. Zoom-in of the jet–ISM interaction region. The bow shock structure is marked by the cyan region, which has been defined by eye. The position of GRS 1915+105 and IRAS 19132+1035, as well as the jet direction, are marked.
In the text

	Fig. 3. Panel A: schematic picturing the jet–ISM interaction region as observed in the MeerKAT image. Panel B: sketch of the model employed to interpret the observed structure.
In the text