© The Authors 2025

1. Introduction

Supermassive black holes (SMBHs; ≥10⁵ M_⊙) at the centres of some active galactic nuclei (AGNs) and stellar-mass black holes (≤100 M_⊙) in binary systems are known to launch relativistic jets. These jets transport a large amount of energy from the nucleus back into the galaxy, thus coupling the growth of the galaxy to the jets (e.g. McNamara et al. 2005; Croton et al. 2006; Nesvadba et al. 2010; Krause 2023). Black hole X-ray binaries (BHXBs), binary systems of stellar-mass black holes in accreting systems with main sequence stars, are considered the lighter analogues of SMBHs. The impact of BHXB jets on the interstellar medium (ISM) is not easily observable, unlike the interactions of AGN jets with the intergalactic medium (IGM). This is due to the different environments that AGNs and BHXBs reside in (Heinz 2002). Unlike AGNs, BHXBs do not reside in dark matter halos and are susceptible to natal kicks (Repetto et al. 2012; Atri et al. 2019). BHXBs can thus be found in the Galaxy with high space velocities and can constantly change their environment. Additionally, the ISM provides a lower barrier to the high-ram-pressured jets from BHXBs, and this results in weaker shock fronts. However, studies have noted the effect of BHXB jets on their local and Galactic-scale surroundings. Heinz et al. (2008) concluded that the plasma released by BHXB jets can fill the disc and halo of our Galaxy. This magnetised plasma can introduce sufficient seed magnetic fields in the ISM to produce the average magnetic field of our Galaxy. The plasma from jets that mixes with the ISM also gives rise to cosmic rays in the ISM (Alfaro et al. 2024). Thus, accurate estimates of the energy injected into the ISM by the relativistic jets of BHXBs are essential for modelling the Galactic magnetic field and the contribution of BHXBs to Galactic cosmic rays.

Active galactic nuclei and BHXBs launch relativistic jets that facilitate the transfer of the accreted power back to the inter-cluster medium and ISM, respectively. One way to understand the relation between accreting matter and the mechanism by which this energy is injected back into the environment is to study the interaction of the radio jets with their surroundings. Estimating various jet properties like the jet speed, the total jet power, the composition of the jet, and the magnetic field at the base of the jet is a challenging task, influenced by a range of assumptions, and has been attempted using a few different techniques (e.g. Russell & Fender 2008; Tetarenko et al. 2018; Fender & Bright 2019). The interaction sites and bow shocks around BHXBs are a probe into many of these jet properties (e.g. Burbidge 1959; Heinz et al. 2006).

Over the decades, gradual improvements have been made in modelling how the jets of AGNs and BHXBs interact with the ISM. These interactions give rise to observable features like radio lobes, radio bow shocks, and large X-ray cavities in AGNs (e.g. Scheuer 1974; McNamara & Nulsen 2007). Early models suggested that radio lobes and hot spots in FR II galaxies require continuous energy input, leading to the idea that jets inflate cocoons that drive bow shocks (Longair et al. 1973; Blandford & Rees 1974). Relativistic jets launched from BHXBs have also been observed to create radio structures as they interact with the ISM (e.g. Gallo et al. 2005; Panferov 2017). The jets launched from a BHXB are expected to sweep up the gas that is in their path, creating a cavity as they propagate from their launch location. The ejected jet plasma can continue into the medium as long as the swept-up mass is comparable to the mass of the ejected plasma and the energy losses of the jet are minimal during propagation. In some cases, the ejecta can propagate farther from the black hole system even if these conditions are not met due to previously ejected material already tunnelling a channel in the ISM. The jet material ultimately reaches the end of this tunnel and interacts with the denser ISM. This creates inflated radio lobes or bubbles far from the base of the black hole jet, and these can lead to bow shocks that expand sideways and forwards from these radio lobes.

Since BHXBs are often found in lower-density environments than AGNs in a dynamical sense, bright jet–ISM interactions require local density amplifications (Heinz 2002; Carotenuto et al. 2024). For example, the nebular structure observed around SS433 is modelled as a combination of a dense supernova remnant and the jets from SS433 creating shocks when interacting with this high-density remnant (e.g. Begelman et al. 1980; Panferov 2017). More recently, a bow shock structure has been discovered near GRS 1915+105 (Motta et al. 2025) adjacent to a large H II region. The bow shock structure detected near Cyg X-1 is close to the tail of a dense H II region, the Tulip Nebula. The Tulip Nebula is considered to be the leading cause of an increase of the local ISM density that has made the jet–ISM interaction near Cyg X-1 detectable as a bow shock structure (Gallo et al. 2005).

Cygnus X-1, or Cyg X-1, is a high-mass BHXB at a distance of 2.2 ± 0.2 kpc consisting of a massive 21.2 ± 2.2 M_⊙ black hole with a blue supergiant of mass 41 ± 8 M_⊙ in a 5.6-day orbit (Miller-Jones et al. 2021). Cyg X-1 was first discovered in 1971 (Hjellming & Wade 1971) and has been the subject of detailed studies and observations owing to its bright nature. BHXBs are jet-dominated in their hard X-ray state, wherein the radio emission is thought to be driven by partially self-absorbed jets (Fender et al. 2004). For Cyg X-1, the hard X-ray state radio emission has been resolved as a collimated jet of size 30 AU (Stirling et al. 2001). BHXBs also have an intermediate state, which is associated with the BHXB leaving the hard X-ray spectral state and moving towards the soft state. During this state, ejecta from the system that is are causally disconnected from the black hole are observable at radio wavelengths (Fender et al. 2006). The energy injected into the ISM by a BHXB is thus due to jet activity in both these states.

Gallo et al. (2005, G05 hereafter) observed Cyg X-1 with the Westerbork Synthesis Radio Telescope (WSRT) at 1.4 GHz and discovered a shock bubble being blown by the jets of Cyg X-1. G05 detected a ring-shaped, diffuse gas structure north of Cyg X-1 with a diameter of 5 pc. The ring was seen in a direction perpendicular to the direction of the motion of Cyg X-1; this, combined with the fact that Cyg X-1 is one of the slowest moving X-ray binaries (XRBs) in our Galaxy (9 ± 2 mas yr⁻¹), allows us to confidently rule out that the structure detected was a shock left in the wake of fast-moving objects. Additionally, the motion of Cyg X-1 cannot be traced back to the centre of the ring, ruling out that the structure could be the supernova remnant of the black hole in Cyg X-1. Later, optical spectroscopic studies of the bow shock structure also confirmed that the origin of the emission could not be related to a supernova remnant (Sell et al. 2015). The direction of the hard state jet of Cyg X-1 is consistent with the scenario that the jet is responsible for creating the bow shock (Stirling et al. 2001). Combined with H_α observations, it was clear that the ring was emitting Bremsstrahlung radiation, and this thermally emitting ring was estimated to be moving with a velocity of 100–360 km s⁻¹. With a well-informed assumption that Cyg X-1 spends 90% of its lifetime in a hard state (Grinberg et al. 2013), the total power of the jet averaged over the lifetime of Cyg X-1 lies in the range 9 × 10³⁵ ≤ P_jet ≤ 10³⁷ ergs s⁻¹ (G05).

