© The Authors 2025

1 Introduction

Bow shocks arising from the hydrodynamic interaction of the interstellar medium (ISM) with stellar wind from runaway massive stars were first detected by Gull & Sofia (1979) using optical nebular emission lines. The detection, however, is difficult because O and B stars are much brighter than their nebulae at optical wavelengths, and only a few optical detections can be found in the literature (e.g. Gvaramadze & Bomans 2008; Gvaramadze et al. 2014).

The IR wavelengths are arguably the best part of the electromagnetic spectrum to observe the bow shocks of massive runaway stars. Meyer et al. (2014) showed that bow shocks are most easily detectable at mid-IR because they radiate the most energy at these wavelengths. This detectability can be seen clearly in the IR surveys from van Buren & McCray (1988), van Buren et al. (1990), Peri et al. (2012, 2015), and Kobulnicky et al. (2016, 2017). The early surveys clearly defined the different morphological and physical classes of objects and their characteristics using results from the Infrared Astronomical Satellite (IRAS) all-sky survey. Peri et al. (2012, 2015) presented results from their Extensive stellar BOw Shock Survey (E-BOSS), which followed up on suitable candidates from the IRAS survey and confirmed many of their candidates to be bow shocks from massive runaway stars. Kobulnicky et al. (2016, 2017) showcased over 600 bow-shock candidates at IR wavelengths, providing an excellent library of references for future observations and studies. Using Gaia DR3 to identify possible runaway stars, Carretero-Castrillo et al. (2025) were able to detect new mid-IR bow-shock candidates with WISE data, some of which have hints of radio emission coincident with the bow shock.

Studies at IR wavelengths are useful, but they have limitations. In particular, since the bow shocks are excellent candidates for high-energy particle acceleration (e.g. del Palacio et al. 2018; del Valle & Romero 2012), it is important to observe at photon wavelengths that are sensitive to both thermal emission and non-thermal emission. This is not possible with IR because of the very strong thermal emission from dust. At high energies (X-rays and γ-rays), even though it was thought that the non-thermal emission can be detected in the arc of the bow shocks, only upper limits have been obtained (e.g. Schulz et al. 2014; Toalá et al. 2016, 2017; H.E.S.S. Collaboration 2018).

Radio wavelengths have proven to be the best part of the electromagnetic spectrum for detecting thermal and non-thermal emission from bow shocks of massive runaway stars. What is more, in radio we can measure gas directly and not dust. This can be seen by the increasing numbers of detections over the past decade, including Benaglia et al. (2010), Benaglia et al. (2021), Moutzouri et al. (2022, hereafter Paper I), Van den Eijnden et al. (2022b), and van den Eijnden et al. (2024).

Non-thermal emission has been detected in at least three bow shocks around the stars BD+43^° 3654 (Benaglia et al. 2010), BD+60^° 2522 (Paper I), and LS 2355 (van den Eijnden et al. 2024). Overall, convincing radio detections have been reported for seven bow shocks around the stars BD+43^° 3654 (Benaglia et al. 2010), RCW 49 S1 (Peri et al. 2015; Van den Eijnden et al. 2022b), NGC 6357 G1 (Van den Eijnden et al. 2022b), NGC 6357 G3 (Peri et al. 2015; Van den Eijnden et al. 2022b), BD+60^° 2522 (Paper I), Vela X-1 (van den Eijnden et al. 2022a), and LS 2355 (van den Eijnden et al. 2024). A tentative detection was reported for EB23: HIP 88652 (Peri et al. 2015; Van den Eijnden et al. 2022b), and Van den Eijnden et al. (2022b) reported likely detections of three further bow shocks embedded in larger-scale diffuse emission.

Radio detections are primarily enabled by upgraded radio infrastructure that provide a higher sensitivity and the capacity to push the limits of observing in the radio with both wide-field and wide-band capabilities. Recognising these new capabilities in radio astronomy, we undertook a survey of bow shocks in 2019 with the National Science Foundation’s (NSF) Karl G. Jansky Very Large Array (VLA) while complementing the dataset with single-dish observations from the 100 m Effelsberg Telescope for ten of the sources in 2020 and 2021. Following the presentation of the first results from our survey in Paper I, here we report on the complete survey of 24 bow shock candidates, including the sources BD+43^° 3654 and BD+60^° 2522 that were discussed in that previous work.

Interest in measuring gaseous emission from bow shocks arises partly because of their potential use as laboratories for collisionless shock physics, including particle acceleration, magnetic field amplification, and non-thermal radiation (del Valle & Romero 2012). They are also excellent laboratories for studying turbulent mixing at shear interfaces between the hot and warm ISM phases (Tan et al. 2021; Mackey et al. 2025). They may also be used to measure the mass-loss rates of O and B stars in cases where uncertainties from stellar spectroscopy may be very large (Gvaramadze et al. 2012). Nearby bow shocks such as that of ζ Oph (Gvaramadze et al. 2012; Green et al. 2022) provide a rare opportunity to investigate hydrodynamic and radiative processes at the wind-ISM interface at high spatial resolution as constrained by observations (Toalá et al. 2016). For these investigations to be successful, the gas density and ram pressure in the wind and ISM (from the reference frame of the driving star) should be well constrained by measurements.

Gaia DR3 (Gaia Collaboration 2023) gives reliable distances and proper motions for a large number of Galactic runaway massive stars and, if combined with a density measurement in the bow shock, these data give the ISM ram pressure. Computer modelling of bow shocks is advancing to the point that realistic 3D magnetohydrodynamic simulations can be run to model specific systems such as ζ Oph (Green et al. 2022) and λ Cep (Scherer et al. 2020) or generic studies of the properties of magnetised bow shocks (Mackey et al. 2025; del Valle et al. 2025). Similarly, synthetic observations of thermal (Green et al. 2022; del Valle et al. 2025) and non-thermal (del Palacio et al. 2018; del Valle & Pohl 2018; del Valle et al. 2025) radio emission based on multi-zone or magnetohydrodynamic models are gaining predictive power. It is therefore timely and valuable to build up a statistical sample of radio-detected bow shocks where the gas density may be measured from observations and the relative contributions of thermal and non-thermal radiation can be quantified by multi-frequency datasets.

The present work is organised as follows: Section 2 presents the selection criteria of the targets, the details of the observations with the VLA and with Effelsberg, and the methods for calculating proper motions of the driving sources of the candidate bow shocks where this was possible. Results are presented in Section 3 and discussed in the context of previous radio detections in Section 4. Section 5 outlines our conclusions and perspectives for future follow-up observations.

2 Methods

2.1 Criteria for target selection

The main aim of our survey was to observe bow shock candidates with the VLA. Therefore, we had to choose targets suitable for observing from that specific location. We used the following catalogues of bow shock candidates: Kobulnicky et al. (2016, 2017) and Peri et al. (2012, 2015). Observations were made during the 2019B Semester when the VLA was in D-configuration, choosing the C band (4–8 GHz) and X band (8–12 GHz).

The biggest limitation was to choose candidates that have declination >−16^° in order to minimise shadowing from the antennas¹ for sources that are low in the sky. Since the D-configuration is the tightest layout for the VLA, there is a high chance the targets with declination ≤−16^° will be overshadowed. This limitation excluded almost half of the candidates.

Next, we had to consider the limitations on angular resolution. The smallest resolvable angle, θ_min, depends highly on the observing wavelength, λ, and the arrangement of the array: θ_min ∼ λ/B_max, where B_max is the length of the largest baseline of the interferometer. The largest observable scale, θ_max, is set by the beamsize of a single component telescope of the interferometer with diameter D_tel, according to θ_max ∼ λ/D_tel, For D-configuration, the observing limits (i.e, observable angular scales of the target²) are between 12^′′ and 240^′′ for C band, and between 7.2^′′ and 145^′′ for the X band. If the target is outside those limits, the array will not be able to resolve the image or will lose flux on large scales. To decide which of the remaining objects we would include in the observation, we measured the size of the bow shocks from the available IR maps online³.

