MeerKAT discovery of radio emission from the Vela X-1 bow shock Free

ABSTRACT

1 INTRODUCTION

Astrophysical bow shocks are extended shock nebulas observed around a wide range of different objects. For instance, such structures have been observed in pulsars, indicating shocks between the pulsar wind and the interstellar medium (ISM; e.g. Cordes, Romani & Lundgren 1993; Gaensler et al. 2000; Stappers et al. 2003), interacting binary systems (e.g. the fast-moving low-mass X-ray binary SAX J1712.6-3739; Wiersema et al. 2009; or the nova-like cataclysmic variable, V341 Ara; Castro Segura et al. 2021). Bow shocks can be created by outflows from the object in question, such as the jet launched by the X-ray binary Cyg X-1 (Russell et al. 2007). Commonly, however, bow shocks are created when the object itself moves through the ISM at a velocity exceeding the local sound speed. The latter is the case for massive runaway stars, where their strong stellar wind sweeps up interstellar material, creating a forward shock moving at the star’s velocity through the ISM, and a reverse shock at the (higher) stellar wind velocity. The arc-like morphology of the bow shock is determined both by stellar properties, such as mass-loss rate, terminal wind velocity, and stellar bulk motion, as well as the local ISM properties (Comeron & Kaper 1998; Meyer et al. 2016).

Stellar bow shocks are most often detected via shock-heated dust emission at infrared (IR) wavelengths: after the first systematic IR searches in IRAS data (see e.g. Noriega-Crespo, van Buren & Dgani 1997), subsequent catalogues using WISE (the Extensive stellar BOw Shock Survey, or E-BOSS; Peri et al. 2012; Peri, Benaglia & Isequilla 2015) and Spitzer (Kobulnicky et al. 2016) have built up a sample of more than 700 IR stellar bow shocks. Optical line emission, such as H α (e.g. Kaper et al. 1997; Gvaramadze et al. 2018) or [O iii] λ5007 (Gull & Sofia 1979), has also been detected in a much smaller subset of bow shocks; a complicating factor in detecting H α line emission may be that it can be screened by the presence of H ii regions surrounding the runaway star (Brown & Bomans 2005; Meyer et al. 2016).

Bow shock emission is expected in bands beyond the optical and IR. In the shocks between the stellar wind and the ISM, charged particles can be accelerated to relativistic energies via the Diffuse Shock Acceleration mechanism (e.g. Bell 1978a, b; Drury 1983; Matthews, Bell & Blundell 2020, and references therein for reviews). The resulting relativistic electron population can emit non-thermally across the entire electromagnetic spectrum (del Valle & Romero 2012; del Valle & Pohl 2018; del Palacio et al. 2018), primarily via synchrotron processes (expected to dominate at radio wavelengths) and inverse Compton scattering of IR and stellar photons (dominant at and/or above X-ray wavelengths). In addition, thermal free–free emission (bremsstrahlung), from the same electron population that causes the optical line emission, is expected at long wavelengths. The relative contribution of these non-thermal and thermal processes depends on the energetics of the stellar winds, the electron density and temperature in the shock, as well as the shock acceleration efficiency of electrons.

Despite this set of emission processes, bow shock detections outside the optical and IR bands remain exceptionally rare around massive runaway stars. At radio wavelengths, only a single such bow shock has been detected: BD+43^o3654 was detected with the Karl G. Jansky Very Large Array (VLA; Benaglia et al. 2010, 2021) and the Giant Metrewave Radio Telescope (GMRT; Brookes 2016) with a spatially variable and generally steep spectrum (defined as α > 0 where S_ν ∝ ν^−α), which has been interpreted as non-thermal synchrotron emission (Benaglia et al. 2010). Despite further searches in other IR-identified bow shocks, no other radio counterparts have yet been identified in surveys or pointed observations (Rangelov et al. 2019; Benaglia et al. 2021, the latter based on private communication with C. Peri). Similarly, the expected non-thermal X-ray or γ-ray emission has not been unambiguously detected, despite several searches (Schulz et al. 2014; Toalá et al. 2016; De Becker et al. 2017; Toalá, Oskinova & Ignace 2017; H. E. S. S. Collaboration 2018).¹ The best candidates of possible high-energy bow shock counterparts are reported by Sánchez-Ayaso et al. (2018), who associate two unidentified Fermi sources with known bow shocks.

In this work, we present the MeerKAT radio detection of the Vela X-1 bow shock. Vela X-1 is a high-mass X-ray binary (HMXB) system, consisting of an accreting neutron star and the supergiant donor HD 77581 separated at ∼1.7 donor-star radii (53.4 |$\rm R_{\odot }$|⁠) in a tight ∼8.96-d orbit. Vela X-1 is a runaway system, travelling at a bulk velocity of ∼54.3 ± 0.6 km s⁻¹ (Gvaramadze et al. 2018, see Section 4.2). The supergiant donor star launches a strong stellar wind, whose mass-loss rate and terminal velocity have been estimated in numerous earlier studies to lie in the range ∼5 × 10⁻⁷ to 5 × 10⁻⁶|$\rm M_{\odot }$| per year and ∼500 to 1700 km s⁻¹, respectively (see e.g. Kretschmar et al. 2021, and references therein for an in-depth review of this HMXB). Similarly, many literature estimates of its distance exist: here, we adopt the Gaia DR2 distance of |$1.99_{-0.11}^{+0.13}$| kpc derived by Kretschmar et al. (2021) and refer to that work for further discussion.

The Vela X-1 stellar wind creates a bow shock, which was first discovered in narrow-band H α imaging by Kaper et al. (1997). To date, Vela X-1 remains one of only two HMXBs with a bow shock (with 4U 1907+09; Gvaramadze et al. 2011). The bow shock of Vela X-1 was later detected at IR wavelengths as well (Peri et al. 2015; Maíz Apellániz et al. 2018), and its morphology suggests that Vela X-1 may have originated in the Vel OB1 association several Myr ago (Kaper et al. 1997). Using the SuperCOSMOS H α Survey,² Gvaramadze et al. (2018) reveal the presence of large-scale diffuse structure in the wake of the binary and its bow shock; their detailed simulations indicate that Vela X-1 interacts with a local ISM overdensity, where an ISM density of approximately three times that of the ambient density can reproduce the observed large-scale morphology and H α surface brightness.

Here, we report the serendipitous detection of radio emission from the Vela X-1 bow shock in MeerKAT observations targeting the HMXB itself. In Section 2, we present the MeerKAT data and analysis methods, while the results are shown in Section 3. Then, in Section 4, we perform initial analytical calculations to assess the origin and nature of the radio emission, comparing a non-thermal and thermal scenario.

2 OBSERVATIONS AND DATA ANALYSIS