The WSRT is an east-west array with a maximum baseline of 2.7 km, and its resolution is twice as high in the east-west direction as compared to the north-south direction for the declination of Cyg X-1 (35°). The map of the ring around Cyg X-1 in G05 has a resolution of 25 × 14 arcsec², and the noise level around the ring is of the order of 70 μJy bm⁻¹, for a total observing time of 60 hours. A more precise measurement of the total jet power, the velocity of the ring, and the average energy input to the ISM requires higher-precision estimations of the observed parameters of the ring-like structure. The MeerKAT radio telescope provides the required resolution and sensitivity to improve the WSRT radio map of Cyg X-1 in the L band (1.4 GHz). With a maximum baseline of 8 km, a total of 64 dishes, and a condensed core, this array provides a great opportunity to image diffuse emission structures like that of the Cyg X-1 bow shock. We obtained a deeper and higher-resolution map of the structure using the MeerKAT telescope in two frequencies, produced the first spectral index map of the bow shock structure, combined these observations with archival optical data to get more accurate estimates of the ISM density, and ultimately modelled the jet parameters of Cyg X-1 that could be responsible for creating the bow shock structure.

The paper is structured as follows: Section 2 contains the description of the observations that were carried out, Sect. 3 lays out the process followed to reduce the data, Sect. 4 presents the results and the image maps obtained, and Sect. 5 discusses the implications of the maps. Finally, we conclude this paper with Sect. 6.

2. Observations

We observed the field of Cyg X-1 using MeerKAT in two frequency bands, the L band and the S band under the ThunderKAT and the XKAT programmes, respectively. The ThunderKAT programme was a Large Survey Program that conducted weekly observations of active transient fields for the years 2018–2023. The XKAT programme is a MeerKAT open call project, which focused solely on observing fields of active XRBs and ran for the years 2023–2024. The L-band observations were conducted on 16 April 2022 (observation block ID 1650088876) with a total bandwidth of 856 MHz centred at a frequency of 1284 MHz. The pointing centre of the observations was right ascension (RA, J2000) ${19}^{\mathrm{h}}{58}^{\mathrm{m}}21\stackrel{\text{s}}{.}67$ and declination (Dec, J2000) of ${35}^{\mathrm{d}}{11}^{\mathrm{m}}50\stackrel{\text{s}}{.}70$. The field was observed for 15 minutes and was bracketed by a 2-minute observation of the phase and amplitude calibrator, J2015+3710 (RA ${20}^{\mathrm{h}}{15}^{\mathrm{m}}28\stackrel{\text{s}}{.}73$ and Dec ${37}^{\mathrm{d}}{10}^{\mathrm{m}}59\stackrel{\text{s}}{.}50$). We also used J1939-6342 (RA ${19}^{\mathrm{h}}{39}^{\mathrm{m}}25\stackrel{\text{s}}{.}03$ and Dec $-{63}^{\mathrm{d}}{42}^{\mathrm{m}}45\stackrel{\text{s}}{.}60$) as a bandpass and flux calibrator, which was observed for 5 minutes at the beginning of the observing block.

The S-band observations were carried out on 14 March 2024 (Observation block ID 1710398476) at a central frequency of 2625 MHz with a total bandwidth of 875 MHz. The pointing centre of the observations was RA ${19}^{\mathrm{h}}{58}^{\mathrm{m}}14\stackrel{\text{s}}{.}23$ and Dec ${35}^{\mathrm{d}}{16}^{\mathrm{m}}22\stackrel{\text{s}}{.}10$, and the field was observed for 30 minutes. The target observations were bracketed by a 2-minute observation of a phase and amplitude calibrator, J2236+2828 (RA ${22}^{\mathrm{h}}{36}^{\mathrm{m}}22\stackrel{\text{s}}{.}47$ and Dec ${28}^{\mathrm{d}}{28}^{\mathrm{m}}57\stackrel{\text{s}}{.}40$). The bandpass and flux calibrator used were the same as for the L-band observations.

3. Data reduction

MeerKAT data observed as a part of the ThunderKAT and XKAT programmes are reduced using oxkat (Heywood 2020), a series of Python code written as wrappers around common radio astronomy software to facilitate a well-calibrated and quick image of the observed field. Depending on the observed field and requirements of the object of interest (e.g. compact or extended emission), further adjustments to the input parameters to these Python scripts are made (e.g. negative or positive Briggs weighting). For both L-band and S-band data reduction on the ILIFU/idia cluster, the process follows three stages.

The first stage (1GC.py) averages the measurement set (MS) into 1024 frequency channels and creates a new MS. This is followed by applying basic flagging, for example, to known bad frequency channels, and edges of the frequency bandpass. flagdata, an autoflagging algorithm in CASA (CASA Team 2022) is also used to identify outliers and flag them in the calibrator data. Once the basic flagging is done, setjy is used to set the flux density scale using the flux calibrator (J1939-6342) and obtain the delay, bandpass and gain solutions. This model is applied to both the phase calibrator (J2015+3710 for the L band and J2236+2828 for the S band) and the target. The amplitude and phase solutions are then obtained using the nearby phase calibrator and applied to the target field (Cyg X-1), which is then split into a new corrected MS.

The second stage (FLAG.py) runs the Tricolor package, which is specifically written for MeerKAT to do automatic Radio Frequency Interference mitigation in the MSs. The first unconstrained image of the fields is then created using wsclean, in that no masks or models are used in creating the image. We used a multi-scale approach for the cleaning algorithm that used Briggs weighting of 0.3 as the best compromise between sensitivity and resolution. An initial cleaning mask that encapsulates most of the bright emission in the field is created using a local rms threshold, and the level of this local rms can be set as required (see the config.py file uploaded on GitHub for the exact input parameters)¹.

The third stage of the scripts (2GC.py) focuses on direction-independent self-calibration of the corrected visibilities from the FLAG.py stage. The initial mask produced by FLAG.py is used to do a masked deconvolution of the data. After this stage, Cubical (Kenyon et al. 2018) is used to perform self-calibration of the MS, which is followed by refining the mask using this self-calibrated data. Finally, the visibilities are modelled with this new, refined mask. To further improve calibration, we ran the 2GC.py two more times by using the final mask produced as the input mask instead of using the mask made by the FLAG.py script. The resulting calibrated image fits file was then assessed using the Cube Analysis and Rendering Tool for Astronomy (CARTA; Comrie et al. 2018).

4. Results

4.1. L band

In Fig. 1 we present the final L-band image obtained by following our data reduction procedure as laid out in Sect. 3. Cyg X-1 is seen as a compact, bright source with a flux density of ≈13 mJy near the centre of the field and a clear U-shaped arc is identified to the north of Cyg X-1 (Fig. 2, left panel), which is similar to the bow shock identified near Cyg X-1 by G05. There is diffuse emission detected to the south of Cyg X-1 too (W2 in Fig. 1; see discussion of ‘The Whale’), although it is unlikely that this is associated with Cyg X-1 given the scale of the emission compared to the size of the northern bow shock.