From those that had the appropriate size, we chose a sample of the brightest ones, and added three extra targets of special interest, namely the bow shocks of BD+43^° 3654, BD+60^° 2522 and ζ Oph. Our final list of target bow shocks include 12 originating from O stars, 4 from B stars, and 8 from stars of currently unknown spectral type. Of these targets, BD+43^° 3654 (Benaglia et al. 2010, 2021) and EB23: HIP 88652 (Peri et al. 2015; Van den Eijnden et al. 2022b) have been previously detected in radio, albeit only a tentative detection for EB23: HIP 88652. Note that for EB31: HD 57682, a typo in Table 8 of Peri et al. (2015) incorrectly labelled the bow shock EB31 as being driven by the star HIP 57862 instead of the correct star, HD 57682 (P. Benaglia, private communication). The list of target bow shocks with coordinates and information on their driving stars is given in Table A.1.

2.2 VLA observations

The extended capabilities of the VLA (Perley et al. 2011), operated by the National Radio Astronomy Observatory, made the telescope an ideal instrument for our observations. We investigated our sources with radio observations at 4 to 12GHz using the VLA (Project ID: 19B-105, PI: M. Moutzouri), applying several techniques to analyse their spectra. The data were reduced and analysed in the same manner as in Paper I, briefly summarised in the following.

The 24 targets in our sample were observed across eight separate observational runs. Each run employed a standard flux density and bandpass calibrator, along with one to four phase and amplitude calibrators. A summary of the observed targets and the associated calibrators for each run is provide in Table 1. The raw data acquired from the observations were reduced using the VLA calibration pipeline included in the Common Astronomy Software Applications (CASA, McMullin et al. 2007). We used the CASA version 5.7 (6.1 for the pipeline) and the National Radio Astronomy Observatory’s (NRAO)⁴ facilities in order to reduce our data.

We imaged each target for the full frequency range of 4 GHz to 12 GHz, comprising all 64 spectral windows of 128 MHz each, using the CASA task ‘tclean’ with the multi-scale multi-frequency deconvolution algorithm (Rau & Cornwell 2011). The images were cleaned for a diameter of four times the FWHM of the primary beam, corresponding to a field of view of 21^′ for central frequency of 8 GHz. We used a cell size of 0.7^′′ and scales sizes of [0, 5, 10, 15, 20, 25, 50, 100, 200] pixels, with a small scale bias⁵ of 0.9. Scales that were too large were ignored by the software.

We also used the multi-scale multi-frequency deconvolution algorithm to calculate a spectral index map by modelling the spectrum of each flux component. These were manually corrected for the primary beam response with the task ‘wide-bandpbcor’ as in Paper I. Next, we applied a mask of 5σ (where σ is the RMS noise level of the intensity map, quoted in Table 2) so that we only plot the spectral index where the diffuse emission is clearly detected above the noise level. These steps were repeated for the spectral index error maps.

2.3 Effelsberg observations

We conducted radio continuum imaging observations in the 4–8 GHz range towards ten selected targets using the 100 m Effelsberg radio telescope in Germany⁶. Nine of the targets were selected on the basis that we detected emission from the direction of the bow shock in a preliminary data reduction of the VLA dataset. This includes all of our claimed and possible detections in this paper except for the sources EB23: HIP 88652 and G075.1730-00.5964, which are marginal cases that did not show up clearly in the preliminary data reduction. To these we added ζ Oph because its large angular extent makes it a good target for Effelsberg, with the possibility to spatially resolve the emission. The observations were performed between 2020 September 16 and 19 (Project ID: 86-20; PI: M. Moutzouri). The 4.5 cm broadband secondary-focus receiver, equipped with two orthogonal linear polarisations, served as the front end, while the SPEctro-POLarimeter (SPECPOL) backend was employed. SPECPOL provides two frequency bands, covering 4–6 GHz (lower band) and 6–8 GHz (upper band), each split into 1024 channels with a channel width of 1.95 MHz.

Observations were carried out in on-the-fly mapping mode, using a scanning speed of 60^′′ s⁻¹ and a step size of 30^′′ to satisfy the Nyquist sampling criterion. Each source was mapped in both right ascension and declination directions to mitigate scanning artefacts using the basket-weaving technique (e.g. Müller et al. 2017). The telescope focus was typically verified after sunrise and sunset, and pointing calibrations were performed every 2–3 hours using nearby bright continuum point sources. The pointing accuracy was better than 10^′′.

The flux density scale was calibrated using observations of 3C 286 and NGC 7027. The radio continuum data were reduced with the toolbox program and the NOD3 software package (Müller et al. 2017). Following the procedure described in Paper I, we generated continuum images at seven central frequencies that are 4.487, 4.743, 4.999, 5.383, 5.639, 6.615, and 7.639 GHz. Radio-frequency interference was visually inspected and manually flagged. Each image covers a bandwidth of about 200–300 MHz with typical 1σ noise levels of ∼6 mJy. The angular resolution is about 140^′′, at 5 GHz. Of the ten sources observed, BD+43^°3654 and NGC 7635 were previously presented in Paper I, while the remaining sources are reported in this work.

2.4 Calculation of peculiar proper motion for stars

We estimated the local proper motions of the targets using surrounding stars as reference benchmarks to identify any peculiar motion relative to the distribution of the nearby objects. These calculations are based on the proper motions and parallaxes provided by Gaia DR3 (Gaia Collaboration 2023). However, some of the faintest stars do not have Gaia counterparts in the current release. Additionally, the parallaxes and proper motions of certain objects exhibit large uncertainties. These targets have values ‘n/a’ in Table A.1 for their proper motion.

The local reference system was established using all stars provided by Gaia within 10^′ of each target. From the projected stars within this 10^′ aperture, we selected those located within a local volume of 250 pc around the target. Parallax was used as the distance estimate for both the target and the local reference volume. The local proper motions of each target, listed in Table A.1 (column Peculiar PM), are with respect to the average proper motions of the selected local volume, also listed in Table A.1 (column Env. PM).

2.5 Estimating the gas density from radio emission and upper limits

Considering thermal bremsstrahlung emission, we can obtain a constraint on the emission measure of the bow shock and thereby on the electron density, n_e, by assuming smoothness and a simple geometry for the bow shock. If we assume all of the detected emission is bremsstrahlung then we obtain an upper limit on n_e because part of the emission could be non-thermal. Non-detections are upper limits by default.

Furthermore, hot stars are surrounded by H II regions of photoionised gas, which also emit bremsstrahlung that contributes to the line-of-sight emission measure, implying that only part of the detected emission is from the bow shock. The degree of contamination depends on the relative size of the bow shock and H II region, as well as the compression factor of the forward shock (because the recombination radiation scaled with n²).

Following van den Eijnden et al. (2022a) and Spitzer (1978), the bremsstrahlung emissivity, j_ν, can be expressed as a function of the gas temperature, T; electron and ion number densities, n_e and n_i; ionic charge, Z_i; and the observation frequency, ν:

$$\frac{j_v}{n_e{n}_i}=5.44\times {10}^{-39}\frac{g_{\text{ff}}{Z}_i}{\sqrt{T}}\exp \left(-\frac{hv}{kT}\right){\mathrm{ergcm}}^3 {\text{s}}^{-1}{\text{sr}}^{-1} {\text{Hz}}^{-1},$$(1)

where g_ff is the Gaunt factor, which depends on T, Z_i, and ν in the radio regime, and where hν/kT ≪ 1. For the case considered below, g_ff ≈ 5.

We also assumed the following:

1. The bow-shock emission is spatially resolved.

2. The density within the layer of the shocked ISM is constant.

3. Within the same shocked layer, there is a constant temperature of T = 8000 K, which is appropriate for isothermal shocks in photoionised gas.

4. The bremsstrahlung emission is dominated by H⁺ ions and electrons so that Z_i = 1 and n_e ≈ n_i.

5. The path length of a ray through the shocked layer near the apex is approximately the standoff distance of the bow shock from the star, R₀.

6. The emission is optically thin.

All of these assumptions are approximately valid for the VLA observations of bow shocks, for stars moving within the diffuse ISM such that the forward shock is radiative but not so dense that radio emission becomes optically thick. The assumptions allowed us to invert the radiative transfer equation for the specific intensity, I_ν = j_νdℓ, (where ℓ is the path of the ray, and absorption is ignored) to obtain a crude estimate of the mean value of n_e within the shocked layer, or an upper limit in case the bow shock is not detected in radio. We obtained