We observed the field around the HMXB Vela X-1 as part of the ThunderKAT Large Survey Project with MeerKAT (Fender et al. 2016), which performs radio observations of active, Southern X-ray binaries, cataclysmic variables, supernovae, and gamma-ray bursts. Two of the three MeerKAT observations included in this work were part of a larger coordinated multiwavelength campaign of the binary; the results of that study, including the MeerKAT results on Vela X-1 itself, will be reported elsewhere (Van den Eijnden et al., in preparation). In this paper, we analyse only the Stokes I data obtained in these observations.

Vela X-1 (J2000 09^h02^m06|${_{.}^{\rm s}}$|86 −40^○33′16|${_{.}^{\prime\prime}}$|9) was observed with the MeerKAT telescope on 2020 September 25, 2020 September 27, and 2020 October 11 (capture block IDs 1600995961, 1601168939, and 1602387062 respectively).The telescope’s L band (856–1712 MHz; hereafter reported at the central frequency of 1.3 GHz) receivers were used, with the correlator configured to deliver 32 768 channels and 8 s integration time per visibility point. For the three observations, there were 59, 61, and 60 antennas used in the array.

Vela X-1 was observed for 30 min in each of the three runs, book-ended by two 2-min scans of the nearby secondary calibrator source J0825−5010. The standard primary calibrator source J0408−6545 was also observed for 5 min at the start of each run.

For each of the three observations, reference calibration was performed using the casa package (McMullin et al. 2007). Bandpass, delay and flux-scale corrections were derived from the scans of the primary, and complex gain and delay corrections obtained from the scans of the secondary, following all flagging of the data. These corrections were applied to the target data, which was then split and flagged using the tricolour³ package. The target data were imaged using wsclean (Offringa et al. 2014), reimaged following the construction of a deconvolution mask, and then self-calibrated using the cubical package (Kenyon et al. 2018) to solve for phase and delay corrections for every 32 s of data.

The self-calibrated data for all three observations were then jointly imaged using ddfacet (Tasse et al. 2018). Direction-dependent gain corrections were derived using killms (Smirnov & Tasse 2015), with the sky partitioned into 15 directions governed by the distribution of bright, compact sources in the field. The data were then re-imaged using ddfacet, applying the directional gain corrections in the process. The final wide-field image presented in Fig. 1 has been primary beam corrected using the katbeam package,⁴ blanking the map beyond the nominal 30 per cent level. Further details on the data reduction process together with the relevant calibration and imaging parameters can be found online⁵ (Heywood 2020).

Figure 1.

Radio continuum emission of the full MeerKAT field of view of Vela X-1, shown by the cross, and its surroundings, created by combining three observing runs. The combined exposure time is 90 min, yielding a 40 |$\mu$|Jy RMS sensitivity. The observations were performed at L-band (1.3 GHz), and reveal both the Vela X-1 bow shock as well as several other large-scale diffuse structures in radio for the first time. The beam size is shown in the bottom left of the image.

Open in new tabDownload slide

3 RESULTS

3.1 MeerKAT observations

In Fig. 1, we show the combined, full-field 1.3-GHz MeerKAT image created from the three observations, extending ∼0.72 degree out from the central pointing towards the HMXB. Multiple large-scale, diffuse structures are visible, including the bow shock immediately above the cross indicating Vela X-1. Similarly, an extended, diagonal structure is observed, which Vela X-1 appears to have crossed through. Both structures, as well as the other extended radio sources in the image, were already known from H α images of this field (i.e. Gvaramadze et al. 2018, using SuperCOSMOS data). In Fig. 2, we show an inset containing the central 0.5 × 0.5 degree square of the field, providing a more detailed view of the bow shock region.

Figure 2.

Inset of the Vela X-1 field shown in Fig. 1, extending 0.25 degree out from the HMXB (indicated by the cross) in RA and Dec. The 12 arcsec circular beam size is shown in the bottom left. The bow shock source region and background region are shown as the black lined and dashed regions; the blue dashed region indicates the comparison overdense region discussed in Section 4.2. The arrow shows the direction of motion of Vela X-1 relative to its surroundings.

Open in new tabDownload slide

Furthermore, in Fig. 3, we overlay radio contours on the SuperCOSMOS image, indicating the remarkable similarity between the observed morphologies in both H α and L-band imaging. This morphological similarity, combined with their relatively similar order-of-magnitude flux densities, suggests a possible connection between the emission mechanisms in the bow shock and the surrounding structures. Such a connection could be explained via thermal processes that are combined with shock-enhanced electron densities in the bow shock, which we will investigate in detail in the next section. While high-spatial-resolution H α images reveal filamentary structures in the bow shock (e.g. Kaper et al. 1997), our L-band MeerKAT image does not reveal similar features. However, this can be attributed to the larger, 12 arcsec beam size of the MeerKAT observations, which is similar to the width of the H α filaments and larger than the separation between them.

Figure 3.

The SuperCOSMOS H α image, overlaid with the contours from the MeerKAT L-band image. The radio contours are shown at levels of 50, 100, and 200 |$\mu$|Jy, calculated after smoothing the image by a factor 4. These levels are chosen to highlight the diffuse structures in the radio image. The large-scale H α and radio structures appear to trace each other closely throughout the field.

Open in new tabDownload slide

A more detailed comparison of the bow shock profile is shown in Fig. 4, where we show the radio flux density from MeerKAT, H α brightness from the SuperCOSMOS Survey, and IR brightness in the W3 band of WISE⁶ along three lines from Vela X-1, angled at 0, 30, and 60 degrees West of North. These comparisons show how all three bands, to their resolution, show compatible bow shock profiles. The H α band, dominated by the stellar PSF at small distances, shows a single peaked profile at the bow shock apex, but a double peaked, filamentary profile for the other directions (note that the third peak in the bottom panel, at ∼1.6–1.7 arcmin, is due to a star). The lower resolution of the IR and radio images (the latter shown by the black bar in the bottom panel), smear out any sub structure that may be present, leaving a broad and smoothly peaked profile.

Figure 4.

The radial flux density or brightness profiles at three directions from Vela X-1, namely 0, 30, and 60 degrees West of North. The black dashed line shows the radio profile, while the red line and dotted line show the logarithmic H α and IR profiles, with arbitrary units, respectively. The black bar in the top right of the bottom panel, shows the beam size of the MeerKAT observation. The difference in resolution between bands is clearly visible.

Open in new tabDownload slide

To measure the radio flux densities of the bow shock, we construct a custom region (black polygon) in ds9 (Joye & Mandel 2003), as shown in Fig. 2. This region was constructed manually, as the combination of increasing flux along shock in the West-ward direction, with the surrounding overdense region on its West-end, prevents using a single flux contour level to define the entire shock. For angles West of North along the shock, the region’s edge is defined where the flux density drops below approximately 100 |$\mu$|Jy, extending West-wards up to point where the H α shock morphology can be identified. On the fainter East side of the shock, un-surrounded by diffuse emission, the shock edge is defined where it drops approximately below the RMS defined below. We discuss the effects of this region definition in Section 4.