Fig. 1. Cyg X-1 field as observed by the MeerKAT telescope in the L band. The orange square marks the location of Cyg X-1, which is seen as a compact source in the middle of the square. The dashed orange arc to the north of Cyg X-1 denotes the outer edge of the bow shock structure. The H II region known as the Tulip Nebula is seen on the eastern side of Cyg X-1. W1 and W2 are two new diffuse emission regions detected in the image to the north-west and south of Cyg X-1, respectively.

Fig. 2. Zoomed-in image of the Cyg X-1 bow shock as observed by MeerKAT (left and middle panel) and by WSRT (G05; right panel). Left panel: Image of Cyg X-1 and the bow shock. The purple arrow represents the direction of the hard state jet of Cyg X-1 as resolved by Stirling et al. (2001). Middle panel: L-band MeerKAT image of the bow shock overlaid with orange contours to show the significant emission region in the bow shock. The contours start at a 2σ level and increase in steps of σ, where σ is the rms. Right panel: Bow shock contours mapped using MeerKAT L-band data overlaid onto the WSRT radio map of the bow shock in G05 to compare the bow shock over a time span of ≈15 years.

To estimate the parameters of the detected bow shock, we made a contour map of the bow shock region as shown in Fig. 2, middle panel. The orange contours in the image start at the 2σ level and increment in steps of σ, where σ is the rms noise in the region around the bow shock. The emission region marked by these contours will be referred to as bow shock region R1 in the rest of the analysis in this paper. The average brightness measured in R1 is 112 ± 24 μJy bm⁻¹, and the integrated flux density of the region (after subtracting compact background emission sources) is 50 ± 2 mJy. The leading edge of the bow shock region is marked as A1 in Fig. 2 (middle panel) and is located at an angular separation of 10.8 arcmin (L_l in the formulations in the following sections) from Cyg X-1, which corresponds to 6.8 pc for the distance of Cyg X-1. The angular separation between the edges of the bow shock region closest to Cyg X-1 (regions marked as A2 and A3) is assumed as the diameter of the radio lobe and is R_l = 7.3 arcmin (4.6 pc). R1 has varying thickness in different parts of the shock region, so we estimated the average thickness by measuring the width at 20 different areas of the bow shock. The average thickness of the bow shock region (ΔR_l) was estimated to be 1.5 ± 0.2 arcmin (0.9 pc). We also overplot R1 contours on the image of the bow shock from G05 to get a comparison of the scales and possible changes in the morphology of the structure (Fig. 2, right panel). This comparison will be discussed in Sect. 5.

4.2. S band

Figure 3 shows the image of our S-band observations’ full field of view. Cyg X-1 is detected as a compact source with a flux density of 14 mJy. The bow shock is detected to the north of Cyg X-1 and a zoom-in into the bow shock is plotted in Fig. 4. W1 and W2 are also detected in the S band; however, W2 is contaminated by side lobes, making it challenging to calculate the brightness of this region. These side lobes can be suppressed by making a uv cut at 500 wavelengths, pointing towards short baseline error that is uncalibrated. The resulting image is reported in Fig. A.1. Since the side lobes do not interfere with the location of the bow shock, we derived the bow shock parameters in the S band before applying the uv cut to maintain optimal sensitivity and resolution. We made a contour map focused on the Cyg X-1 bow shock as seen in Fig. 4, middle panel. The orange contours start at 2σ, and increase by σ to a value of 1.27 × 10⁻³ Jy beam⁻¹, where σ is 10 μJy bm⁻¹ and is the rms. The average brightness of the bow shock is estimated by considering the emission within the faintest contour level around the bow shock as the edge of the bow shock and is 38 ± 10 μJy bm⁻¹. The integrated flux density of the bow shock is 51 ± 1 mJy. The angular separation between the leading edge of the bow shock (B1 in Fig. 4, middle panel) and Cyg X-1 is 10.5 arcmin (L_s in the formulations in the following sections). This separation translates to 6.7 pc for the distance of Cyg X-1 and agrees with L_l. The edges of the bow shock closer to Cyg X-1 are marked as B2 and B3 in Fig. 4, middle panel. The angular separation between these regions is the diameter of the bow shock (R_s) and is measured to be 7 arcmin (4.4 pc). The width of the bow shock (ΔR_s) is calculated in the same way as the L-band observations, where we averaged the width of the bow shock in 20 areas of the bow shock and the mean width is 1.4 ± 0.2 arcmin, which is consistent with ΔR_l.

Fig. 3. Cyg X-1 field as observed by the MeerKAT telescope in the S band. Like Fig. 1, the location of Cyg X-1 is marked with an orange square and the bow shock structure is marked with a dashed orange arc. The Tulip Nebula, W1, and W2 are all also detected in the S band.

Fig. 4. Zoomed-in image of the Cyg X-1 bow shock as observed by MeerKAT in the S band. Left panel: Image of Cyg X-1 and the bow shock. Middle panel: Contours of the S-band emission overlaid on the image, with B1, B2, and B3 marked to show the emission regions used to estimate various bow shock parameters. B4 and B5 mark the location of the smaller arcs within the larger bow shock structure. We also trace these structures with black arcs to draw the reader’s eye to the structure. Right panel: L-band bow shock contours overlaid on the S-band image of the bow shock to compare the morphology of the bow shock in the two frequency bands.

5. Discussion

5.1. Spectral index

Particles accelerated due to interaction with the ISM have a power-law energy distribution at the shock front that is seen as a power-law radio spectrum (S_ν  ∝  ν^α, where α is the spectral index). The electron population in these regions can emit non-thermally, primarily via the synchrotron process or thermally via Bremsstrahlung processes. The electron density, the particle acceleration process and its efficiency, and the temperature of the shock determine the emitting process of the shock, and in turn the spectral index. A flat spectral index indicates thermal emission processes whereas a steep spectral index suggests that synchrotron emission processes contribute to the bow shock emission (Longair 2011). Here, we investigate the integrated spectral index and the spectral index distribution of the bow shock to understand the origin of the emission. We compared the measured flux densities of the bow shock detected in the two frequency bands of the MeerKAT data. The uv plane in the L band has a denser coverage of the shorter spacings as compared to the S band, and thus we forced the same uv-minimum cut in both the L- and S-band image to recover similar angular scales at the two frequencies. We then convolved the S-band image to the resolution of the L-band image (14 arcsec) so we can compare the images with the same resolution and minimise resolving out flux in the S band. We also only map the spectral index of emission that was greater than 2σ, where σ is 15 μJy bm⁻¹ and 20 μJy bm⁻¹ for the L-band and low-resolution S-band images, respectively.