$$n_e=1.80\times {10}^{-21}\sqrt{\frac{n_e{n}_i}{j_v(T,Z,v)}\frac{\text{pc}}{R_0}\frac{sr}{{\theta }_{\text{B}}}\frac{I_v}{{\text{Jy beam}}^{-1}}} {\text{cm}}^{-3},$$(2)

where θ_B is the solid angle of the telescope beam, which is in the range θ_B ∈ [1, 3.5] × 10⁻⁹ sr for our VLA observations.

For the stars with a known spectral type (17 of the 24 targets), all of which also have a measured distance, we could also estimate the pre-shock gas density from the measured standoff distance, R₀, given by

$$R_0=\sqrt{\frac{\dot{M}{v}_{\infty }}{4\pi {\rho }_0\left(v_{\star }^2+c_s^2\right)}},$$(3)

where v_∞ is the terminal velocity of the stellar wind, ρ₀ is the ISM gas density, and c_s ≈ 12 km s⁻¹ is the ISM sound speed in photoionised gas. The ISM electron density, n_e,0 can be obtained from ρ₀ assuming 1.1 electrons per H nucleus (if He is singly ionised) and the ISM mass fraction of H of X_H ≈ 0.7. We have measured R₀ and v_⋆ for these 17 targets using Gaia DR3 distances and proper motions, so an estimation of $\dot{M}{v}_{\infty }$ will give n_e,0.

We estimated the stellar luminosity, mass, temperature, and radius using the tables of Martins et al. (2005) for the O stars, and from these we obtained the mass-loss rates by using the prescription of Vink et al. (2001) (assuming solar metallicity). Wind velocities were obtained from the escape velocity following Eldridge et al. (2006). For the B stars, we used the tables of Krtička (2014) to get mass-loss rates and wind velocities for main-sequence stars and those of Krtička et al. (2021) to get the same quantities for supergiants. G108.0269-00.3497 has no luminosity class, so we assumed a supergiant, and for both it and G107.1371-00.6974, we took $\dot{M}$ and v for the 25 M_⊙ calculations. For ζ Oph, EB23: HIP 88652, BD+43° 3654, and BD+60^° 2522 we took mass-loss rates from literature estimates (Gvaramadze et al. 2012; Gvaramadze & Bomans 2008; Moutzouri et al. 2022; Green et al. 2019). The mass-loss rates and wind velocities we obtained are listed in Table B.1.

3 Results

The 24 targeted bow shocks and bow-shock candidates are listed in Table A.1, sorted by increasing RA of the driving star. They are identified according to their name in the E-BOSS II (Peri et al. 2012, 2015) and Kobulnicky et al. (2016) catalogues in the first column, and the ID of the driving star is listed in the second column. Of the 24 targets, six were convincingly detected as extended radio sources with the VLA, five show some hints of emission, and 13 are undetected.

3.1 Detections

Intensity images from the VLA for the six clear detections are shown in Figure 1. Given the small number of radio-detected bow shocks currently known in the literature, this is an excellent return of positive detections for the survey. They are discussed one by one in the following.

3.1.1 G018.2660-00.2988

This bow shock candidate is entry 72 in the catalogue of Kobulnicky et al. (2016), driven by the O6V-O5V (Ji et al. 2012) star 2MASS J18251808-1309427. This star is part of a young star cluster creating the interstellar bubble N 22 (Ji et al. 2012) and the larger H II region RCW 166. The bow shock is observed with a width ∼50^′′ (shown in the third row of Figure 1), although it is quite possible that its features are obscured by the more extended bubble that surrounds it. Indeed, the radio emission that we see in Figure 1 seems to correlate better with the nebular emission in the DSS Red channel than with the IR arc, implying that we are detecting diffuse emission from the H II region as well as the putative bow shock. The peculiar proper motion of ≈25 km s⁻¹ to the NE is approximately what one would expect based on the orientation of the bow shock, but the uncertainties in the direction are quite large.

The spectral index map shows that the source emits a mixture of thermal emission and non-thermal emission (third row, second panel of Figure 1) although the uncertainties are again rather large (third row, third panel of Figure 1), particularly in regions with nominally very negative spectral index. At the position of the IR arc, the spectral index could be consistent with both thermal and non-thermal emission, with −0.5 ≲ α ≲ 0. Further investigation on the north part of the source could be made, although the IR arc itself is so small and embedded in the bubble that a definitive conclusion on its nature may not be easy to obtain. The Effelsberg image (third row, fourth panel of Figure 1) shows that the bubble is slightly resolved spatially, with the brighter north side clearly also brighter with Effelsberg. Unfortunately the bow shock is much too small to be seen in the Effelsberg image and is superimposed on the H II region emission.

With the VLA we detect a clear lack of emission from the central part of the H II region, where most of the O stars in the young star cluster are located (Ji et al. 2012). The most likely explanation is that the O stars have evacuated a low-density cavity through the combined actions of their stellar winds. In that case the IR arc could represent the edge of the wind bubble (Mackey et al. 2016) and may not be associated with a significant gas overdensity. Alternatively, we may be seeing a young O star that has just been ejected from its birth cluster and is driving a bow shock at the edge of the cluster wind-cavity.

This young star cluster is located in a region of bright gamma-ray emission associated with pulsars and pulsar-wind nebulae (H.E.S.S. Collaboration 2020). Our possible detection of non-thermal emission indicates that this young star cluster could be efficiently accelerating particles (cf. Peron et al. 2024), although high-energy radiation will be difficult to detect from this cluster on account of source confusion from the bright nearby sources.

Fig. 1 Presentation of the detected targets listed in Table A.1. Each row shows a single target, labelled in the title of the first column. From left to right, the panels show (1) the VLA intensity map in mJy beam-1 with black contours at [–10, –5,5,10,20,40,80,160] σ; (2) the map of the powerlaw spectral index, α, with 10 σ intensity contour overlaid; (3) the map of the l-σ absolute uncertainty in α, with a contour showing α = –0.5; and (4) the Effelsberg intensity map at 7.639 GHz (colour scale in mJy beam-1) with VLA intensity contours overlaid and the beamsize shown for cases where it is not too large (in all cases, about 3 pixels in diameter). Also shown in the first panel are a line showing 1' scale at the bottom right, the synthetic beam FWHM at bottom left, blue dotted contours of WISE 22 pm emission, and a black arrow showing the direction of peculiar proper motion, where reliably estimated.

3.1.2 G026.1437-00.0420

This object is number 123 in the catalogue of Kobulnicky et al. (2016), with the star 2MASS J18390969-0603279 of unknown spectral type at the focal point of the IR arc. It is in projection located in the highly reddened star cluster Stephenson 2 (Stephenson 1990; Ortolani et al. 2002). It is also close to the centre of the MW bubble MWP1G026137-000346 (Simpson et al. 2012), which is probably associated with the bow shock and not a separate structure. The bow shock is clear in our observations (Figure 1, first row, left panel), with a width of ∼ 113^′′. To the south of the bow shock a point source is detected (not the driving star), which has been detected by Bruzewski et al. (2021) in their VLA follow-up of unidentified Fermi-LAT sources (source VLASS J183908.71-060359.6).

Interestingly Figure 1 shows that we predominantly detect one side of the bow shock, with the eastern side undetected although it is just as bright in IR as the northern side. This impression is artificially enhanced by the relatively strong gradient in background emission across the image from bottom-left to top-right, which we were not able to satisfactorily remove (possibly because of bright point sources to the south). So, while this apparent one-sided detection is intriguing, it should be confirmed either through a re-analysis of the dataset or from further observations.

An apparently significant peculiar proper motion is measured opposite to the direction of the bow shock, but the parallax error of the Gaia source is large and so our method is likely unreliable. For this reason we do not plot the proper-motion vector. Unfortunately, little is known about the driving star of the bow shock or its stellar wind. Our detection in radio implies that this is a true bow shock and not a bow wave because there is enhanced gas emission as well as dust emission.

The spectral index map and its uncertainty, for the bow shock and the point source are shown in the second and third panels of the first row of Figure 1, respectively. The point source has a negative spectral index, typical of active galactic nuclei, while the bow shock shows nominally a negative spectral index but with very large uncertainty. We cannot draw any conclusions about the nature of the emission (thermal or non-thermal) because the uncertainties are too large.