The peak radio flux density of the entire bow shock is |$S_{\nu , \rm max} = 270 \pm 40$||$\mu$|Jy beam⁻¹, where the error equals the RMS in the nearby region (dashed circle in Fig. 2) devoid of point sources and strong diffuse emission. We also apply the radioflux⁷ script to measure the integrated radio flux across the entire bow shock, finding |$S_{\nu , \rm total} = 5.27 \pm 0.07$| mJy, corresponding to a mean flux density of |$S_{\nu , \rm mean} = 130 \pm 40$||$\mu$|Jy beam⁻¹. For our later modelling (Section 4), we also note that the projected bow shock area assumed in these calculations is 10 100 arcsec². While the bow shock structure can also be identified in the sub-bands that are not strongly affected by RFI, we do not calculate an in-band spectral index; even when averaging across the entire shock, the S/N is not sufficient for a reliable measurement.⁸

3.2 Existing ASKAP and Chandra data

To further constrain the spectral shape, we therefore instead opt to use observations with different observatories, at different frequencies. First, we consulted the data release of the Rapid ASKAP Continuum Survey (RACS; McConnell et al. 2020), which provides complete sky coverage at declinations below +41^o at a spatial resolution similar to MeerKAT (∼ 15 arcsec). All images in the first data release are taken at a central frequency of 887.5 MHz with a 288 MHz bandwidth, just below (but overlapping with) the MeerKAT L-band. The bow shock location is covered by three pointings – 0911-37A, 0840-37A, and 0855-43A – with median RMS sensitivities of 285, 292, and 222 |$\mu$|Jy. No hints of the bow shock are visible in any of the three images. Combining the latter non-detection with the mean MeerKAT L-band flux density of 130 ± 40 |$\mu$|Jy beam⁻¹, implies an upper limit on the spectral index of α ≥ 1.4 (where S_ν ∝ ν^−α) – which is not constraining for any realistic emission process in the bow shock. Furthermore, the different array configurations of MeerKAT and ASKAP, with the former’s core-heavy distribution especially suitable for the detection of extended emission, further complicates the comparison between the two radio observatories.

If the observed radio emission originates from non-thermal synchrotron processes, it indicates the presence of a population of relativistic electrons accelerated in the shock. As discussed in more detail in Section 4, such an electron population can generate X-ray emission either as the high-energy part of its synchrotron spectrum, or inverse-Compton scattering of dust and stellar photons. Whether the former contributes in the X-ray band depends on physically interesting parameters such as the maximum electron energy and the magnetic field. Hence, we also investigated the X-ray properties of the bow shock.

As Vela X-1 is a persistently-accreting HMXB, its X-ray luminosity (typically ∼10³⁶ erg s⁻¹) will overpower any diffuse emission from the bow shock. Given its X-ray luminosity and relatively small (less than 1 arcmin) angular distance between the binary and the shock, the point spread function of Vela X-1 overlaps with the shock for all X-ray observatories. However, Vela X-1 is viewed almost edge on and consistently shows eclipses between orbital phases ∼0.9 and ∼0.1 (Falanga et al. 2015), when the Chandra count rate drops by two orders of magnitude and most of the X-ray flux is concentrated in a handful of narrow emission lines (Watanabe et al. 2006). Chandra observed the field of Vela X-1 twice during eclipses (ObsIDs 102 and 1926, ∼28 and ∼83 ks exposures, respectively), both times with the High-Energy Transmission Grating (HETG) employed. While this instrumental set-up is optimized for high-resolution X-ray spectroscopy, the observation provides an image of the observed field from the 0^th order data collected on chip S3. Therefore, we downloaded both observations, reprocessed the data using the task chandra_repro in ciao version 4.13 (Fruscione et al. 2006),⁹ and finally generated an image by using dmcopy, selecting only photons between 0.5–10 keV and screening cosmic ray-induced events below 2 keV.

In Fig. 5, we show the resulting 0^th-order image of ObsID 1926 on the left, alongside ObsID 1928, which was taken outside eclipse, for comparison on the right. During eclipse, despite the faintness of Vela X-1, the HMXB is detected as the point source, and faint gratings arms are visible in their characteristic X-shaped pattern. The red region in both panels shows the MeerKAT L-band bow shock region, adjusted slightly to avoid the arms. While the binary PSF only leaks into the shock region outside of eclipse, no evidence for an X-ray bow shock can be seen in the eclipse observation. For the two eclipse observations 102 and 1926, we measure background count rates in the full shock region (2.72 ± 0.05) × 10⁻² and (3.2 ± 0.1) × 10⁻² cts s⁻¹. These rates are slightly lower than the approximate background rate for the entire chip, scaled by the surface of the region, confirming the above conclusion that no significant bow shock X-ray emission is detected. To convert these count rates to a flux upper limit, we use the S3 background spectrum,¹⁰ finding that the former count rate corresponds to a flux of approximately 7 × 10⁻⁹ Jy at 1 keV. Therefore, we set a 3σ upper limit to the integrated bow shock X-ray flux at 1 keV of 2 × 10⁻⁸ Jy.

Figure 5.

The 0^th order S3 chip data taken by Chandra during eclipse (ObsID 1926; left) and out of eclipse (ObsID 1928, right). The red region indicates the same bow shock region as in Fig. 2, adjusted slightly to avoid the X-shaped grating arms. Despite the lack of counts from the HMXB leaking into the shock region during eclipse, no bow shock counterpart can be identified.

Open in new tabDownload slide

4 DISCUSSION

In this work, we report the 1.3-GHz radio discovery of the bow shock of the HMXB Vela X-1 with the MeerKAT telescope. We do not detect a counterpart at either lower radio frequencies (using RACS) or X-ray energies (Chandra). Therefore, we lack the most direct observable that can help distinguish between a thermal and non-thermal origin of the radio emission: a spectral index measurement.¹¹ Instead, in this discussion, we will perform simple, analytical calculations to compare these two scenarios, and investigate their compatibility with the energetics of the system and the constraints from other wavebands. To increase the reproducibility of our results, these calculations can be repeated in the Jupyter notebook accompanying this paper (see Data Availability). We stress, however, that these two processes are not mutually exclusive. As we discuss briefly at the end of Section 4.3, both processes are likely to co-exist and contribute to the total radio emission.

4.1 General considerations