In Fig. 5, we present the spectral index distribution of the Cyg X-1 bow shock (left panel) and the right panel represents the error in the spectral index measurement. We have also plotted the histogram of the spectral index and spectral index error distribution. The spectral index is mostly smooth across the bow shock and the 68% confidence interval of the values is −0.9 ≤ α ≤ 0.4. The spectral index error distribution has a 68% confidence interval of 0.6 ≤ Δα ≤ 1.5. Thus, the error bars on the spectral index map are too large to conclude whether the emission in different regions of the bow shock is flat or negative spectrum as it appears to be consistent with both. Even though the absolute spectral index values have large uncertainties, we can see that there are no noticeable spectral index gradients, pointing to the same emission mechanism for the complete structure. The integrated flux density of the bow shock is measured as 35 ± 3 mJy and 27 ± 3 mJy, leading to an integrated spectral index of −0.4 ± 0.5. It is worth noting that these integrated flux densities are significantly different than the integrated flux density measured before applying the uv-cut and forcing the same beam size in both bands. We find that the integrated spectral index is also inconclusive to determine the origin of the bow shock emission. To determine the origin of the emission from the radio spectral index alone, we need a more accurate spectral index measurement, which will be possible in the future with our Ultra High-Frequency band observation of this region that covers a frequency range 544–1087 MHz using the MeerKAT.

Fig. 5. Spectral index and error distribution map of the bow shock structure between L- and S-band MeerKAT detections. Bottom panels: Distribution in the RA and Dec plane. Top panels: Histogram representations of the distribution of value and errors of the spectral index. The spectral index does not vary noticeably across the bow shock, and the values lie between −0.9 and 0.4, with the error bars consistent with 1.

5.2. Modelling the jet–ISM interaction site

The bright jet–ISM interaction site as seen in the MeerKAT L-band and S-band images provides an opportunity to estimate the time-averaged power being pumped into the region by the jets of Cyg X-1. Due to the lack of models that directly explain BHXB jet and ISM interactions, we followed a relatively simplistic model of extending the AGN jet–IGM interaction expectations to the BHXB case. This follows the model explored by Kaiser & Alexander (1997) and Kaiser et al. (2004) for AGNs and BHXBs, respectively. They use the surface brightness of the jet–ISM interaction site to determine various parameters of the jet and the radio bubbles, extending it to the specific case of a BHXB, GRS 1915+105. A similar, though simplified, model was applied to the bow shock detected around Cyg X-1 by Gallo et al. (2005). In this paper, we extend that model and use it to explain the new L-band and S-band detections of the bow shock. There are a few notable differences between the modelling input parameters in G05 and our work, which have been described in the following section. Here, we present detailed calculations of the ISM density of the surroundings, the expansion velocity of the bow shock, the age of the jet flow and the energy transport rate that has been able to create the bow shock. We have also provided the Jupyter Notebook and thus the Python code that we used to do these calculations (see Sect. 6). We begin by summarising the input parameters and assumptions of the model. We explain the choices when these parameters are mentioned in the relevant sections:

• The timescale of the variability of the BHXB is small as compared to the dynamical change timescale of the large-scale radio lobes.

• The bow shock radio emission is of Bremsstrahlung origin and the assumed temperature range of the emitting region is 10⁴ < T < 3 × 10⁶ K (G05).

• The ISM gas producing the radio emission is completely ionised (ionisation fraction is 1).

• The bow shock is as deep in our line-of-sight as it is wide in the plane of the sky (ΔR).

• The jet opening angle has been assumed to be 0°–30° to obtain generalised estimations.

• We assumed that the bow shock brightness is the same everywhere within the estimated error bars and thus all estimations are homogeneous across the bow shock.

The first step in modelling the jet–ISM interaction site is to determine the shocked and unshocked ISM density. For the case of the Cyg X-1 bow shock, we used the radio emission detected from the bow shock in the L- and S-band MeerKAT observations and the H_α surface brightness of the bow shock (Russell et al. 2007; Sell et al. 2015).

5.2.1. Using radio emissivity

When the radio jets from Cyg X-1 plunge into its environment, it creates an over-pressured radio lobe that expands sideways. The over-pressured lobe allows the jet to propagate efficiently with minimal energy losses until it reaches an over-dense medium. As the end of the jet interacts with this over-dense region, it shocks the ISM and transfers energy to the medium. This results in a shock-compressed ISM at the edge of an over-pressured lobe. Due to limb brightening effects, we detect the edges of the shocked fronts in our radio and optical maps. We present a schematic of our model and a simplification that explains what we observe in our radio observations in Fig. 6. As presented in the schematic, the only signature of the jet–ISM interaction observed in the case of Cyg X-1 is the bow shock structure. We do not detect a distinct end of the jet as it interacts with the ambient medium. We also do not detect the over-pressured lobe. The existence of these signatures is a result of the interplay between the jet power, the environment in which the jet is dissipating energy and the particle content of the plasma. It is not straightforward to disentangle the contribution of these parameters to single out the reason for not detecting hotspots or the over-pressured lobe. The same issue has long been an open question in the extragalactic equivalent of these systems, namely radio galaxies. For instance, the radio jets of Cygnus A are associated with radio bubbles, hotspots and bow shock structure (Carilli & Barthel 1996) whereas no hotspots are detected in the case of Centaurus A (Feain et al. 2011). For the case of Cyg X-1, we speculate that the ISM and the jet-induced cavity are well mixed, such that the jet is not very collimated at the distance where the bow shock structure is detected. This would lead to a more gradual dissipation of kinetic energy than what is needed for creation of a distinct, bright hotspot. This idea is also supported by the detection of a broader, non-uniform bow shock structure instead of a sharp-edged bow shock. The analogy with radio galaxies is furthermore seen when comparing Cyg X-1 with the recently detected bow shock driven by the Galactic system GRS 1915+105. GRS 1915+105 has been associated with a hotspot, an over-pressured bright lobe (Chaty et al. 2001; Rodriguez & Mirabel 1998), and a bow shock structure (Motta et al. 2025), making it distinctly different from the behaviour of jets from Cyg X-1 at large scales.

Fig. 6. Panel A: Simplified schematic of the observations of the jet–ISM interaction site near Cyg X-1 as detected in the MeerKAT image. Panel B: Schematic of the model used to derive the jet parameters of Cyg X-1.

We attempted to measure the integrated radio spectral index and its distribution across the bow shock with the data presented in this paper; however, the measurements are not accurate enough to determine whether the emission can be attributed to Bremsstrahlung processes. Thus, here we relied on the conclusions of Gallo et al. (2005), who suggested that the surface brightness of the bow shock in their WSRT L-band detection and the H_α detection points to the radio emission of the ring-like structure being of Bremsstrahlung origin. The minimum temperature of the bow shock is thus > 10⁴ K, below which ionised hydrogen gas is unable to efficiently emit braking radiation. Gallo et al. (2005) could not find any soft X-ray emission from the bow shock region and thus put an upper limit of 3 × 10⁶ K. For our calculations, we considered the temperature range 10⁴ < T < 3 × 10⁶ K.

The Bremsstrahlung radio emissivity (ϵ_ν) and monochromatic luminosity (L_ν) of an ionised hydrogen region (Longair 1994) is related as

$$\begin{array}{c}\hfill {ϵ}_{\nu }=\frac{L_{\nu }}{V}=C_{\mathrm{radio}}g(\nu ,T)\frac{n_e^2}{\sqrt{T}}exp\left(\frac{h\nu }{k_{b}T}\right),\end{array}$$