The spatial resolution of the Effelsberg data is comparable to the long axis of the bow-shock emission seen with the VLA and is therefore sensitive to larger-scale diffuse emission (fourth panel in first row of Figure 1). Any emission from the bow shock is diluted by diffuse emission within the large beam. There are bright point sources to the south that are also contributing significantly to the emission at the position of the bow shock in the Effelsberg map (see Figure D.1 for a larger-scale view). The spectral index in the Effelsberg map at the position of the bow shock is consistent with thermal emission.

3.1.3 G039.2486-00.0647

This candidate bow shock, entry 214 in the catalogue of Kobulnicky et al. (2016), is driven by the star 2MASS J19031363+0535230, of unknown spectral type. There is no Gaia source associated with it, and no proper motion information is available. The region is catalogued as an ultra-compact H II region by Bronfman et al. (1996), id [BNM96] 039.251–0.072. In WISE 22 µm images it appears as a small arc around the star, surrounded by a horseshoe-shaped nebula at 12 µm, presumably the compact H II region with diameter ≈45^′′. The morphology that we detect in the VLA image shown in Figure 1 (fifth row, first panel) matches closely the larger horseshoe-shaped nebula (with a width of ∼38^′′) and we can conclude that we detect the H II region. There is a hint of a double-shelled structure in the VLA image, but the spatial resolution is insufficient to clearly distinguish any putative bow shock from the H II region. This could therefore be a candidate for either a bow shock or for an IR arc around a star moving subsonically within its H II region in a dense ISM (Mackey et al. 2016).

A spectral-index map and its error are plotted in the second and third panels of Figure 1 (fifth row), respectively. the spectrum is almost flat (−0.25 ≲ α ≲ 0) where the emission is bright, indicating that the VLA emission is dominated by thermal emission, while any potential non-thermal emission is too weak to detect. The nebula is a point source for Effelsberg (fifth row, fourth panel, Figure 1).

3.1.4 G049.7065-00.1722

Entry 264 in the catalogue of Kobulnicky et al. (2016) is a semi-circular arc of mid-IR emission around the star 2MASS J19232226+1448087, of unknown spectral type. In Gaia DR3 this star has unknown distance (negative parallax) and so its significant proper motion cannot be converted into a peculiar velocity with respect to its surroundings. Like G039.2486-00.0647, the putative bow shock is surrounded by a larger ring of IR emission indicating an H II region with a somewhat cometary shape. The H II region is catalogued as [WAM82] 049.704-0.172 in Wink et al. (1982) (also possibly the radio source WSRTGP 1921+1442 mentioned in the Westerbork Synthesis Radio Telescope 327 MHz survey; Taylor et al. 1996), and in IR it is somewhat reminiscent of RCW 120 (Zavagno et al. 2007).

The radio map is shown in the sixth row of Figure 1, where it appears that we detect both the inner arc and an outer incomplete ring of emission in the VLA intensity map (first panel). The outer edge of the radio emission lines up well with the bright outer arc of emission from WISE 12 µm, suggesting that we are detecting emission from the photoionised gas and photodissociation region within the H II region, as well as the bow shock.

The second and third panels of Figure 1 (sixth row) show the spectral index map and its uncertainty, respectively. Nominally most of the region has a strongly negative spectral index, although the 1σ uncertainties are relatively large over most of the emitting region and do not allow us to draw a decisive conclusion. The uncertainties are smallest at the inner-arc, and the non-thermal spectral index derived from this region is intriguing and should be followed up with further observations at different frequencies. Figure 1 (sixth row, fourth panel) shows that we detect the region clearly with Effelsberg, but the nebula is too small to detect any substructure and the emission is most likely dominated by the H II region. The spectral index derived from Effelsberg data is flat, consistent with thermal emission.

3.1.5 BD+43^° 3654 and its bow shock EB27: BD+433654

This source has been detected at a number of radio frequencies (Benaglia et al. 2010, 2021) and was reported on extensively in our previous work (Paper I), where we also detected radio emission from the bow shock. The negative spectral index of our measurements support the previous detection of non-thermal emission from this bow shock. It is a largely extended bow shock around a massive O supergiant that has a large peculiar proper motion to the north. For completeness, the VLA intensity map, spectral index map, spectral-index uncertainty map and Effels-berg intensity map are shown in the second row of Figure 1, and we refer the reader to Paper I for further discussion. The peculiar proper motion re-calculated with Gaia DR3 data is plotted on the VLA intensity image, and is consistent with our previous calculation within the uncertainties.

3.1.6 BD+60^° 2522 and its bow shock NGC 7635

The VLA detection of the Bubble Nebula, NGC 7635, was reported on extensively in our previous work (Paper I). It is a large bow shock around a massive O giant/supergiant that has a peculiar proper motion to the north. We detected both thermal and non-thermal emission from the nebula. The radio detection is reproduced here for completeness in the fourth row of Figure 1, and we refer the reader to Paper I for further discussion.

3.2 Marginal cases

For five of our sources there is some indication of radio emission, but not sufficient that we would claim a detection of a bow shock. The VLA intensity map for each of them is shown in Figure 2 and each source is discussed in turn below. As is evident from Figure 2, in three of these cases (EB23: HIP 88652, G107.1371-00.6974 and G108.0269-00.3497) we measure faint emission that is not sufficiently above the noise level of the map to claim a detection, and further investigations with deeper surveys are strongly recommended. The other two cases show weak evidence for emission with position and morphology consistent with the IR arc, but this is overlaid on brighter diffuse emission from a surrounding H II region and so cannot be confirmed. Further higher-sensitivity observations with higher spatial resolution may be able to disentangle the two components and confirm a detection.

3.2.1 EB23: HIP 88652

This star (aka BD-14^° 4880 and HD 165319) was identified as driving a bow shock in the E-BOSS catalogue (Peri et al. 2012), and consists of a 5^′-long arc of mid-IR emission to the W of the star. It is surrounded by a H II region with various names including [LPH96] 015.128+3.310 (Lockman et al. 1996). The bow shock and H II region are also detected by MSX at 8.3 µm and in Hα nebular emission (Gvaramadze & Bomans 2008). Peri et al. (2015) detected an excess of radio emission at the position of the bow shock in NVSS 1.4 GHz emission, and Van den Eijnden et al. (2022b) reported a tentative detection in RACS data at 887.5 MHz. Considering their reported diffuse emission at a level of ≈1 mJy beam⁻¹ and the quoted beamsize, the surface brightness of their tentative detection is very similar to ours. We detect emission from the region of the bow shock, nominally above 10σ at the apex of the bow shock, but it lies at the point of constructive interference of two stripes of noise that make a faint X-shape on the image and so we cannot claim a secure detection. The image is heavily affected by imaging artefacts from the bright point source to the WNW of the star, as also reported by Van den Eijnden et al. (2022b). We measure a small peculiar proper motion of 19 km s⁻¹ in a direction consistent with the orientation of the bow shock (Table A.1). This is clearly a promising target for follow-up observations, given the consistent results at three different observing frequencies with two different telescopes.

3.2.2 G031.6770+00.1775

This bow-shock candidate, entry 159 in the catalogue of Kobulnicky et al. (2016), is actually a region of diffuse mid-IR emission, not clearly arc-shaped in WISE emission, and surrounds the star 2MASS J18483100-0102094 with unknown spectral type. In Spitzer 24 µm the higher spatial resolution reveals a faint arc of emission at a standoff distance of about 10^′′ to the N of the star. The star and bow-shock candidate are at the centre of the H II region HRDS G031.677+00.179 (Anderson et al. 2011). We clearly detect diffuse emission around the star in Figure 2, presumably free-free emission from the H II region. There is no clear arc of emission to the NNE that correlates with the IR arc, but the VLA intensity map is dominated by the larger-scale diffuse emission and we cannot claim to detect the bow shock. The peculiar proper motion we obtain is unreliable because of the large uncertainties in the Gaia parallax (Table A.1).

Fig. 2 Intensity plots of the targets listed in Table A.1 for which there is a hint of possible emission but not enough to claim a detection. The colour bar values are in mJy beam−1, with negative values shown in blue. In the bottom right, a line that corresponds to 1′ is shown, while in the bottom left the cross-section of the synthetic beam’s FWHM can be seen. Contours of the most interesting features from the WISE 22 m counterpart are overlaid with blue lines. Black contours are radio emission at levels [−10, −5, 5, 10, 20, 40, 80, 160] σ, with negative contours using a dashed line.