Before evaluating the two scenarios mentioned above, we briefly discuss some general considerations for these calculations. First, we should consider whether the system is in a steady state, given its prior interactions with the overdense region (Gvaramadze et al. 2018). Given the separation of ∼5 arcmin between the current position of Vela X-1 and this overdensity, tracing back Vela X-1’s direction of motion (shown by the arrow in Fig. 2), and its space velocity of ∼54.3 km s⁻¹ and 1.99 kpc distance, this interaction likely took place ∼5.5 × 10⁴ yr ago. Simulations of bow shock formation in, for instance, Betelgeuse show how a steady state is achieved at similar time-scales, even though no bow shock is present at the start of the simulations (Mohamed, Mackey & Langer 2012). A similar conclusion follows from comparing the standoff distance of the bow shock to its distance from Vela X-1 at a 90 degree angle from the apex: in radio, this ratio R(0)/R(90) = 0.50 ± 0.12, consistent with the expected value of 0.57 for the shape of a steady state bow shock (Wilkin 1996). Based on these two lines of reasoning, we will thus assume the bow shock is in a steady state.

For our calculations, we will assume that the bow shock is homogeneous throughout, which enables us to perform simple analytical calculations and estimates. In reality, the radio image, as well as the filamentary H α image, show that this assumption does not hold. However, we prioritize performing analytic calculations at this stage, and will leave the spatially resolved analyses to future works using simulations beyond the scope of this work. In addition, the radio flux density varies slightly on the scale of the MeerKAT beam: the peak and average flux density differ by a factor ∼2. We expect that our physical estimates will trace the variations across the shock to order of magnitude, on the scales of the beam size (roughly 4 per cent of the bow shock region in Fig. 2). In future studies, a spatially resolved analysis would be interesting to understand this spatial flux density variations and its relation to the emission’s origin. Similarly, future studies of the spectral shape may probe whether the bow shock region defined for this work, fully traces the region where the shock’s radio emission is dominant.

For geometrical purposes, we will follow the common assumption that the bow shock has a depth that is of the order of its width, Δ. In the radio image, we measure Δ ≈ 35 arcsec, comparable to the standoff distance R₀ ≈ 65 arcsec. Given the surface area of the bow shock region (10 100 arcsec²) and distance to Vela X-1, this translates to an approximate bow shock volume of V_bowshock = 0.28 pc³. The bow shock volume factor, i.e. the fraction of the wind power passing through the bow shock, can be estimated as |$\eta _{\rm vol} \equiv V_{\rm bowshock}/(4\pi R_0^3/3) \approx 0.36$|⁠. If we instead use the approach by De Becker et al. (2017), which compares surface fluxes, we find η_vol ≈ 0.15. This difference arises from the non-circular shape of the bow shock arc. In our following calculation, using a higher η_vol implies lower injection efficiencies (Section 4.3); therefore, our assumption on η_vol makes those calculations more conservative with respect to the conclusions that we draw. In Table 1, we list all stellar (wind) and bow shock parameters that are relevant to the equations shown in the main text. Those that show up in the Appendix only, are listed there. Finally, we stress that the literature contains many estimates of |$\dot{M}_{\rm wind}$|⁠, v_∞, and v_*; we have opted for either representative or recent values and find that none of our qualitative conclusions depend on these exact choices.

4.2 The thermal/free–free scenario

Thermal free–free emission (bremsstrahlung) from a bow shock will depend, fundamentally, on two physical properties of the electron population in the shock: the electron number density n_e and temperature T. These parameters will not only set the radio luminosity (and spectrum) of the system, they also determine the surface brightness of H α line emission. As the Vela X-1 bow shock is detected in both H α and radio images, we can combine the wavebands to infer what electron density and temperature are required to explain the observations.

In the left-hand panel of Fig. 6, we plot electron density as a function of temperature using the observed MeerKAT flux and both literature S_H α constraints. The dashed blue and black lines indicate the factor 2 uncertainty on the S_H α from Kaper et al. (1997), and a |$10{{\ \rm per\ cent}}$| assumed uncertainty for Gvaramadze et al. (2018), who do not report an uncertainty in their work. Measuring the intersections with the red MeerKAT curve, we find |$n_e = 18^{+5}_{-2}$| cm⁻³ and |$T = 1.6_{-1.0}^{+4.0}\times 10^3$| K using Kaper et al. (1997), and |$n_e = 33^{+1}_{-2}$| cm⁻³ and |$T = 3.3^{+0.9}_{-0.6}\times 10^4$| K using Gvaramadze et al. (2018).

Figure 6.

Left: the relation between post-shock electron density n_e and temperature T consistent with the observed radio (red line, this work) and H α (blue dash–dotted line, Isaac Newton Telescope, Kaper et al. 1997 and black dashed line, SuperCOSMOS H α Survey, Gvaramadze et al. 2018). The intersections between the radio and H α curves indicate solutions that can account for both bands. The narrow blue and black lines indicate the (assumed) uncertainty on the H α surface brightness. Right: the same relation between n_e and T for both the bow shock (red line) and the overdense region (red dashed line) to the bottom-left of Vela X-1 in Fig. 3. The blue and grey bands indicate the bow shock temperatures where the radio and H α curves intersect in the left-hand panel.

Open in new tabDownload slide

To assess whether these inferred values are realistic, we can perform several checks. First, we can briefly consider the inferred temperature. Assuming mass, momentum, and energy conversion, we can obtain a post-shock temperature estimate from |$kT \sim (3/16)\mu m_p v_*^2$|⁠, where μ ≈ 0.6 for cosmic abundances and m_p is the proton mass (Helder et al. 2009). Using the stellar velocity of Vela X-1, we infer T ≈ 4 × 10⁴ K, indeed close to the value inferred using the Gvaramadze et al. (2018) H α maps. Secondly, the inferred temperature values are also similar to the range of 6 × 10³–1.4 × 10⁵ K found for a sample of H α-detected bow shocks in Brown & Bomans (2005).

Thirdly, we can use the detection of other large scale, diffuse emission structures with MeerKAT. For instance, the ridge of diffuse emission to the South East of Vela X-1, clearly visible in Fig. 1, is also detected in the SuperCOSMOS H α image. In fact, Gvaramadze et al. (2018) suggest through simulations that this diffuse structure corresponds to an overdense region in the ISM, finding that a factor 3 increase in density compared to the surrounding ISM can explain the observed H α morphology of both the field and the bow shock. In both the MeerKAT radio and SuperCOSMOS H α images, the bow shock shows a higher flux than the ridge, although both are of similar order of magnitude. Therefore, one can imagine a scenario where the radio and H α emission of both structures originates from the same (free–free) emission mechanism, with an enhanced bow shock flux caused by a shock-enhanced density increase.

To test the above scenario, we define a region across the radio-brightest segment of the overdense ridge in the radio image, as shown by the blue dashed contour in Fig. 2. In similar fashion to the bow shock region, we manually defined this arced box region such that it is dominated by ridge emission without containing significant contributions for radio point sources. A slightly different region morphology would, naturally, yield slightly different calculations below; however, this used region samples the mean flux density in the main part of the overdense ridge. We repeat the bow shock calculation for this overdense region, measuring a radio flux density of |$S_{\nu , \rm overdensity} \approx 28.8$| mJy, integrated over the region’s area of ∼ 42 400 arcsec². As we don’t know the 3D geometry of this overdense region, we will follow our rationale for the bow shock and assume its depth is similar to its width of ∼1.8 arcmin ≈2.1 parsec (assuming a similar distance as Vela X-1, supported by the evidence for the interaction between Vela X-1 and the overdense region in Gvaramadze et al. 2018).