where L_ν = F_ν × 4π × d² is the monochromatic luminosity, F_ν is the integrated radio flux density of a flat spectrum source at a central frequency ν, d is the distance to the flat spectrum source, V is the volume of the emitting region, T is the temperature of the shocked ISM, C_radio = 6.8 × 10⁻³⁸ erg s⁻¹ cm⁻³ Hz⁻¹, $g(\nu ,T)\approx \frac{\sqrt{3}}{2\pi }[ln\left(\frac{128{ϵ}_0^2{k}_b^3{T}^3}{m_e{e}^4{\nu }^2{Z}^2}\right)-{\gamma }^2]$ is the Gaunt factor, ϵ₀ = 1, k_b is the Boltzmann constant, m_e is the mass of electron, e is the electric charge, ν is the central observing frequency, Z is the hydrogen atomic number, and n_e is the electron density of ISM.

The bow shock structure appears as an arc or a U-shaped curve with thickness ΔR. Knowledge of the true shape of the bow shock is sensitivity limited and suffers from projection effects of the three-dimensional bow shock as a two-dimensional bow shock structure in the image plane. For this paper, we only use the detected radio emission in our observations to estimate the volume and thus emissivity for the bow shock structure. We summarise these values for both the L-band and S-band detection of the bow shock in Table 1. Rearranging the equation for Bremsstrahlung emissivity, we can derive the electron density as

$$n_e=\sqrt{\frac{L_{\nu }}{Vg(\nu ,T)\sqrt{T}{C}_{\mathrm{radio}}}}·$$

Using the values of the measurements as reported in Table 1, we can estimate the electron density of the shocked ISM for both the L-band and S-band bow shock parameters. Since there are error bars associated with the measurements of the ring thickness, the integrated flux density, distance to Cyg X-1, and a range of temperatures that can emit the radiation, we use a Markov Chain Monte Carlo (MCMC) approach to derive a distribution of the electron density of the shocked ISM. The resulting probability density distribution using the bow shock parameters derived from the L-band detection suggests that the electron density of the shocked ISM is $n_{\mathrm{eL}}(T)={24}_{-8}^{+12}$ cm⁻³, representing the median, 16 and 84 percentiles of the distribution. Similarly, using the bow shock parameters from the S-band detection we estimate the electron density of the shocked ISM to be $n_{\mathrm{eS}}(T)={28}_{-9}^{+14}$ cm⁻³. Such high shocked ISM densities have been observed in regions around two other XRBs; Vela X-1 (van den Eijnden et al. 2022) and the newly reported GRS 1915+105 (Motta et al. 2025).

At high temperatures of > 10⁴ K, almost all the hydrogen will be ionised at such densities and thus the ionisation fraction will be unity. The total number density of the shocked ISM will then be the same as the electron density, n_eL(T) and n_eS(T). This is a notably different assumption as compared to the one used in G05, resulting in a different total number density of the shocked ISM as compared to the one reported in G05. For strong shocks and assuming monatomic gas, shocked ISM is four times denser than the pre-shocked ISM (Landau & Lifshit’s 1991), the pre-shocked ISM number density will be ≈6 cm⁻³ when considering the L-band bow shock parameters and ≈7 cm⁻³ when considering S-band bow shock parameters.

5.2.2. Using the H_α surface brightness

The Cyg X-1 bow shock has been the subject of radio and optical campaigns to understand the energetics responsible for creating the bow shock structure. Russell et al. (2007) reported the detection of the bow shock structure using an H_α filter, which supported a thermal plasma origin for the emission seen in the radio and optical bands. Similar to radio surface brightness, the H_α surface brightness can also be used to determine the electron density of the emitting region. We used the following relation between H_α surface emissivity (j_Hα), the temperature of the emitting region (T_h), and the electron density (n_eh) to estimate the possible range of the electron density of the shocked gas that can explain the detected H_α emission (Osterbrock 1989; Mackey et al. 2013; Gvaramadze et al. 2018; van den Eijnden et al. 2022):

$$\begin{array}{c}\hfill j_{\mathrm{H}\alpha }=2.85\times {10}^{-33}\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-3}{\mathrm{arcsec}}^{-2}{\left(\frac{T}{\mathrm{K}}\right)}^{-0.9}{\left(\frac{n_{\mathrm{eh}}}{{\mathrm{cm}}^{-3}}\right)}^2.\end{array}$$(1)

The surface brightness (S_Hα) that is measured from H_α maps is a result of integrating the surface emissivity along the line of sight in the bow shock and is given by

$$S_{\mathrm{H}\alpha }={\displaystyle \int j_{\mathrm{H}\alpha }}(x)dx=j_{\mathrm{H}\alpha }\mathrm{\Delta }.$$

Rearranging these equations, we can use the surface brightness and temperature range to get an ISM density distribution:

$$n_{\mathrm{eh}}(T)=\sqrt{S_{\mathrm{H}\alpha }\frac{{(T/\text{K})}^{0.9}}{C_{\mathrm{H}\alpha }\mathrm{\Delta }}},$$

where C_Hα = 2.85 × 10⁻³³ erg s⁻¹ cm⁻³ arcsec⁻² and Δ is the width of the bow shock in cm. The background subtracted H_α surface brightness of the Cyg X-1 bow shock is ≈1.3 × 10⁻¹⁶ ergs s⁻¹ cm⁻² arcsec⁻² (Sell et al. 2015). We limit temperature to a conservative range of 5 × 10³–2 × 10⁴ K motivated by the range of temperatures for which Eq. (1) is valid. We note that the temperature range is significantly lower than the consideration for the Bremsstrahlung emitting region to ensure that all hydrogen is not ionised and can still emit H_α. Using the above temperature range and a 15% error budget for S_Hα and Δ, we perform MCMC simulations to obtain a distribution of n_eh. We find that the ISM densities that can explain the detected H_α emission can be constrained to 11 ± 2 cm⁻³, where the error bars are the 16th and 84th percentile of the distribution.

5.2.3. Combining radio and Hα maps

The shocked ISM density determined above using the radio detections in the L and S bands is 2–3 times higher than the ISM density that can explain the emission detected in the H_α band. The possible temperature ranges combined with different ISM densities suggest that the bow shock structure could be a result of multiple shocks, with a contact discontinuity between the radio and H_α emitting region to be a factor of 2–3 (n_eL/n_eh to n_eS/n_eh). This discontinuity factor is moderate given that strong shocks could lead to enhancements by a factor of 4, and even up to 9 as has been shown for the case of Vela X-1 (Gvaramadze et al. 2018). Density and temperature jumps were also needed to explain the X-ray bow shock observed in Cygnus A (Clarke et al. 1997). Another extreme example of a single object showing spatially distinct regions due to the different shock regions is in the case of a high velocity pulsar B1957+20, where the X-ray and H_α emission trace the termination shock and the forward shock, respectively (Stappers et al. 2003). In the case of Cyg X-1, a smoking gun of this would be a clear spatial differentiation between the radio and optical emission regions, which is currently not possible due to the resolution limitations of both the radio and optical images of the bow shock structure. The presence of substructure, background/foreground point-sources also makes it complicated to compare the radio and optical bow shock maps.