3.2.3 G075.1730-00.5964

This bow-shock candidate, entry 323 in Kobulnicky et al. (2016), is surrounded by the H II region [LPH96] 075.167–0.575 (Lockman et al. 1996) and driven by the central star 2MASS J20235070+3624265, of undetermined spectral type. Table A.1 shows that we could not measure a reliable peculiar proper motion for the central star because the Gaia parallax has very large uncertainties, and so we do not plot the proper motion in Figure 2. There is a 5σ contour of weak emission coinciding with the IR arc but, because the emission is largely dominated by the H II region, we cannot claim this as a detection. Indeed the background emission has many positive and negative 5σ fluctuations, implying that the region we used to extract the RMS noise level may be underestimating the true noise level of the map. We also detect significant diffuse emission that corresponds to the H II region, brighter to the N than to the S in IR images, and spatially correlated with the WISE 12 µm emission. The presented dataset does not allow us to conclude whether the star is a runaway nor if there is gaseous emission from the bow shock.

3.2.4 G107.1371-00.6974

This bow-shock candidate, entry 359 in the catalogue of Kobulnicky et al. (2016), is located to the S of the driving star (aka BD+57^° 2606 and HD 215806), reported as an Oe star by Chick et al. (2020). It is located on the periphery of the blister H II region [ABB2014] WISE G107.156-00.988 (Anderson et al. 2014) driven by the open cluster NGC 7380, with gas expansion in the direction opposite to the orientation of the bow shock. As such it could be an example of the ‘interstellar weather vanes’ found by Povich et al. (2008) in M17 and RCW 49, where the orientation of the bow shock traces not the motion of the star, but a bulk flow of gas escaping from a high-pressure region. Supporting this possibility, we detect no significant peculiar proper motion (≈5 km s⁻¹) for the driving star (Table A.1). We detect weak emission in the region of the IR arc, shown in Figure 2, but again the emission is located at the superposition of noise stripes across the image and so we cannot claim a detection.

3.2.5 G108.0269-00.3497

Entry 360 in the catalogue of Kobulnicky et al. (2016), this bow shock candidate is located to the SW of the driving star BD+58^° 2492 (aka HD 216411), a B0.7Ia supergiant (Negueruela et al. 2024). Figure 2 shows that we detect faint emission spatially coincident with the bow shock, but it is superimposed on a stripe of noise from a nearby bright point source that enhances the signal artificially above the 5σ threshold. We also measure a large and well-constrained peculiar proper motion of the star of ≈85 km s⁻¹ in a direction appropriate for the orientation of the bow shock.

3.3 Non-detections

Of the 24 targeted sources, 13 were not detected by the VLA and these fields are shown in the Appendix in Figure C.1. Some are clearly affected by bright sources outside the field of view, whereas others simply have no detection down to the RMS noise level of the map, either because we are not sensitive to the spatial scales involved (which should not occur because of our selection criteria) or because the gaseous emission is much weaker than expected based on the bright IR emission from dust. Discussion of each individual target is in the Appendix.

Fig. 3 Comparison of observations of the bow shock of ζ Oph with the VLA (left), Effelsberg (centre), and the mid-IR image from Spitzer Space Telescope MIPS 24 µm. All three images are plotted in units of megaJansky per steradian. The VLA and Effelsberg beams are plotted at the lower-left of the respective panels. For the VLA and Effelsberg plots, IR emission contours at levels 70 and 100 MJy sr−1 are overplotted in dotted blue lines, and for Effelsberg and Spitzer plots the VLA intensity contours are overlaid in black, showing only the detected point sources.

3.4 Limits on electron density from radio emission and upper limits

Following the method outlined in Section 2.5, we can obtain upper limits on the electron density, n_e, for the detections and upper limits. The results of this estimate for our 24 sources are shown in Table 2, assuming an observing frequency ν = 6 GHz. For the three well-studied bow shocks in our sample, around the stars ζ Oph, BD+43^° 3654 and BD+60^° 2522, we can compare our estimated n_e to literature measurements.

ζ Oph: Gvaramadze et al. (2012) estimated n_e ≈ 3 3 cm⁻ for the ISM upstream of the bow shock, and our upper limit for the shocked layer (n_e ≲ 39.7 cm⁻³) corresponds to a compression factor ≤13, although it should be noted that for this very extended bow shock the VLA may be missing flux from large scales, weakening this constraint. Green et al. (2022) simulated the bow shock of ζ Oph and made synthetic 6 GHz bremsstrahlung emission maps, predicting a peak brightness of I_ν ≈ 0.035 MJy sr⁻¹, corresponding to 0.10 mJy beam⁻¹ for the VLA beam. This is twice the 3σ upper limit from Table 2, so either the ISM density considered in Green et al. (2022) was significantly too large, or we are losing a lot of flux in the VLA image. To test this, ζ Oph is the only bow shock in our sample that should be spatially resolved with Effelsberg, and so we can perform the same analysis with Effelsberg data. Using maps at 4.951 GHz and 6.599 GHz, with RMS noise levels of 3.8 mJy beam⁻¹ and 3.5 mJy beam⁻¹, beamsizes of 141^′′ and 110^′′, respectively, we obtain 3σ upper limits of n_e(ν = 4.951 GHz) ≲ 44.7 cm⁻³ and n_e(ν = 6.599 GHz) ≲ 54.6 cm⁻³. The constraints are similar to (but higher than) the VLA limit, and are significantly more robust because there is no loss of large-scale flux. A comparison of the VLA and Effelsberg data is shown in Figure 3, where we plot VLA, Effelsberg, and Spitzer MIPS 24 µm data all in units of MJy sr⁻¹. It is notable that the mid-IR emission is at least 10⁴× more luminous than the radio emission.

BD+60^° 2522: Our estimate of n_e ≲ 250 cm⁻³ is larger than estimates of n_e ≈ 100–4200 cm⁻³ from optical nebular spectroscopy (Esteban et al. 2016). The similarity of these numbers nevertheless implies that thermal emission should be contributing measurably to the total emission at 6 GHz and should start to dominate at higher frequency. It should be noted that the crude assumptions of constant n_e and the simple geometry likely introduce uncertainties of at least a factor of roughly two in the estimates of n_e, so we cannot draw any clear conclusions about the relative contributions of thermal and non-thermal emission to the flux.

BD+43^° 3654: The gas density in the pre-shock ISM was estimated by Comerón & Pasquali (2007) to be n_e ≈ 6 cm⁻³, and so our upper limit of n_e ≲ 69.5 cm⁻³ implies a compression factor of ≲ 12. We estimate the space velocity of the star to be v_⋆ ≈ 45 km s⁻¹ (Table A.1), which would give a Mach number of ≈4 and compression factor of ≈16 for an isothermal shock in the photoionised ISM. As noted for BD+60^° 2522, these crude estimates have rather large and poorly quantified uncertainties, but imply that thermal emission should be contributing much of the emission at and above 6 GHz.

In both cases of BD+43^° 3654 and BD+60^° 2522, further observations at lower and higher frequencies should be performed to measure the transition from a spectrum dominated by synchrotron emission (at low ν) to one dominated by thermal emission (at high ν). The lower frequencies have been measured by Benaglia et al. (2021) for BD+43^° 3654, but higher frequency measurements would be valuable to measure the thermal contribution. A thermal contribution to the flux from the bow shock would reduce the large energy-requirements for high-energy particles noted in Paper I. They assumed that all of the emission was synchrotron and found that it was difficult to inject enough energy into the region to explain the emission, particularly for BD+60^° 2522 where both newly accelerated particles and compressed interstellar cosmic-ray electrons struggled to explain the level of emission.

It should be noted that the H II regions around BD+60^° 2522 and BD+43^° 3654 should contribute to the emission measure in addition to thermal and non-thermal emission from the bow shock, although the VLA observations may not be sensitive to any large-scale uniform emission. Indeed Fig. 1 shows no large-scale diffuse emission for either of these sources, indicating that The VLA has resolved out the H II region and it does not contribute to the estimated emission measure.