Using these geometrical assumptions and the flux measurement in equation (3), we plot the radio constraints on n_e and T for both the bow shock and the overdense region in the right-hand panel of Fig. 6. The blue and grey shaded areas indicate the temperature ranges derived earlier for the bow shock. Under the assumption that the temperature is not significantly changed in the shock, we find that a shock density increase by a factor ∼2.3 would explain the data. Since in reality, the temperature in the shock will change as well, the exact density increase will be slightly different – although we stress that the exact factor scales linearly with the assumed depth of the overdense region along the line of sight, which has not been measured.

In adiabatic shocks, the Rankine–Hugoniot equations imply a density enhancement by a factor of 4 (Landau & Lifshitz 1959). In the analysis by Gvaramadze et al. (2018), a higher density increase, by factors up to ∼9, is inferred from simulations that provide good matches to the observed SuperCOSMOS images. With such enhancement factors, the inferred shock electron densities would imply ISM densities in the range ∼2.0–4.5 cm⁻³ (using Kaper et al. 1997) or ∼3.7–8.3 cm⁻³ (using Gvaramadze et al. 2018). Literature estimates of the ISM density have been made either using measurements of the standoff distance and stellar (wind) properties, or the optical and IR emission of large scale, close-by H ii regions. The former estimates are typically of the order of ∼1–2 cm⁻³ (e.g. Kaper et al. 1997; Peri et al. 2015; Gvaramadze et al. 2018), with higher estimates up to ∼10 cm⁻³ using different stellar wind parameters (e.g. Gvaramadze et al. 2011). The latter type tends to suggest higher ISM densities as well, of the order of ∼5–15 cm⁻³ (Kaper et al. 1997; Lequeux 2005). Therefore, given this range in estimates, the ISM density inferred by combining radio and H α maps appears reasonable, as does the inferred shock density increase (Fig. 6, right). We also note that our inferred densities scale with bow shock depth as Δ^−0.5; therefore, a larger depth would imply lower inferred shock and ISM densities.

Finally, we note that the thermal emissivity scales as |$n_e^2$|⁠; therefore, for shock density enhancements of at least 4 to 9, we expect thermal radio flux densities from the ambient ISM (i.e. outside the overdense region) that are one to two orders of magnitude fainter (assuming, simplistically, similar temperatures and depths). Such flux densities are below our radio RMS sensitivity, which can explain why the diffuse radio emission so closely traces the overdense H α region but no other diffuse radio emission is detected.

The optically thin free–free spectrum has a spectral index of α ∼ 0.1, up to a break frequency where hν ∼ kT. For the inferred temperatures, we therefore expect this cut-off to lie in the optical bands or at lower frequencies, consistent with the non-detection of X-ray emission. However, directly detecting either this cut-off, or bow shock continuum emission, in optical bands appears unfeasible; at the low-frequency end of the optical band (ν ∼ 4 × 10¹⁴ Hz), the expected average flux is only ∼0.08 |$\mu$|Jy arcsec⁻², or ∼26.6 AB mag arcsec⁻². At lower frequencies, the spectrum is either dominated by dust (IR) or too faint and extended to be detectable with current sub-mm observatories (i.e. ALMA).

Finally, we can briefly turn to the only other radio-detected bow shock, BD+43^o3654, and consider it in the same thermal framework. While its integrated flux density is reported in Benaglia et al. (2010) (4.8 GHz) and Benaglia et al. (2021), the total surface area is not. Instead, we can use the flux of brightest beam in the bow shock and consider only that segment of the structure. The radio contours at 4.8 GHz suggest a peak flux density of the order of ∼4 mJy beam⁻¹ (Benaglia et al. 2010), while at 3 GHz, the peak flux is ∼5.5 mJy beam⁻¹ (Benaglia et al. 2021). At 4.8 GHz, the 12 arcsec circular beam, 1.7 kpc distance, and known bow shock width imply a volume of ∼4.4 × 10⁻² pc³. At 3 GHz, the beam size is 20.2 arcsec × 12.5 arcsec, implying instead a volume of ∼7.7 × 10⁻² pc³. At both frequencies, we can then combine the flux, volume, and distance estimates to constrain the electron density and temperature. As shown in Fig. 7, both estimates are quite similar, given the simplifying assumptions in this approach. Furthermore, the implied electron densities are between 1.7 and 1.9 times larger than for Vela X-1. The shock density enhancement, ranges between 3 and 10 for temperatures below ∼10⁵ K, compared to the ISM density 9 cm⁻³ around BD+43^o3654 measured by Benaglia et al. (2021); at a temperature of ∼2.6 × 10⁴ K, corresponding to a stellar velocity of 43.6 km s⁻¹ (Benaglia et al. 2021) and assuming |$kT \sim (3/16)\mu m_p v_*^2$|⁠, this factor lies between 6 and 7.

Figure 7.

Same as Fig. 6, now showing the n_e(T) relation based on the radio detection of Vela X-1 (red) and BD+43^o3654 (black). For the latter, the two lines indicate estimates based on observations at two different radio frequencies (see text for details). The grey band indicates the inferred ISM electron density ranges discussed in Benaglia et al. (2021).

Open in new tabDownload slide

4.3 The non-thermal/synchrotron scenario

Alternatively, the radio emission from the Vela X-1 bow shock could originate from synchrotron emission by a population of relativistic electrons. Such a scenario was first invoked for the radio bow shock in BD+43^o3654 (Benaglia et al. 2010) and has since been developed analytically and numerically to predict the radio and X-ray/γ-ray detectability of other bow shocks (e.g. del Valle & Romero 2012; del Valle & Pohl 2018; del Palacio et al. 2018). In a nutshell, the stellar wind provides the energy budget to accelerate electrons via diffusive shock acceleration, creating a relativistic electron population. This population loses energy via synchrotron emission, inverse-Compton scattering of dust and stellar photons, and relativistic bremsstrahlung. Alternatively, particles can escape the bow shock region via diffusion. The relative importance of the processes depends on geometry, magnetic field, and the stellar and IR radiation fields.

For a subset of the calculations in this Section, we refer the reader to the Appendix for all details. We can start to assess the non-thermal/synchrotron scenario by estimating the shock magnetic field via equipartition arguments, minimizing the combined energy in particles and magnetic fields. In this calculation, and in the remainder of this section, we do not know the radio spectral index α. Therefore, we assume that α ≈ 0.5, as found for BD+43^o3654 (Benaglia et al. 2010). In the Appendix, and several of the figures in this calculation, we will also show the results assuming α ≈ 0.7. However, we note up front that wherever our qualitative conclusions depend on the exact value of α, we will discuss this explicitly.