The density values of the unshocked ISM (6–7 cm⁻³) are consistent with the complex environment of Cyg X-1 given that there is a massive, bright, star-forming region namely the Tulip Nebula very close to Cyg X-1. Additionally, if the newly discovered extended emission of The Whale turns out to be at a similar distance as that of Cyg X-1, this would be another over-dense region close to Cyg X-1. The unshocked ISM density values are in agreement with high-density ISM surroundings of Vela X-1 (van den Eijnden et al. 2022) and GRS 1915+105 (Motta et al. 2025). Conversely, if the integrated surface brightness of the Cyg X-1 bow shock was at the same level as the rms around the region, we would not have been able to detect it in our observations. Thus, if the bow shock existed but was undetected in our L-band observations, it would imply that the electron density of the shocked ISM was between 1–3 cm⁻³, and the unshocked ISM would be ≈0.5 cm⁻³. This provides an upper limit of the ISM densities of the environments around other XRBs where such jet–ISM interaction sites have not yet been detected in the ThunderKAT/XKAT data.

5.3. Velocity and age of the bow shock

Based on Rankine–Hugoniot conditions, which describe the parameters of the medium on both sides of a shock front, the expansion velocity ($\dot{L}$) of a strong shock (i.e. a very high Mach number) in a mono-atomic gas can be approximated as (Landau & Lifshit’s 1991; Kaiser et al. 2004)

$$\dot{L}\approx \sqrt{\frac{16k_b}{3m_p}T},$$

where m_p is the mass of the proton and T is the temperature of the gas. The velocity of the bow shock can be used to estimate the velocity of the leading edge of the jet that is driving into the ISM and creating the bow shock, since Kaiser & Alexander (1997) showed that these two velocities are roughly similar. Using this expression, we estimate the velocity of the bow shock to be 21 km s${}^{-1}<\dot{L}<364$km s⁻¹ for a temperature range 10⁴ K  < T < 3 × 10⁶ K. This velocity is consistent with results from G05, as the same temperature range was used in both cases.

Since the expansion velocity of a bow shock driven by a jet is nearly equal to the velocity of the jet’s driving end, the bow shock’s age can provide constraints on the age of the jet itself (Kaiser & Alexander 1997; Kaiser et al. 2004). We adopted the methodology and assumptions of Kaiser et al. (2004), who modelled the jet–ISM interaction sites near GRS 1915+105 (also Motta et al. 2025). Both the approaches mentioned above use the basic foundation built in Kaiser & Alexander (1997) to model the formation and evolution of hotspots and bow shocks around FR II type galaxies and were later applied to the features around the BHXB GRS 1915+105. The large-scale structures around FR II galaxies are approximated as self-similar features, a result attributed to an ambient density profile that decreases as a function of the square of the distance to the bow shock (Scheuer 1974). Falle (1991) updated this model by suggesting that the self-similar evolution of a cocoon needed the density profile of the ambient medium to be a function of $\frac{1}{d^2}$ or steeper. Finally, Kaiser & Alexander (1997) were able to show that both the cocoon and the bow shock would evolve showing self-similarity if the ISM density around the source was less steep than $\frac{1}{d^2}$. This helped them conclude that the length of the jet that powers the bow shock grows with time as L_j ∝ t^3/5, assuming a uniform environment of the source. The above models are all approximate since the ambient density is not always as simple as a power law, nor are large-scale structures always due to the jets being self-similar. However, in the absence of a complete physical model, these approximations will provide a good starting model that can be applied to XRBs too. Hence, we assume that in our model (like for AGNs and in the case of Kaiser et al. 2004), the lifetime of the large-scale structure of the jet can be estimated as a ratio of the length of the jet to the velocity of the apex of the bow shock and is given by

$$t=\frac{3L_j}{5{\dot{L}}_j}·$$

Using the range of $\dot{L_j}$ derived in Sect. 5.3 and the distance of the apex of the bow shock from the Cyg X-1 as L_j, we calculate the lifetime of the bow shock of Cyg X-1 to be in the range 0.04–0.3 Myr, where the limits are the 16 and 84 percentiles of the resulting distribution. We assume that the duty cycle of the jets of Cyg X-1 is 90% based on the observation that the radio jets of Cyg X-1 are in a hard state and power collimated jets for 90% of its lifetime (Gallo et al. 2005). Thus, we can extrapolate this assumption in claiming that the lifetime of the large-scale structure and the creation of jets from Cyg X-1 is the same, and hence the age of the jets is estimated to be 0.04–0.3 Myr. The jet that is responsible for the bow shock structure is orders of magnitude younger than the lifetime of Cyg X-1 itself, which is estimated to be between 5–7 Myr (Mirabel & Rodrigues 2003).

5.4. Energy transport rate

It has been shown for FR II radio galaxies that a constant supply of energy is required to sustain the radio cocoons and bow shocks around these AGNs On large scales, it is expected that over the lifetime of the jet, the energy transport rate of the jet is constant. Falle (1991) suggested that the evolution of large-scale structure is self-similar beyond the characteristic length scale, and Kaiser et al. (2004) showed that this was indeed the case for XRBs like GRS 1915+105. The characteristic length is calculated as

$$L_0={\left[\frac{Q_0^2}{{\rho }_0^2{c}^6({\gamma }_j-1)}\right]}^{1/4},$$

where Q₀ is the energy transport rate, ρ₀ is the density of the unshocked ISM and γ_j is the Lorentz factor of the jets. We assume a γ_j of 1.1 for this calculation to maximise the characteristic length estimation. Here, we show that the characteristic scale length of the jet (8 × 10⁻⁵ pc) is much smaller than the large scale jet length that is powering the bow shock in Cyg X-1, which suggests that the modelling used to estimate the time-averaged energy of the jets in GRS 1915+105 applies to Cyg X-1.

Thus, we can estimate the power of the jet (Q_jet) by using our calculations of the density of the ISM (n_eS), the length of the jet (L_j), the velocity of the bow shock ($\dot{L_j}$) and a constant (C₁). C₁ depends on the opening angle and thermodynamical properties of the jet. The opening angle of the low-hard state jet of Cyg X-1 was inferred to be 0.3°–2° using high-resolution imaging capabilities of the Very Long Baseline Array. However, the opening angle of the jet has not been directly measured in the intermediate X-ray state of Cyg X-1 and it is possible that the opening angle changes between different X-ray spectral states of a BHXB. Thus, here, we use a wider uniform distribution prior of the opening angle to be 0°–30° for estimating the energy transport rate. The calculation and variation of C₁ with opening angle is presented in the accompanying Jupyter notebook. The equation used to estimate the energy transport rate is below for reference:

$$Q_{\mathrm{jet}}={\left(\frac{5}{3}\right)}^3\frac{{\rho }_0}{C_1^5}{L}_j^2{\dot{L_j}}^3.$$