For the other sources we obtain the upper limits n_e ∈ [12– 125] cm⁻³, which is appropriate for shock-compressed gas in the diffuse ISM, considering compression factors ≈10 and bearing in mind that bow shocks are more likely to be detected for stars moving in higher density environments. The exception is G018.2660-00.2988, which is located at the edge of a young star cluster in a massive star-forming region, and so the measured n_e ≲ 411 cm⁻³ is much larger than for the stars in the diffuse ISM. Alternatively it could be an indication of strong non-thermal emission from the bow shock. Considering also the size of the bow shock,_RR₀, we estimate upper limits on the emission measure, $EM={\displaystyle \int n_e^2}d\ell$, in the range 57 to 5200 cm⁻⁶ pc, with typical values in the range 100 to 500 cm⁻⁶ pc.

We can compare our upper limits with those derived from the measured R₀ and v_⋆, together with estimated mass-loss rates and wind velocities in Table B.1. Re-arranging Equation (3) gives the unshocked ISM density, and the associated electron number densities, n_e,0, are quoted in Table B.1. In most cases n_e,0 < n_e with the notable exception of EB23: HIP 88652, for which we obtain a large n_e,0 = 120.5 cm⁻³ on account of the large mass-loss rate $\left(\dot{M}=5.3\times {10}^{-6}{\text{M}}_{\odot }{\text{yr}}^{-1}\right)$ suggested by Gvaramadze & Bomans (2008). O9.5I supergiants have a large range of measured $\dot{M}$, however, and if we use the lowest value from Mokiem et al. (2007) of $\dot{M}=7\times {10}^{-7}{\text{M}}_{\odot }{\text{yr}}^{-1}$ then the ISM density is reduced to n_e,0 = 15.9 cm⁻³. This is still quite large, but no longer in conflict with the radio-derived upper limit on n_e, and consistent with the density estimate of the surrounding H II region by Gvaramadze & Bomans (2008).

4 Discussion

This work constitutes the largest targeted radio survey of bow shocks to date, and the detection rate of 25% for clear detections, and 46% when also including possible detections, demonstrates the success of our targeted approach. The only research on a similar scale was by Van den Eijnden et al. (2022b), who used the publicly released Rapid ASKAP Continuum Survey (McConnell et al. 2020) data to search for radio emission from 50 IR-selected bow shocks, finding 3 clear detections and 3 tentative detections, a success rate of 6% (12%) for clear (tentative) detections. The lower success rate of Van den Eijnden et al. (2022b) may be attributed to the fact that their dataset was a very wide-field survey rather than targeted deep radio observations. Using a less sensitive and lower resolution dataset (the NVSS), Peri et al. (2015) could find associations of radio emission with a number of bow shocks, some of which were subsequently confirmed by Van den Eijnden et al. (2022b). This highlights the necessity of deep radio observations to detect emission from bow shocks.

Despite this success, of our six detections only 3 are clearly bow shocks whereas the other three appear to have strong contaminating emission from the surrounding H II region, which limits the degree to which we can interpret our results. Similarly, for the tentative and possible detections that require follow-up confirmation, two of the five sources have significant contamination from surrounding diffuse emission. Only approximately half of our firm and possible detections appear as a clean bow-shock structure in radio emission, whereas all are clearly arc-shaped bow shocks in mid-IR imaging. Similar to the survey of Van den Eijnden et al. (2022b), we find tantalising evidence that slightly deeper observations would give clear detections of even more bow shocks. Our results should be used to plan future observations, including consideration of the locations of bright background sources (particularly for EB23: HIP 88652, G107.1371-00.6974 and G108.0269-00.3497). Assuming these are background active galactic nuclei with steep spectra, they could be even more problematic for imaging at lower frequencies where non-thermal emission from the bow-shock candidates would be stronger.

If we do not detect gas emission from a bow-shock candidate, it could indicate that the IR emission comes from pile-up of dust grains that are dynamically decoupled from the ISM by the radiation pressure of the star – a dust wave (Henney & Arthur 2019) in which the overdensity of dust is not coincident with a gas overdensity. Alternatively, it could merely indicate insufficient sensitivity to the gaseous emission in radio. For the non-detections, the 3σ upper limits on the emission can be converted to upper limits on n_e in the shocked layer (with some symmetry assumptions) in the range ≲[12–125] cm⁻³, with a few tens being typical. For an adiabatic shock in a low-density environment (n_e ∼ [0.1–1] cm⁻³) this is not very constraining and deeper observations are required (or analysis techniques that better remove the striping artefacts from bright sources within and outside the field of view). For an isothermal shock, typical of denser environments with short cooling timescales, we expect a compression factor of approximately ten for runaway stars with the typical space velocity v_⋆ ≈ 30 km s⁻¹ for a runaway star, meaning that our upper limits constrain the ISM density to n_e ≲ [1.2–12] cm⁻³. The upper limits are therefore not unreasonable for any of the undetected sources and we cannot conclude that gas-dust decoupling is taking place.

For ζ Oph our upper limits on the emission with the VLA and Effelsberg are already close to being in conflict with predictions from magnetohydrodynamic simulations (Green et al. 2022), and so a modest increase in sensitivity of a factor of three should be sufficient to either detect the bow shock or rule out much of the parameter space in wind power and ISM density considered by Green et al. (2022). These results underline how much more luminous bow shocks are in mid-IR dust emission than in tracers of gas emission, as noted by predictions from simulations (Meyer et al. 2014; Mackey et al. 2016; Acreman et al. 2016). Our upper limit for radio emission from ζ Oph shows that the mid-IR emission is at least 10⁴× more luminous than at 6 GHz.

Nearly all of our detections and possible detections show strongly negative spectral index, α, in the VLA analysis, but all have flat spectra consistent with optically thin bremsstrahlung at the source position in the Effelsberg data. Furthermore, the uncertainties on the VLA-derived α are large for most sources, except for BD+43^° 3654, BD+60^° 2522 and the very bright sources G018.2660-00.2988 and G039.2486-00.0647, both of which show α consistent (at least marginally) with thermal emission from the surrounding H II region. G049.7065-00.1722 shows good evidence for non-thermal emission from the bow shock, which should be followed up with further observations, although for this object our lack of knowledge about the driving star limits what can be learned from the nebula.

Green (2022) noted that interferometric observations can give spuriously negative spectral indices for faint and diffuse sources and so, given the relatively large uncertainties in the VLA-derived spectral indices for the fainter sources, we do not consider them reliable enough to make any claims about the emission mechanism. We attempted combining the single-dish and interferometric datasets using the CASA task ‘feather’ (Rau et al. 2019), as was done in Paper I for BD+43^° 3654 and BD+60^° 2522, but the outcome was less convincing for these other sources because the resulting spectral index obtained was strongly dependent on the choice of weighting. Attempting to plot the spectrum in the 4–12 GHz range using the VLA and Effelsberg sub-bands was not attempted for these sources. With the current dataset we cannot claim to have convincingly discovered any more non-thermal-emitting bow shocks, but nor can we exclude that the bow shocks have non-thermal emission. Further observations, ideally at similar sensitivity and spatial resolution but at lower (or higher) frequency, are required to determine a robust spectral index for the faint diffuse emission. It may be possible to detect polarised emission, at least for the bright bow shocks and nebulae in our sample, which would be a clear signature of synchrotron emission, independent of the spectral-index measurement (e.g. Benaglia et al. 2021).

Kobulnicky & Chick (2022) also measured the 2D space velocity of massive stars with bow shocks. Of our 24 targets, only BD+43^° 3654 is contained in their Table 1, for which our peculiar velocity is almost identical. Their Table 3 contains numbered nebulae but it is not clear how to relate these numbers to our sources and so we could not make a comparison. Similarly, Carretero-Castrillo et al. (2023) and Carretero-Castrillo et al. (2025) studied kinematics of stars driving bow shocks. There are four common targets with our survey, EB31: HD 57682, EB23: HIP 88652, BD+43^° 3654 and BD+60^° 2522, and for all of these we find consistent peculiar velocities within 1σ uncertainties by comparing with their $V_{\text{PEC}}^{2\text{D}}$ values.