Further following the analysis performed for BD+43^o3654 and other IR bow shocks (del Valle & Romero 2012; del Valle & Pohl 2018; del Valle, Romero & De Becker 2013; del Palacio et al. 2018), we now turn to estimating the relevant time-scales for electron acceleration, emission losses, and escape. We consider diffusive shock acceleration as the acceleration process. As radiative processes, we include synchrotron radiation, inverse-Compton scattering of electrons in the stellar and IR radiation fields (using the analytical approximations for the Thompson and Klein–Nishina regimes from Khangulyan, Aharonian & Kelner 2014), and relativistic bremsstrahlung. The full calculations are detailed in the Appendix.

The resulting time-scales are plotted, as a function of electron energy and assuming the equipartition magnetic field, in Fig. 8. The maximum electron energy is given by the intersection between the acceleration and shortest loss time-scale, the latter being dominated by convective escape. We find a maximum energy around E_max ∼ 3 × 10¹² eV, similar to our previous assumption. Radiative losses would only become dominant at energies roughly one order of magnitude higher. One notable feature is the similarity of the escape time-scale and the equipartition time-scale, the former only a factor ∼2 larger than the latter for this magnetic field. In other words, roughly half of all available kinetic wind energy should be injected into particle acceleration to ensure a steady-state scenario – a required efficiency that further increases when radiative losses become relevant close to E_max. We will return to this issue of exceptionally high electron acceleration efficiency at the end of this section, fully considering its magnetic field dependence.

Figure 8.

The time-scales of electron acceleration, radiative cooling (synchrotron, inverse Compton, relativistic bremsstrahlung), and convective escape in the bow shock, assuming the equipartition magnetic field B_eq = 30 |$\mu$|G. The full calculations are detailed in the Appendix.

Open in new tabDownload slide

Next, we can consider the broad-band synchrotron spectrum to match the observed radio flux density and Chandra X-ray flux upper limit. The latter non-detection might observationally constrain the maximum electron energy and magnetic field, as the synchrotron cooling break depends on both those quantities. As derived in the Appendix, E_max is set by electron escape for magnetic fields up to a transitional value of ∼84 |$\mu$|G, above which synchrotron losses take over and determine E_max instead. Importantly, this implies that the synchrotron cooling break frequency quadratically increases with magnetic field strength for B < 84 |$\mu$|G, while it remains constant at ν_cool ≈ 6 × 10¹⁶ Hz for stronger magnetic fields. This maximum cooling break frequency lies a factor 2 below the minimum energy in the Chandra band (0.5 keV). Given the steepening of the synchrotron spectrum above the cooling break, this small factor implies that the magnetic field may indeed be constrained by the Chandra non-detection.

In Fig. 9, we show the expected synchrotron spectra assuming either α = 0.5 or α = 0.7, and the magnetic field from either equipartition or the transition between the escape and synchrotron dominated regime (both consistently recalculated in the case of α = 0.7). From these modelled SEDs, we can conclude that the synchrotron spectrum violates the Chandra non-detection for magnetic field strengths in the synchrotron-dominated regime, combined with relatively shallow spectral indices (here, α ∼ 0.5). In other words, the magnetic field should either be closer to equipartition or the spectral index should be steeper. Conversely, for α ∼ 0.7, any arbitrarily high magnetic field will remain consistent with the Chandra non-detection.

Figure 9.

The modelled synchrotron SED for four combination of magnetic field and radio spectral index α, normalized to match the observed MeerKAT radio flux. Only for magnetic fields exceeding the maximum field strength implied by the requirement of compressibility of the stellar wind (equation 10) and relatively shallow spectral indices (i.e. the black dashed line), the SED becomes inconsistent with the Chandra non-detection.

Open in new tabDownload slide

Figure 10.

The injection efficiency of electrons, defined in equation (9), as a function of magnetic field (left) and radio spectral index (right). In the left-hand panel, the efficiency is plotted for different assuming spectral indices, while in the right-hand panel, two different magnetic fields are assumed. The minimum efficiency is obtained for α = 0.5 and B = 46 |$\mu$|G (i.e. the maximum field strength that leaves the wind compressible), yielding η_e ≈ 13 per cent.

Open in new tabDownload slide

These values for the injection efficiency are significantly higher than inferred in comparable systems (del Palacio et al. 2021): for instance, del Palacio et al. (2018) infer values between 16 per cent and |$0.4{{\ \rm per\ cent}}$| assuming ζ_B = 0.01 and ζ_B = 1, respectively, for BD+43^o3654. Similarly applying two extreme values of ζ_B, Benaglia et al. (2021) find efficiencies of 10 per cent (ζ_B = 0.03) and 1.5 per cent (ζ_B = 1.0) for the same source. Indeed, we can apply equation (9) to the brightest beam of BD+43^o3654 in Benaglia et al. (2010) (see the end of Section 4.2) assuming the stellar wind parameters reported in that work: this confirms that the required injection efficiencies in BD+43^o3654 are one order of magnitude lower for the same spectral index α. Further comparisons can be made with Stappers et al. (2003), who report efficiencies of ∼4 per cent and <9 per cent for the Crab nebula and the pulsar wind nebula around B1957+20, respectively.

These high injection efficiencies required in the Vela X-1 bow shock are not the result of assuming incorrect values of α either. While α might be different than 0.5 or 0.7, the injection efficiency does not depend monotonically on spectral index. In the right-hand panel of Fig. 10, we plot η_e versus α for two different magnetic field (equipartition for α = 0.5 and assuming ζ_B = 1). In both cases, the efficiency is minimized around α = 0.5. One can also see how a higher magnetic field lowers the efficiency, implying that η_e(B = B_max, α = 0.5) = 13 per cent is the lowest possible value for Vela X-1. Finally, for significantly lower α, the non-detection of X-ray emission requires a lower magnetic field (e.g. Figs A2 and 9), further increasing the inferred injection efficiency.

We can conclude that only a fine-tuned combination of the electron density distribution slope and magnetic field yields an injection efficiency approaching the maximum values inferred for BD+43^o3654 and other types of nebula. Instead, without such fine-tuning, significantly higher values up to 100 per cent are required to explain the Vela X-1 bow shock radio emission fully through a non-thermal scenario.

Alternatively, we can turn this argument around: if the injection efficiency is in reality a factor 10 lower than those derived above (i.e. ∼1.3 per cent at minimum), one-tenth of the observed radio flux could be attributed to synchrotron processes. As the electron densities in the thermal scenario scale with |$S_\nu ^{0.5}$|⁠, this would imply only a ∼5 per cent decrease in the inferred n_e. Such a decrease does not invalidate the arguments in Section 4.2, and is smaller than the systematic differences between the densities inferred using the Kaper et al. (1997) and Gvaramadze et al. (2018) H α maps. Therefore, a dual thermal + non-thermal scenario, with the observed flux dominated by the former, is a realistic possibility.