We determine that the energy transport rate of the jet that powers the bow shock of Cyg X-1 is of the order of 2 × 10³⁴ ergs s${}^{-1}<Q_{\mathrm{jet}}^{\mathrm{b}}<1\times {10}^{37}$ ergs s⁻¹ for the opening angles 0°–30°. For the more stringent range of opening angles (0.3°–2°), the energy transport rate is in the range 2 × 10³¹ ergs s${}^{-1}<Q_{\mathrm{jet}}^{\mathrm{a}}<1\times {10}^{35}$ ergs s⁻¹. The constraints on the energy transport rate of the jet averaged over the jet lifetime for Cyg X-1 is measured to be lower than that estimated for GRS 1915+105. This is expected given that the ISM in which Cyg X-1 is embedded is less dense (6–7 particles cm⁻³) than that of the region surrounding GRS 1915+105 (100–160 particles cm⁻³; Motta et al. 2025), and the scale length of the jet is larger for GRS 1915+105. To use this hierarchy of jet power to represent the true jet powers of the two XRBs, we need an estimate of the age of the XRBs. Cyg X-1 is a high mass XRB whose accretion disc is fed by the strong winds from its massive companion (Čechura et al. 2015), with an estimated age of 10⁶ yr (Mirabel & Rodrigues 2003). On the other hand, in the case of GRS 1915+105, the accretion is sustained by Roche lobe overflow, with the system being reported to have near-Eddington accretion for a large fraction of its active phases (Done et al. 2004), and it is also older than Cyg X-1 (Podsiadlowski et al. 2003). The estimations for jet power for both Cyg X-1 and GRS 1915+105 are lower than the lower limit for the power distribution of jetted AGNs (10⁴³ ergs s⁻¹; Fan & Wu 2019).

We can also estimate the pressure inside this bow shock and it is given by

$$p=0.0675\frac{C_1^{10/3}}{R^2}{\left(\frac{{\rho }_0{Q}_0^2}{L_j^4}\right)}^{1/3}.$$

Thus, the pressure is estimated as 3 × 10⁻¹³ ergs cm⁻³ < p < 2 × 10⁻¹¹ ergs cm⁻³. Hence, the Cyg X-1 bubble has a lower pressure as compared to the GRS 1915+105 lobe (1.9 × 10⁻¹¹–3 × 10⁻⁹ erg cm⁻³).

Most of the above analysis was based on the bow shock being a thermal emitter. However, some regions of the bow shock are also consistent with −1 spectral index, although the error bars are large in those regions. So, for completeness, we also estimated the minimum magnetic field associated with non-thermal emission of the bow shock. This minimum magnetic field occurs at close to the equipartition of energy between the electron population and the magnetic fields present in the emitting region. The equipartition magnetic field is given by

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

where f is a filling factor (assumed 1), c₁₂ is a constant that relates the energy stored in the electrons to the magnetic field, L is the total luminosity of the emitting region across the frequency band and V is the volume of the emitting region. For an assumed spectral index of −0.9 informed by the uncorrected values of our spectral index map, the minimum magnetic field is > 1 × 10⁻⁶ G for L-band and > 8 × 10⁻⁷ G for S-band luminosity. These values are lower than those estimated for the bow shock around GRS 1915+105 (Motta et al. 2025).

5.5. Bow shock morphology

The velocity of bow shock expansion can also be measured directly in the plane of the sky by measuring the change in the position of the bow shock over time. However, such measurements of faint diffuse emission like that of the Cyg X-1 bow shock would require either an extremely long time baseline or a very fast bow shock. As a test, we compared the bow shock location and morphology in G05 and the bow shock location as measured by our observations (see Fig. 2, right panel). The bow shock observations although separated by 17 years, are almost aligned with each other as is seen by the contours of the L-band MeerKAT bow shock overlaid on the G05 image of the bow shock with WSRT. There is some apparent change in the morphology of the bow shock. The leading edge of the bow shock is seen connected leading to a complete U-shaped arc in G05. On the other hand, the leading edge in the most recent L-band observations with MeerKAT shows that the edge is as faint as the noise in the image with some connecting structure visible, which is at a 2σ level. However, the leading edge of the bow shock that is barely detected in the L band is significantly detected in the S band. A longer observation of the bow shock in the L band will help improve the uv coverage and also increase the sensitivity, aiding in assessing whether the leading edge of the bow shock persists today in the L band.

As a direct consequence of the S-band MeerKAT image being more sensitive than the L-band image, we can identify new features in the bow shock, like the two smaller arcs (marked as B4 and B5 in Fig. 4, middle panel, and seen clearly in the Fig. 4, left panel) within the larger bow shock. Hints of these smaller structures are also seen in the L-band MeerKAT image, but are confusion-limited. We have overlaid the contours of the L-band bow shock on the S-band image in Fig. 4 (right panel). These smaller arcs were seen as a solid extended structure between the arms of the U-shaped arc of the bow shock near the leading edge in G05. Here, we present three tentative interpretations of the detection of these smaller arcs. The first one is that the smaller arcs are the result of different transient jet-ejection episodes that might have occurred over the lifetime of Cyg X-1. This would indicate that some jets have been launched misaligned from the usual, perpendicular jet axis to the accretion disc. These misaligned jets end up sweeping different, smaller regions as compared to the more stable, compact jets that are responsible for the larger U-shaped structure. A second possibility is that depending on the relative location of the black hole and the massive companion star, the stellar winds can bend the extended jet even when it is in the low-hard state creating the smaller arcs (Bosch-Ramon & Barkov 2016), and at larger distances from these smaller arcs, the shocked front could be sweeping up gas and creating an even larger interference shock front pattern. The third possibility is that the smaller arcs are seen due to projection effects, whereby small under (or over) densities in the local ISM could enable the shocks to break out of the overall structure. The two-dimensional image shown in Fig. 4 is a three-dimensional structure, so these smaller arcs may be localised in regions that are seen in projection on the overall bow shock structure. However, we note that the scenarios presented above are currently only hypotheses as the smaller arcs are too faint to be able to do detailed modelling.

5.6. The Whale

Other than the main subject of these observations, the Cyg X-1 bow shock, there is another noticeable large-scale structure in Figs. 1 and 3. The combined emission of W1 and W2 appears to be shaped like a whale and hence we have named the emission structure ‘The Whale’. We added annotations in Figs. 1 and 3 to showcase the structure. Currently, there are no distance estimates for these emission regions; thus, we only present a qualitative discussion of these large-scale structures, as multi-wavelength follow-up and deeper radio observations of this region are needed to understand the nature and driving force of this emission. The Whale could be a single large diffuse structure, or two independent regions that appear to be connected here due to the limited dynamic range of the image. The tail of The Whale (W1) could be a shock interface or a H II region and has an average brightness of 109 μJy bm⁻¹ in L band and 40 μJy bm⁻¹ in the S band. We report a bright (6 mJy in L band and 11 mJy in the S band), compact object at a position of RA ${19}^{\mathrm{h}}{57}^{\mathrm{m}}35\stackrel{\text{s}}{.}687$ and Dec ${35}^{\mathrm{d}}{23}^{\mathrm{m}}38\stackrel{\text{s}}{.}374$ and at an angular distance of 3 arcmin from the apex of the arc-like structure, W1. We used Vizier² to check if this source had any optical counterparts, and there were three Gaia matches within the 7 arcsec radius of this radio source. Two of these Gaia sources have negative parallax, whereas one had a parallax of 0.3 ± 0.2, which hints at its possible Galactic nature. There are three more radio sources near W1 on the same side of this radio source, but none of them had compelling Gaia counterparts. We will conduct a detailed analysis of the possible association and interaction site W1 in a follow-up study.