Similar to Kobulnicky & Chick (2022) and Carretero-Castrillo et al. (2025), we find the majority of our stars are moving slowly through their environment. Out of 18 stars with well-determined v_⋆, we find only 3 with v_⋆ ≥ 30 km s⁻¹, 6 more with 15 km s⁻¹ ≤ v_⋆ < 30 km s⁻¹, and 9 with v_⋆ < 15 km s⁻¹. This means we have only 3 stars satisfying the classical definition of a runaway star (v_⋆ ≥ 30 km s⁻¹), with 15 being classified as walkaway stars (v_⋆ < 30 km s⁻¹), at least when considering only the motion in the plane of the sky.

5 Conclusions

We have carried out a survey of 24 IR-selected bow shocks in the radio C and X bands with the VLA, following up on ten of these targets with Effelsberg single-dish C band observations. We clearly detected arc-shaped emission from six of the 24 targets. In three of these sources, the interpretation is somewhat complicated by the surrounding diffuse radio emission (presumably) from the H II region that the bow shocks are embedded in. For five more sources, we observed evidence of emission that could be a low-significance detection of the bow shock, but noise in the map and imaging artefacts from bright point sources prevents a convincing detection. In two of these cases there is significant emission from a surrounding H II region that further obscures the picture.

Excluding the bright sources BD+43^° 3654 and BD+60^° 2522 with synchrotron spectra reported on in Paper I, two of the bright sources (G018.2660-00.2988 and G039.2486-00.0647) have a relatively flat spectral index that may be consistent with thermal emission. G049.7065-00.1722 may be a non-thermal emitter, and G026.1437-00.0420 has such a large uncertainty on the spectral index map that we cannot draw any conclusions. Therefore, we do not add to the three confirmed non-thermal-emitting bow shocks that have been reported in the literature, but we identify strong candidates.

Interpreting the entirety of the radio emission as thermal bremsstrahlung emission in the optically thin limit, we estimated upper limits on the gas density within the shocked-ISM layer of the bow shock, finding that all upper limits are consistent with expectations for shock-compressed bow shocks of stars moving through the diffuse ISM in the Galactic Plane. For the case of BD+60^° 2522, the well-constrained gas density obtained from optical spectral lines is lower than our upper limits, which is consistent with the conclusion of Paper I that a significant part of the emission is non-thermal. For this case it appears that both thermal and non-thermal emission contribute at a similar level to the emission at ≈6 GHz, so broader spectral coverage should be able to constrain both components.

Our attempts to detect the bow shock of ζ Oph were unsuccessful both with the VLA and Effelsberg observations. Comparison of the derived upper limits with simulation results suggests that improved sensitivity by a factor of three with Effelsberg should yield a detection or would very strongly constrain possible models for the bow shock.

Follow-up observations of our confirmed and possible detections at both higher and lower frequencies are strongly encouraged. These would allow the contributions of thermal and non-thermal emission from the sources to be disentangled and give useful constraints on the gas density and the population of relativistic electrons.

Data availability

FITS files used to generate Figures 1, 2, 3, C.1 and D.1 are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/704/A268

Acknowledgements

JM acknowledges support from a Royal Society-Science Foundation Ireland University Research Fellowship. NC gratefully acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG) – CA 2551/1-2. PJH acknowledges support from Universidad Nacional Autónoma de México Postdoctoral Program (POSDOC). YG was supported by the National Natural Science Foundation of China (Grant No. 12427901), the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDB0800301), and the Ministry of Science and Technology of China under the National Key R&D Program (Grant No. 2023YFA1608200). This publication results from research conducted with the financial support of Taighde Éireann – Research Ireland under Grant numbers 14/RS-URF/3219, 17/RS-EA/3468, 20/RS-URF-R/3712, 22/RS-EA/3810, IRCLA\2017\83. M.R.R. is a Jansky Fellow of the NSF National Radio Astronomy Observatory, USA. This work is based on observations with the National Science Foundation’s (NSF) Karl G. Jansky Very Large Array (VLA), operated by the National Radio Astronomy Observatory (Project ID: 19B-105; PI: M. Moutzouri). The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This research made use of the Common Astronomy Software Applications (CASA, McMullin et al. 2007). This work is based on observations with the 100 m telescope of the MPIfR (Max-Planck-Institut für Radioastronomie) at Effelsberg (Project ID: 86-20; PI: M. Moutzouri). We thank the operations team at the Effelsberg 100 m telescope for their assistance with our observations. This research has made use of the SIM-BAD database, operated at CDS, Strasbourg, France (Wenger et al. 2000). This research has made use of the Astrophysics Data System, funded by NASA under Cooperative Agreement 80NSSC25M7105. This research made use of the python packages APLpy, (Robitaille & Bressert 2012), Astropy (Astropy Collaboration 2018), Numpy (Harris et al. 2020), matplotlib (Hunter 2007).

References

Appendix A Targets

Table A.1 shows a list of the 24 bow-shock candidates targeted with VLA observations, including literature information on the driving star’s name/ID, spectral type, parallax and peculiar proper motion, where these could be reliably estimated, as discussed in Section 2.

Appendix B Estimated mass-loss rates and wind velocities for the bow-shock driving stars

For 17 of the 24 targets we found spectral types in the literature and obtained mass-loss rates, $\dot{M}$, and wind velocities, v_∞, using the method described in Section 2.5. These values are listed in Table B.1, together with the derived ISM electron density, n_e,0, from Equation 3. Note that we do not consider a possible radial velocity component to v_⋆, and so v_⋆ may be underestimated, leading to an overestimate of n_e,0. We also assume the bow shock is oriented perpendicular to our line of sight, which may not be true and introduces some uncertainty in the observational estimate of R₀. Some stars have v_⋆ < c_s such that a bow shock will not form because stellar motion is not supersonic in the surrounding ISM. In this case n_e,0 is not quoted because Equation 3 does not apply.

Appendix C Sources not detected with the VLA

Looking at each of the non-detected targets in turn, and comparing with the literature and the source catalogues, Peri et al. (2015) and Kobulnicky et al. (2016), we can make the following observations on each source:

1. G119.4436-00.9208: Entry 364 of the candidate bow shocks in Kobulnicky et al. (2016) is driven by the O9V star BD +60^°39 and shows an IR arc to the SW of the star. Figure C.1 shows that we do detect a moderate peculiar proper motion of the star (indicating a walkaway star) in the direction suggested by the bow shock’s orientation, but we do not detect any radio emission coinciding with the IR emission. Again this may be due to a lack of sensitivity, or the IR emission could arise from a bow wave that is not accompanied by an increase in the gas density.

2. EB01: HIP 2036: Entry 5 of the E-BOSS catalogue of Peri et al. (2012) is a bow shock candidate around the B0.2V star HD 2083 (aka BD+71^° 16 and HIP 2036) with a faint IR arc to the north of the star. The emission is undetected in our VLA observations (Figure C.1), and nor do we detect a significant peculiar proper motion for the star.

3. EB29: SER6: This bow shock candidate from group 7 of the E-BOSS II catalogue (Peri et al. 2015) is a bright mid-IR arc of emission around the star TYC 3677-890-1 (aka 2MASS J01112094+5733282), which does not have a spectral type reported in Simbad. Figure C.1 shows that we detect no diffuse radio emission down to the noise level of the map (which is very noisy compared with most of the other targets because of the bright source to the south). There is a small peculiar proper motion for the star (indicating a walkaway star), but in a direction opposite to what would be expected if this source were a bow shock.

4. G133.1567+00.0432: Entry 367 of the candidate bow shocks in Kobulnicky et al. (2016) is driven by the late O star LS I +60 226. We do not detect the bow shock in radio (Figure C.1), nor do we measure a significant peculiar proper motion for the star (≈7 km s⁻¹), although in the IR images it appears to be an excellent candidate.

5. G134.3552+00.8182: Entry 368 of the candidate bow shocks in Kobulnicky et al. (2016) is driven by the late O-type variable star V^⋆ KM Cas, and appears to be within a large H II region, with the IR arc in the direction of the closest bright nebulosity. It is not detected by the VLA survey (Figure C.1). As such, it could also be an IR arc tracing the edge of a wind bubble around a star moving subsonically with respect to its surroundings (Mackey et al. 2016) and we would not expect to see an excess emission in free-free emission. Consistent with this we do not detect large proper motion for the star (≈12 km s⁻¹, indicating a walkaway star).