4.4 Why are only the bow shocks of BD+43^o3654 and Vela X-1 detected in radio?

The MeerKAT detection of the Vela X-1 bow shock is only the second radio bow shock of a runaway massive star after the discovery of the prototype in BD+43^o3654 (Benaglia et al. 2010), and the first radio bow shock around an HMXB. Radio searches for other bow shocks have since been performed using the NRAO VLA Sky Survey (NVSS; Condon et al. 1998) for HIP 16518, HIP 34536, HIP 78401, HIP 97796 (De Becker et al. 2017) and pointed VLA observations (as mentioned by Benaglia et al. 2021, based on private communication with C. Peri). However, these searches have all remained unsuccessful. Considering both the thermal or non-thermal scenario for the radio emission (or their combination), we will now briefly discuss the future prospects for radio bow shock detections and whether exceptional physical conditions are required.

First, we can consider the non-thermal scenario. The four stars that De Becker et al. (2017) did not detect in radio, can be compared to Vela X-1 and BD+43^o3654 in the frame work introduced in Section 4.3. For these sources, De Becker et al. (2017) collect distances, bow shock widths, stellar mass-loss rates, and wind velocities. For these, we can assume a maximum electron energy of 10¹² eV and a high injection efficiency of η_e = 10 per cent. Using these numbers, we can estimate the expected flux density in a 1.4-GHz, 45 arcsec beam of the NVSS. For these estimates, we assume a best-case scenario of α = 0.5, corresponding to the highest S_ν for a given η_e (Fig. 10, right), and consider two high values of B (50 and 100 |$\mu$|G). In all cases, the bow shock is larger than the beam size.

For all four stars, non-thermal flux densities in the range of 0.4–1.5 mJy beam⁻¹ and 1.0–4.3 mJy beam⁻¹ can be expected, assuming B = 50 and 100 |$\mu$|G, respectively. Given the typical RMS of the survey of ∼0.5 mJy beam⁻¹, the non-detection of radio emission from these bow shocks is unsurprising: only in the best-case circumstances in terms of spectral shape, magnetic field, and injection efficiency, only HIP 34536 might be 3σ detectable for the strongest considered magnetic field. However, even then, at such low significance, likely no bow shock morphology will be identifiable. These sources are not outliers within the E-BOSS bow shock samples, nor are they similar in their relevant properties: they span wind velocities between 500 and 2500 km s⁻¹, mass-loss rates between 6 × 10⁻⁹ to 5 × 10⁻⁷|$\rm M_{\odot }$| yr⁻¹, and distances between 0.22 and 2.2 kpc. Therefore, we deem it reasonable to extend this conclusion to the wider bow shock sample.

Similarly, we can investigate a thermal scenario. For all four systems considered above, the first E-BOSS catalogue lists inferred ISM densities. Using the same geometrical assumptions as in the non-thermal scenario, and assuming a temperature, we can again predict the expected radio flux in the NVSS. For these four systems, the ISM densities range between 0.01 and 2 cm⁻³. We can assume a shock density enhancement of 4 and, more extremely, 10 (Gvaramadze et al. 2018). However, we estimate flux densities of only 0.2 |$\mu$|Jy beam⁻¹ to 0.8 mJy beam⁻¹ in the latter case, assuming a temperature of T = 10⁴ K. Even at a lower temperature of 10³ K, the maximum flux density is 1.85 mJy beam⁻¹, for HIP 78401, which has the highest inferred n_e and smallest distance. Therefore, we also do not expect any significantly detectable radio emission in the NVSS in the thermal scenario.

In the thermal scenario, the most effective band to search for emission depends on the electron temperature. Comparing equations (1) and (4), it is clear that the ratio of radio flux density and H α surface brightness is independent of the electron density. However, their ratio will depend strongly on temperature, i.e. S_ν/S_H α ∝ T^0.4g(ν, T) (ignoring the exponent in equation 1). In the temperature range between 10² and 10⁵ K, this ratio increases by a factor 50, implying that the radio continuum emission may be easier to pick up at higher temperatures. However, the absolute flux density/surface brightness depends very strongly on electron density as well; a relatively hot but low density medium will remain undetectable in both bands.

What then makes BD+43^o3654 and Vela X-1 stand out? For Vela X-1, the high sensitivity and core-heavy configurations of MeerKAT are vital to detect the bow shock’s radio emission. However, as discussed in Section 4.2, the presence of a local overdense ISM structure (Gvaramadze et al. 2018) is as essential in a free–free scenario. In BD+43^o3654, on the other hand, both a relatively high local ISM density and a high kinetic wind power could be the crucial factor: in our injection efficiency estimates, we assumed a stellar mass-loss rate of 1.6 × 10⁻⁴|$\rm M_{\odot }$| yr⁻¹ and a velocity of 2300 km s⁻¹(Benaglia et al. 2010; Kobulnicky, Gilbert & Kiminki 2010), leading to significantly higher wind powers than in the four stars considered above. It may be that, indeed, either exceptional stellar wind or ISM properties are essential for a radio bow shock detection.

The prospect of recent and upcoming all-sky radio surveys has several consequences for radio bow shock studies. One example of such a survey is the aforementioned RACS (McConnell et al. 2020) that covers the entire sky South of +40^o declination at 887.5 MHz down to a typical RMS sensitivity of 0.25 mJy beam⁻¹. For RACS, the ∼10 times smaller beam area combined with the ∼2 times better sensitivity implies that the flux detectability of individual bow shocks (as considered above) does not greatly improve compared to NVSS. However, given the uniform sky coverage, especially in the Galactic plane and centre regions, a large number of known bow shocks can be investigated further in the radio band. The ongoing Galactic plane survey with MeerKAT (Goedhart et al., in preparation) will provide similar sky coverage but at higher sensitivity and, as our detection of the Vela X-1 bow shock shows, with an array configuration especially suitable for bow shock studies. Crucially, with both surveys, any bow shock that is radio detected will be better resolved than it would have been with the NVSS, helping to distinguish between bow shocks and other morphologies. Especially in a thermal scenario with a high surrounding ISM density, comparing the radio morphology to IR maps may be vital to identify radio counterparts. However, based on our considerations, detecting radio bow shocks in these large-scale surveys may still require either high stellar wind kinetic powers (such as BD+43^o3654) or locally overdense or complex ISM environments (such as Vela X-1).

5 SUMMARY AND CONCLUSIONS

In this paper, we have presented the discovery of 1.3-GHz radio emission from the Vela X-1 bow shock with MeerKAT. The MeerKAT data also reveal other large-scale structures of diffuse radio emission, tracing the known H α structures in the field. Our analysis of publicly available X-ray and lower frequency radio observations does not reveal bow shock emission at these other frequencies. These results present only the second radio bow shock around a massive runaway star, compared to over 700 known bow shocks observed in the IR band. It also presents the first radio bow shock detected around an X-ray binary.