6. Conclusions

Jet–ISM interaction sites are a promising laboratory for estimating the energetics of the jet powering the site. In this study we extrapolated models developed for AGNs to the BHXB Cyg X-1. We observed the Cyg X-1 region using the MeerKAT telescope in two frequency bands, the L band (central frequency 1.28 GHz) and the S band (central frequency 2.64 GHz). We detected the jet-powered bow shock as a U-shaped arc in both frequencies, with an average brightness of 110 ± 24 μJy bm⁻¹ and 38 ± 13 μJy bm⁻¹ in the L and S bands, respectively. We were able to resolve smaller arcs and finer bow shock structure in the S-band image, which shows that the morphology of the bow shock is more complicated than a single U-shaped edge brightened bow shock structure. We report the spectral index map of the bow shock and find that the distribution of spectral indices is uniform throughout the structure. In the future we will conduct observations in other frequency bands to accurately measure the absolute spectral index.

We combined the observed bow shock parameters from our L- and S-band observations to estimate the density of unshocked ISM (6–7 cm⁻³) within a temperature range (1 × 10⁴–3 × 10⁶ K) that could explain the emissivity in the radio bands. We modelled the bow shock around Cyg X-1 and determined the age of the Cyg X-1 jets to be 0.04–0.3 Myr and the power of the jet feeding the bow shock structure to be 2 × 10³¹ ergs s${}^{-1}<Q_{\mathrm{jet}}^{\mathrm{a}}<1.0\times {10}^{35}$ ergs s⁻¹ for the case of opening angles of 0.3°–2.0°. This well-constrained jet power is lower than the earlier estimate for this bow shock and also lower than the only other BHXB for which this has been determined, GRS 1915+105. The success of these observations suggests that a systematic search for bow shocks around XRBs is a great way to advance studies of black hole jet powers, outputs, and jet–ISM interaction and densities, and to tie inflows and outflows of accreting black hole systems. The future generation of radio telescopes like SKA, SKA-mid, and ngVLA will be able to make drastic advancements as they will be able to probe bow shocks even in lower-ISM-density environments.

Data availability

The datasets used in this paper are available in the public archive of SARAO at https://archive.sarao.ac.za/. The code used to conduct the analysis, create plots and maps is available here https://github.com/pikkyatri/CygX-1-bow-shock.

Acknowledgments

PA is supported by the WISE fellowship program, which is financed by NWO. SEM acknowledges support from the INAF Fundamental Research Grant (2022) EJECTA. JvdE acknowledges a Warwick Astrophysics prize post-doctoral fellowship made possible thanks to a generous philanthropic donation, and was supported by funding from the European Union’s Horizon Europe research and innovation programme under the Marie Skłodowska-Curie grant agreement No 101148693 (MeerSHOCKS) for part of this work. ASTRON, the Netherlands Institute for Radio Astronomy, is an institute of the Dutch Research Council (De Nederlandse Organisatie voor Wetenschappelijk Onderzoek, NWO). 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). The authors are grateful to Elena Gallo for providing the archival radio map of the Cyg X-1 bow shock. The authors would like to thank Fraser Cowie, David Russell, Payaswini Saikia for their inputs on the manuscript of the paper, and Andrew Hughes and Joe Bright for scheduling the MeerKAT observations.

References

Appendix A: Additional figure

Fig. A.1. Cyg X-1 field as observed by the MeerKAT telescope in the S band after applying a minimum uv cut of 500 wavelengths. A black circle marks the location of Cyg X-1 and the bow shock. The locations of W1 and W2 as detected in the L-band image are marked in green. The W2 has become fainter, indicating that this emission was being picked up by the short baselines that have been flagged out here but are present in the L-band image.

All Tables

All Figures

Fig. 1. Cyg X-1 field as observed by the MeerKAT telescope in the L band. The orange square marks the location of Cyg X-1, which is seen as a compact source in the middle of the square. The dashed orange arc to the north of Cyg X-1 denotes the outer edge of the bow shock structure. The H II region known as the Tulip Nebula is seen on the eastern side of Cyg X-1. W1 and W2 are two new diffuse emission regions detected in the image to the north-west and south of Cyg X-1, respectively.

In the text

Fig. 2. Zoomed-in image of the Cyg X-1 bow shock as observed by MeerKAT (left and middle panel) and by WSRT (G05; right panel). Left panel: Image of Cyg X-1 and the bow shock. The purple arrow represents the direction of the hard state jet of Cyg X-1 as resolved by Stirling et al. (2001). Middle panel: L-band MeerKAT image of the bow shock overlaid with orange contours to show the significant emission region in the bow shock. The contours start at a 2σ level and increase in steps of σ, where σ is the rms. Right panel: Bow shock contours mapped using MeerKAT L-band data overlaid onto the WSRT radio map of the bow shock in G05 to compare the bow shock over a time span of ≈15 years.

In the text

	Fig. 3. Cyg X-1 field as observed by the MeerKAT telescope in the S band. Like Fig. 1, the location of Cyg X-1 is marked with an orange square and the bow shock structure is marked with a dashed orange arc. The Tulip Nebula, W1, and W2 are all also detected in the S band.
In the text

Fig. 4. Zoomed-in image of the Cyg X-1 bow shock as observed by MeerKAT in the S band. Left panel: Image of Cyg X-1 and the bow shock. Middle panel: Contours of the S-band emission overlaid on the image, with B1, B2, and B3 marked to show the emission regions used to estimate various bow shock parameters. B4 and B5 mark the location of the smaller arcs within the larger bow shock structure. We also trace these structures with black arcs to draw the reader’s eye to the structure. Right panel: L-band bow shock contours overlaid on the S-band image of the bow shock to compare the morphology of the bow shock in the two frequency bands.

In the text

Fig. 5. Spectral index and error distribution map of the bow shock structure between L- and S-band MeerKAT detections. Bottom panels: Distribution in the RA and Dec plane. Top panels: Histogram representations of the distribution of value and errors of the spectral index. The spectral index does not vary noticeably across the bow shock, and the values lie between −0.9 and 0.4, with the error bars consistent with 1.

In the text

	Fig. 6. Panel A: Simplified schematic of the observations of the jet–ISM interaction site near Cyg X-1 as detected in the MeerKAT image. Panel B: Schematic of the model used to derive the jet parameters of Cyg X-1.
In the text

Fig. A.1. Cyg X-1 field as observed by the MeerKAT telescope in the S band after applying a minimum uv cut of 500 wavelengths. A black circle marks the location of Cyg X-1 and the bow shock. The locations of W1 and W2 as detected in the L-band image are marked in green. The W2 has become fainter, indicating that this emission was being picked up by the short baselines that have been flagged out here but are present in the L-band image.

In the text