6. G137.4203+01.2792: Entry 369 of the candidate bow shocks in Kobulnicky et al. (2016) is driven by the late O7Vz star BD +60^°586 (Sota et al. 2011), and again appears to be located within a large H II region and with the IR arc in the direction of the nearest bright nebular emission to the east. There is a vague hint of some emission but it is well within the noise level of the map and the WISE contour in Figure C.1 leads the eye somewhat to see some emission. It may be worth following up with deeper observations, although the lack of any significant peculiar proper motion (< 5 km s⁻¹) argues against a bow-shock interpretation of the IR emission.

7. G151.0318-00.8271: Entry 371 of the candidate bow shocks in Kobulnicky et al. (2016) is driven by the O5V((f)) star TYC 3339-851-1 (Roman-Lopes & Roman-Lopes 2019) (aka 2MASS J04055303+5106580) and surrounded by mid-IR nebular emission brightest to the SW. Figure C.1 shows that we do not detect any radio emission from the bow shock, but we do measure a small peculiar proper motion for the star appropriate for the bow-shock interpretation (indicating a walkaway star).

8. EB09: HIP 26397: The B0.5V star (aka HD 37032 and HIP 26397) drives bow shock 63 from the E-BOSS catalogue of Peri et al. (2012), which is undetected in our VLA observations (Figure C.1). Interestingly, the IR nebulosity in the cutout image of Peri et al. (2012) looks more similar to an Archimedean spiral or a ring nebula than an arc. We do detect a small (and uncertain) peculiar proper motion for the central star of ≈7 km s⁻¹ that is consistent with the orientation of the brightest part of the IR emission.

9. EB12: HIP 31766: this star (aka HD 47432 and HIP 31766) is an O9.7I supergiant (Sota et al. 2011) and listed as a pulsating variable in Simbad. It has what appears to be a double-arc structure in mid-IR emission, pointing to the SW, and is entry 78 in the E-BOSS catalogue of Peri et al. (2012). It is undetected in our VLA observations (Figure C.1), although there is a hemisphere of low-level positive emission that correlates with the position of the IR arc, and so this may be worth following up with deeper observations. We also do not detect any significant peculiar proper motion of the star (≈4 km s⁻¹), and so it is possible that this is a dusty ejecta nebula and not a bow shock.

10. EB31: HD 57682⁷: this star is surrounded by an IR arc detected as entry 103 in the E-BOSS II catalogue of Peri et al. (2015). The star (aka HD 57682) is an O9.7IV sub-giant (Sota et al. 2014) and we do not detect any radio emission. This could be due to lack of sensitivity or it could be an indication that the IR arc is a bow wave driven by radiation pressure, which leads to an overdensity in dust but not in gas (van Buren & McCray 1988; Henney & Arthur 2019). We measure a large peculiar proper motion for the star of ≈90 km s⁻¹ in the direction consistent with the orientation of the IR arc.

11. ζ Oph: the well-known bow shock around the closest O star to Earth, ζ Oph, has been detected in optical line emission (Gull & Sofia 1979), IR emission (van Buren & McCray 1988) and diffuse X-rays from the shocked stellar wind (Toalá et al. 2016; Green et al. 2022). It is, however, a very large bow shock in terms of angular scale and, possibly because of its large angular size, we do not detect it with the VLA (Figs. 3 and C.1). As discussed in section 2.5, the non-detection with Effelsberg gives a slightly weaker constraint statistically, but is significantly more robust because large-scale flux is not lost. Our proper motion calculation is in good agreement with other values from the literature (e.g. Gvaramadze et al. 2012; Green et al. 2022).

12. G026.5272+00.3808: In this field of view we detect bright radio emission from the H II region IRAS 18355-0532 (Bronfman et al. 1996), but the arc-shaped bow shock (candidate 121 in Kobulnicky et al. 2016, labelled G026.5272+00.3808) in the direction WSW of the star is not detected (Figure C.1). The star driving the bow shock appears to be 2MASS J18382147-0531233. The proper motion that we estimate from Gaia DR3 is unreliable because the parallax is negative.

13. G078.2889+00.7829: this bow shock candidate, entry G078.2889+00.7829 in Kobulnicky et al. (2016), surrounds the O8V(n)((f)) star LS II 39 53 (Maíz Apellániz et al. 2016), which is located within a large H II region, and we do not detect any radio emission with the VLA (Figure C.1). It may be an IR arc from a slowly moving star within a H II region that does not produce a bow shock, in which case the emissivity of the arc is enhanced in IR but not in tracers of gas density such as radio free-free emission (Mackey et al. 2016). Consistent with this interpretation, we do not detect a large peculiar stellar motion (≈9 km s⁻¹) for the star, and the direction is opposite to what would be expected if the IR arc were a bow shock.

Fig. C.1 Intensity plots of the non-detected targets listed in Table A.1. We use a saturated colour scale extending to negative values to show the noise level in the maps. The colour bar values are in milliJansky per beam. In the bottom right, a line that corresponds to 1′ is shown, while in the bottom left, the cross-section of the synthetic beam’s FWHM can be seen. Contours of the most interesting features from the WISE 22 μm counterpart are overlaid with blue lines.

Appendix D Effelsberg maps of larger-scale emission

Figure D.1 presents the Effelsberg 6.615 GHz radio continuum images for the nine selected sources (all except ζ Oph which is shown in Figure 3). While continuum maps were created at seven frequency bands (see Sect. 2.3), we show only the 6.615 GHz images here as a representative example to illustrate the data quality and morphological features. The naming convention of the sources is given in Table D.1

Fig. D.1 Effelsberg 6.615 GHz continuum images of the nine sources. The source name is labelled in the top left corner of each panel.

All Tables

All Figures

Fig. 1 Presentation of the detected targets listed in Table A.1. Each row shows a single target, labelled in the title of the first column. From left to right, the panels show (1) the VLA intensity map in mJy beam-1 with black contours at [–10, –5,5,10,20,40,80,160] σ; (2) the map of the powerlaw spectral index, α, with 10 σ intensity contour overlaid; (3) the map of the l-σ absolute uncertainty in α, with a contour showing α = –0.5; and (4) the Effelsberg intensity map at 7.639 GHz (colour scale in mJy beam-1) with VLA intensity contours overlaid and the beamsize shown for cases where it is not too large (in all cases, about 3 pixels in diameter). Also shown in the first panel are a line showing 1' scale at the bottom right, the synthetic beam FWHM at bottom left, blue dotted contours of WISE 22 pm emission, and a black arrow showing the direction of peculiar proper motion, where reliably estimated.

In the text

Fig. 2 Intensity plots of the targets listed in Table A.1 for which there is a hint of possible emission but not enough to claim a detection. The colour bar values are in mJy beam−1, with negative values shown in blue. In the bottom right, a line that corresponds to 1′ is shown, while in the bottom left the cross-section of the synthetic beam’s FWHM can be seen. Contours of the most interesting features from the WISE 22 m counterpart are overlaid with blue lines. Black contours are radio emission at levels [−10, −5, 5, 10, 20, 40, 80, 160] σ, with negative contours using a dashed line.

In the text

Fig. 3 Comparison of observations of the bow shock of ζ Oph with the VLA (left), Effelsberg (centre), and the mid-IR image from Spitzer Space Telescope MIPS 24 µm. All three images are plotted in units of megaJansky per steradian. The VLA and Effelsberg beams are plotted at the lower-left of the respective panels. For the VLA and Effelsberg plots, IR emission contours at levels 70 and 100 MJy sr−1 are overplotted in dotted blue lines, and for Effelsberg and Spitzer plots the VLA intensity contours are overlaid in black, showing only the detected point sources.

In the text

Fig. C.1 Intensity plots of the non-detected targets listed in Table A.1. We use a saturated colour scale extending to negative values to show the noise level in the maps. The colour bar values are in milliJansky per beam. In the bottom right, a line that corresponds to 1′ is shown, while in the bottom left, the cross-section of the synthetic beam’s FWHM can be seen. Contours of the most interesting features from the WISE 22 μm counterpart are overlaid with blue lines.

In the text

	Fig. D.1 Effelsberg 6.615 GHz continuum images of the nine sources. The source name is labelled in the top left corner of each panel.
In the text