Lacking a radio spectral index measurement, we turn to the bow shock’s energetics and brightness to assess the underlying emission mechanism. We first consider a thermal scenario, dominated by optically thin free–free emission. In such a scenario, the radio and H α emission of both the bow shock and diffuse structures, originates from the same process and particle population. Combining their constraints, we find reasonable estimates for the bow shock and ISM density and temperature, consistent within their systematics with earlier, independent estimates. Alternatively, a non-thermal scenario, dominated by synchrotron emission form a shock-accelerated population of electrons, is harder to reconcile with the observed bow shock properties. In particular, it requires high energy injection efficiencies of ≳ 13 per cent, depending on the exact electron distribution and magnetic field. The observed emission may, of course, in reality originate from a combination of both processes, dominated by the thermal emission.

Finally, we consider why the great majority of stellar bow shocks have escaped detection at radio frequencies so far. Our considerations regarding the thermal and non-thermal scenarios, suggest that either high density/complex ISM environments or exceptionally energetic stellar winds, respectively, are required for such a detection. However, as our MeerKAT observations show, the advent of SKA precursors and the future SKA may allow for the detection of a significant number of other stellar bow shocks. Building up such a sample will test the hypothesis above, and may reveal new examples of synchrotron-dominated bow shocks as evidence of electron acceleration to very high energies.

ACKNOWLEDGEMENTS

We thank the anonymous referee for their helpful report that improved the quality and clarity of this work. We thank the staff at the South African Radio Astronomy Observatory (SARAO) for scheduling these observations. The MeerKAT telescope is operated by the South African Radio Astronomy Observatory, which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation. This work was carried out in part using facilities and data processing pipelines developed at the Inter-University Institute for Data Intensive Astronomy (IDIA). IDIA is a partnership of the Universities of Cape Town, of the Western Cape and of Pretoria. We acknowledge the use of data obtained from the High Energy Astrophysics Science Archive Research Center (HEASARC), provided by NASA’s Goddard Space Flight Center. This paper includes archived data obtained through the CSIRO ASKAP Science Data Archive, CASDA (https://data.csiro.au). The Australian SKA Pathfinder is part of the Australia Telescope National Facility (grid.421683.a) which is managed by CSIRO. Operation of ASKAP is funded by the Australian Government with support from the National Collaborative Research Infrastructure Strategy. ASKAP uses the resources of the Pawsey Supercomputing Centre. Establishment of ASKAP, the Murchison Radio-astronomy Observatory and the Pawsey Supercomputing Centre are initiatives of the Australian Government, with support from the Government of Western Australia and the Science and Industry Endowment Fund. We acknowledge the Wajarri Yamatji people as the traditional owners of the Observatory site. This work makes use of several python packages, namely numpy (Oliphant 2006), astropy (Astropy Collaboration 2013, 2018), matplotlib (Hunter 2007), and aplpy (Robitaille & Bressert 2012). JvdE is supported by a Lee Hysan Junior Research Fellowship awarded by St. Hilda’s College, Oxford. TDR acknowledges financial contribution from ASI-INAF n.2017-14-H.0, an INAF main stream grant. GRS is supported by NSERC Discovery Grants RGPIN-2016-06569 and RGPIN-2021-04001. SEM acknowledges financial support from the Violette and Samuel Glasstone Research Fellowship programme, the UK Science and Technology Facilities Council (STFC), and Oxford Centre for Astrophysical Surveys, which is funded through generous support from the Hintze Family Charitable Foundation. SEM also acknowledges financial contribution from the agreement ASI-INAF n.2017-14-H.0 and the INAF mainstream grant, and from PRIN-INAF 2019 n.15. SM acknowledges funding from the South African National Research Foundation (NRF) and University of Cape Town VC2030 Future Leaders Award. JCAM-J is the recipient of Australian Research Council Future Fellowship (project number FT140101082).

DATA AVAILABILITY

A Jupyter notebook to generate the figures in this work and redo the analysis in the Discussion, can be found on this link upon publication: https://github.com/jvandeneijnden/VelaX-1-MKAT. There, the underlying MeerKAT and Chandra images are also included. The raw MeerKAT and Chandra data, as well as the RACS images, are publicly available in their respective online observatory repositories.

Footnotes

REFERENCES

APPENDIX A: NON-THERMAL CALCULATIONS

In this appendix, we will provide addition derivations of the calculations presented in Section 4.3 of the main paper. As noted there, the below calculations are build upon the pioneering work of Benaglia et al. (2010), del Valle & Romero (2012), del Valle & Pohl (2018), and del Palacio et al. (2018). We refer the reader to those works for more details underlying the calculations.

A1 Equipartition and particle energies

Figure A1.

The energy contained in relativistic particles and the magnetic field, as a function of the magnetic field strength, assuming an equal number of protons and electrons. For a spectral index α = 0.5, we measure an equipartition magnetic field of B_eq ≈ 30 |$\mu$|G.

Open in new tabDownload slide

A2 Loss and acceleration time-scales

Finally, we calculate the relativistic bremsstrahlung via equation (11) in del Valle & Romero (2012), assuming that the wind density is enhanced by a factor 4 in the shock. The time-scale of escape, the final loss mechanism, is energy independent and is estimated as t_escape ≈ Δ/v_∞. Following common assumptions, we model the acceleration time-scale assuming diffusive shock acceleration, given in Terada et al. (2012) as (where we correct the factor c to yield the correct units) t_acc = (20ξ/3) × (E/eB) × (1/v_∞)². We also assume Bohm diffusion, i.e. ξ = 1.

A3 The broad-band synchrotron spectrum

To model the broad-band synchrotron spectrum, we (i) assume a spectral index, (ii) match the observed MeerKAT L-band radio flux density to the SED, and (iii) calculate the low- and high-frequency behaviour due to synchrotron-self-absorption and the maximum electron energy, respectively. We can then also compare this to the X-ray upper limit: this non-detection may specifically constrain the maximum electron energy and magnetic field, as the high-frequency break depends on both those quantities. However, E_max and B are themselves also related: at equipartition, Fig. 8 shows that electron escape sets the maximum energy. For different magnetic field strengths, the acceleration time-scale will decrease as B⁻¹, while the synchrotron cooling time-scale decreases as B⁻². As escape and inverse-Compton processes are all independent of the magnetic field, at some B, synchrotron losses will set the maximum energy and, in turn, the cooling break.

Figure A2.

The maximum electron energy (top) and resulting synchrotron cooling break frequency (bottom) as a function of magnetic field strength. Escape and synchrotron indicate the dominant energy loss mechanism, setting the maximum electron energy. At no point does the cooling break frequency extend into or above the Chandra band.

Open in new tabDownload slide

A4 The injection efficiency